! compareglo2.f90
! f90 main program written on 13.01.00 by Tim Mitchell
! last modification on 25.10.01
! pgf90 -o ./../goglo/compareglo2 filenames.f90 initialmod.f90 glofiles.f90 sortmod.f90 
!    ./../goglo/compareglo2.f90
! purpose is to permit us to avoid multiple saves to .glo and truncations to zero   

program CompareGlo2

use FileNames
use InitialMod
use GloFiles
use SortMod

implicit none

real, pointer, dimension (:,:)		:: GloLoaded
real, pointer, dimension (:)		:: GloSaved, ThisGlo, BaseGlo, NewTot, NewEn, GloMin, GloMax, OpVec
real, pointer, dimension (:)		:: PertRaw, PertTot, BaseTot

integer, pointer, dimension (:) 	:: WorkRegSizes, Indices, NewRegSizes
integer, pointer, dimension (:,:) 	:: WorkMapIDLReg, NewMapIDLReg

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

real, allocatable, dimension (:,:)	:: UseClass

integer, allocatable, dimension (:)	:: UseGlo, FreqClass

real, parameter :: MissVal = -999.0

real :: RunTot, RunEn, RunSqTot, RunMin, RunMax
real :: Aye, Bee, Are, Num, Den
real :: SumX, SumY, SumXX, SumYY, SumXY, En
real :: OpTot, OpSqTot, OpAbsTot, OpEn, OpMean, OpStDev, OpDiff
real :: RegSqTot, RegTot, RegEn, RegMean, RegStDev, RegRMSE
real :: CutOffVal, TestVal, ChosenKay, RealConstant, AvMag, XMedianReal
real :: MinVal, MaxVal, Median, Mode
real :: ClassValue, ClassBound, LastBound, RangeMin, RangeMax, RangeGap
real :: ValueToChange, ValueToLeave

integer :: WorkGrid, WorkLongN, WorkLatN, WorkDataN
integer :: WorkRegN, WorkGloN, NewRegN
integer :: SelectGloN, ChosenGlo
integer :: AllocStat, ReadStatus
integer :: XReg, XGlo, XGloFree, XGloSave, XSelect, XClass, XBound, XLat, XLong
integer :: QAnom, QMain, QSampPop, QConvert, QWhichOp, QOutput, QMinMax, QClassMethod
integer :: CutOffSign, CutOffType, RegKept, RegCutOff
integer :: ChosenExp, Operator, IntConstant, XMedianInt, OpMissN, ExitNow
integer :: ClassN, ValidN

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

!*******************************************************************************
! initialise

open (99,file="./../../../scratch/log-compareglo2.dat",status="replace",action="write")

print*
print*, "  > ***** CompareGlo *****"
print*

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

call GrabGlo

Blank = ""
QMain = -1

do
  if      (QMain.EQ.-1) then
  else if (QMain.EQ.1) then
    SelectGloN = 2
    call SelectGloOp
    call TwoFromOne
  else if (QMain.EQ.2) then
    call WhichOp
    call GrabGloN
    call SelectGloOp

    if (QWhichOp.EQ.1) call CalcMean
    if (QWhichOp.EQ.2) call CalcMedian
    if (QWhichOp.EQ.3) print*, "  > ##### CompareGlo: Not set up to calc mode. #####"
    if (QWhichOp.EQ.4) call CalcStDev
    if (QWhichOp.EQ.5) call CalcRange
  else if (QMain.EQ.5) then
    SelectGloN = 2
    call SelectGloOp
    call CalcLSR
  else if (QMain.EQ.6) then
    SelectGloN = 3
    call SelectGloOp
    call PosABC
  else if (QMain.EQ.7) then
    SelectGloN = 2
    call SelectGloOp
    call MaskAB
  else if (QMain.EQ.8) then
    SelectGloN = 1
    call SelectGloOp
    call OperateAConstant
  else if (QMain.EQ.9) then
    SelectGloN = 2
    call SelectGloOp
    call OperateAB
  else if (QMain.EQ.10) then
    call GrabGloN
    call SelectGloOp
    call ClassifyABox
  else if (QMain.EQ.11) then
    call GrabGloN
    call SelectGloOp
    call ClassifyAGroup
  else if (QMain.EQ.12) then
    SelectGloN = 2
    call SelectGloOp
    call CompareAB
  else if (QMain.EQ.13) then
    SelectGloN = 1
    call SelectGloOp
    call RegionaliseA
  else if (QMain.EQ.14) then
    SelectGloN = 1
    call SelectGloOp
    call RegPerCentA
  else if (QMain.EQ.15) then
    SelectGloN = 1
    call SelectGloOp
    call TransformVal
  else if (QMain.EQ.16) then
    SelectGloN = 2
    call SelectGloOp
    call CompMasked
  else if (QMain.EQ.20) then
    call SaveToGlo
  else if (QMain.EQ.30) then
    call PrintToScreen
  else if (QMain.EQ.31) then
    call PrintToScreenExtra
  else if (QMain.EQ.32) then
    SelectGloN = 1
    call SelectGloOp
    call WriteRegLog
  else if (QMain.EQ.33) then
    SelectGloN = 1
    call SelectGloOp
    call WriteRegScreen
  else
    print*, "  > Select the comparison to make: "
    print*, "  >   1. A minus B"
    print*, "  >   2. Summary operations on A,B,..."
    print*, "  >   5. LSR between A and B"
    print*, "  >   6. Position of A relative to B, C"
    print*, "  >   7. Mask A using B"
    print*, "  >   8. Operate on A with k"
    print*, "  >   9. Operate on A with B"
    print*, "  >  10. Classify A: assign value to each box"
    print*, "  >  11. Classify A: group values and report"
    print*, "  >  12. Compare A and B"
    print*, "  >  13. Convert A to new regions (save/print)"
    print*, "  >  14. As 13. For percentage units"
    print*, "  >  15. Convert one value to another"
    print*, "  >  16. Compare masked columns"
    print*, "  >  20. Save A to .glo"
    print*, "  >  30. Print stats for all to screen"
    print*, "  >  31. Print extra stats"
    print*, "  >  32. Write a column to log"
    print*, "  >  33. Print A's regional values to screen"
    print*, "  >  99. Exit."
  end if
  
  if (allocated(UseGlo)) deallocate (UseGlo, stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: Main: Deallocation failure #####"
  
  print*
  print*, "  > Main menu: select option: "
  do
	read (*,*,iostat=ReadStatus), QMain
	if (ReadStatus.LE.0) exit
  end do
  
  if (QMain.EQ.99) exit
end do

deallocate (GloLoaded)
deallocate (WorkMapIDLReg,WorkRegSizes,WorkRegNames)

print*

close (99)

contains

!*******************************************************************************
! choose operation on A,B,...

subroutine WhichOp

print*, "  > Select the operation to perform:"
print*, "  > 1=mean, 2=median, 3=mode, 4=stdev, 5=range"
do
	read (*,*,iostat=ReadStatus), QWhichOp
	if (ReadStatus.LE.0.AND.QWhichOp.GE.1.AND.QWhichOp.LE.5) exit
end do  

end subroutine WhichOp

!*******************************************************************************
! load .glo

subroutine GrabGlo

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

XGloFree = WorkGloN + 1
WorkGloN = WorkGloN + 30

allocate (GloLoaded(WorkGloN,WorkRegN), &
	  ThisGlo  (WorkRegN),          stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: GrabGlo: Allocation failure #####"
GloLoaded = MissVal

do XGlo = 1, (WorkGloN-30)
  print*, "  > load .glo number ", XGlo
  
  ThisGlo   = MissVal
  call LoadGloVec (Blank, WorkMapIDLReg, ThisGlo)
  
  do XReg = 1, WorkRegN
    GloLoaded (XGlo,XReg) = ThisGlo (XReg)
  end do
end do

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

end subroutine GrabGlo

!*******************************************************************************
! select number of .glo

subroutine GrabGloN

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

end subroutine GrabGloN

!*******************************************************************************
! select particular .glo to use

subroutine SelectGloOp
  
allocate (UseGlo(SelectGloN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SelectGloOp: Allocation failure #####"
UseGlo = 0

print*, "  > Select particular .glo to use in op, totalling: ", SelectGloN
do XSelect = 1, SelectGloN
    print*, "  > Select .glo number: ", XSelect
    do
	read (*,*,iostat=ReadStatus), ChosenGlo
	if (ReadStatus.LE.0.AND.ChosenGlo.GE.1.AND.ChosenGlo.LE.(XGloFree-1)) exit
    end do
    UseGlo (XSelect) = ChosenGlo
end do  

end subroutine SelectGloOp

!*******************************************************************************
! subtract Two from One

subroutine TwoFromOne

  print*, "  > Calc abs diff (=1), % diff (=2), fraction (=3), reduced mag (=4) ?"
  do
	read (*,*,iostat=ReadStatus), QAnom
	if (ReadStatus.LE.0.AND.QAnom.GE.1.AND.QAnom.LE.4) exit
  end do

  if (QAnom.EQ.1) then
   do XReg = 1, WorkRegN
    if (GloLoaded(UseGlo(1),XReg).NE.MissVal.AND.GloLoaded(UseGlo(2),XReg).NE.MissVal) then
      GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg) - GloLoaded(UseGlo(2),XReg)
    end if
   end do
  else if (QAnom.EQ.2) then
   do XReg = 1, WorkRegN
    if (GloLoaded(UseGlo(1),XReg).NE.MissVal.AND.GloLoaded(UseGlo(2),XReg).NE.MissVal) then
     if (GloLoaded(UseGlo(2),XReg).NE.0) then
      GloLoaded (XGloFree,XReg) = 100.0 * (GloLoaded(UseGlo(1),XReg) - &
      				GloLoaded(UseGlo(2),XReg)) / abs (GloLoaded(UseGlo(2),XReg))
     end if
    end if
   end do
  else if (QAnom.EQ.3) then
   do XReg = 1, WorkRegN
    if (GloLoaded(UseGlo(1),XReg).NE.MissVal.AND.GloLoaded(UseGlo(2),XReg).NE.MissVal) then
     if (GloLoaded(UseGlo(2),XReg).NE.0) then
      GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg) / GloLoaded(UseGlo(2),XReg)
     end if
    end if
   end do
  else if (QAnom.EQ.4) then
   print*, "  >   abs(A)<B-->0 else A>0-->A-B else A<0-->A+B"
   do XReg = 1, WorkRegN
    if (GloLoaded(UseGlo(1),XReg).NE.MissVal.AND.GloLoaded(UseGlo(2),XReg).NE.MissVal) then
     if (abs(GloLoaded(UseGlo(1),XReg)).LE.GloLoaded(UseGlo(2),XReg)) then
       GloLoaded (XGloFree,XReg) = 0.0
     else
      if (GloLoaded(1,XReg).GT.0) then
       GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg) - GloLoaded(UseGlo(2),XReg)
      else
       GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg) + GloLoaded(UseGlo(2),XReg)
      end if
     end if
    end if
   end do
  end if

print*, "  > We have stored the result in column: ", XGloFree
if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1

end subroutine TwoFromOne

!*******************************************************************************
! calc Mean

subroutine CalcMean
   
  do XReg = 1, WorkRegN    
   RunTot = 0.0
   RunEn  = 0.0
   
   do XGlo = 1, SelectGloN
    if (GloLoaded(UseGlo(XGlo),XReg).NE.MissVal) then
      RunTot = RunTot + GloLoaded(UseGlo(XGlo),XReg)
      RunEn  = RunEn  + 1.0
    end if
   end do
   
   if (RunEn.GE.1) then
     GloLoaded (XGloFree,XReg) = RunTot / RunEn
   end if
  end do

print*, "  > We have stored the result in column: ", XGloFree
if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1

end subroutine CalcMean

!*******************************************************************************
! calc Median

subroutine CalcMedian
   
allocate (OpVec (SelectGloN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcMedian: Allocation failure #####"
OpVec = MissVal

do XReg = 1, WorkRegN    
   do XGlo = 1, SelectGloN
     OpVec(XGlo) = GloLoaded(UseGlo(XGlo),XReg)
   end do
   
   GloLoaded (XGloFree,XReg) = MedianValue (SelectGloN,OpVec)
end do

print*, "  > We have stored the result in column: ", XGloFree
if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1

end subroutine CalcMedian

!*******************************************************************************
! calc StDev

subroutine CalcStDev
   
  if (WorkGloN.GE.2) then
    print*, "  > Calc sample (=1) or population (=2) stdev ?"
    do
	read (*,*,iostat=ReadStatus), QSampPop
	if (ReadStatus.LE.0.AND.QSampPop.GE.1.AND.QSampPop.LE.2) exit
    end do
  else
    print*, "  > n=1 therefore calc sample stdev"
    QSampPop = 1
  end if

  do XReg = 1, WorkRegN    
    RunSqTot = 0.0
    RunTot   = 0.0
    RunEn    = 0.0
   
    do XGlo = 1, SelectGloN
     if (GloLoaded(UseGlo(XGlo),XReg).NE.MissVal) then
      RunSqTot = RunSqTot + (GloLoaded(UseGlo(XGlo),XReg) * GloLoaded(UseGlo(XGlo),XReg))
      RunTot   = RunTot   +  GloLoaded(UseGlo(XGlo),XReg)
      RunEn    = RunEn    +  1.0
     end if
    end do
   
    if (RunEn.GE.2) then
     if (QSampPop.EQ.1) then
       GloLoaded (XGloFree,XReg) = (RunSqTot/RunEn) - ((RunTot/RunEn)*(RunTot/RunEn))
       GloLoaded (XGloFree,XReg) = sqrt (GloLoaded (XGloFree,XReg))
     else
       GloLoaded (XGloFree,XReg) = (RunSqTot/RunEn) - ((RunTot/RunEn)*(RunTot/RunEn))
       GloLoaded (XGloFree,XReg) = GloLoaded (XGloFree,XReg) * RunEn / (RunEn - 1)
       GloLoaded (XGloFree,XReg) = sqrt (GloLoaded (XGloFree,XReg))
     end if
    end if
  end do

print*, "  > We have stored the result in column: ", XGloFree
if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1

end subroutine CalcStDev

!*******************************************************************************
! calc range

subroutine CalcRange
   
allocate (GloMin (WorkRegN), &
	  GloMax (WorkRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcRange: Allocation failure #####"
GloMin = MissVal
GloMax = MissVal

do XReg = 1, WorkRegN    
   RunMin   =  1000000.0
   RunMax   = -1000000.0

   do XGlo = 1, SelectGloN
    if (GloLoaded(UseGlo(XGlo),XReg).NE.MissVal) then
      if (GloLoaded(UseGlo(XGlo),XReg).LT.RunMin) then
        GloMin (XReg) = GloLoaded(UseGlo(XGlo),XReg)
        RunMin        = GloLoaded(UseGlo(XGlo),XReg)
      end if
      
      if (GloLoaded(UseGlo(XGlo),XReg).GT.RunMax) then
        GloMax (XReg) = GloLoaded(UseGlo(XGlo),XReg)
        RunMax        = GloLoaded(UseGlo(XGlo),XReg)
      end if
    end if
   end do
end do

do XReg = 1, WorkRegN			! store min values
  GloLoaded (XGloFree,XReg) = GloMin (XReg)
end do
if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1

do XReg = 1, WorkRegN			! store max values
  GloLoaded (XGloFree,XReg) = GloMax (XReg)
end do
if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1

print*, "  > We have stored the min and max in columns: ", (XGloFree-2), (XGloFree-1)

deallocate (GloMin,GloMax, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CalcRange: Deallocation failure #####"

end subroutine CalcRange

!*******************************************************************************
! calc LSR

subroutine CalcLSR
   
SumX  = 0.0
SumY  = 0.0
SumXX = 0.0
SumYY = 0.0
SumXY = 0.0
En    = 0.0

Aye = MissVal
Bee = MissVal
Are = MissVal

do XReg = 1, WorkRegN    
 if (GloLoaded(UseGlo(1),XReg).NE.MissVal.AND.GloLoaded(UseGlo(2),XReg).NE.MissVal) then
   SumX  = SumX  + GloLoaded(UseGlo(1),XReg)
   SumY  = SumY  + GloLoaded(UseGlo(2),XReg)
   SumXX = SumXX + GloLoaded(UseGlo(1),XReg) ** 2
   SumYY = SumYY + GloLoaded(UseGlo(2),XReg) ** 2
   SumXY = SumXY + GloLoaded(UseGlo(1),XReg) *  GloLoaded(UseGlo(2),XReg)
   En    = En    + 1.0
 end if
end do
   
if (En.GE.2) then
    Num = (En*SumXY)-(SumX*SumY)
    Den = (En*SumXX)-(SumX*SumX)
    if (Den.NE.0) Aye = Num / Den

    if (Aye.NE.MissVal) Bee = (SumY-(Aye*SumX))/En
    
    Num = (SumXY/En)-((SumX/En)*(SumY/En)) 
    Den = sqrt((SumXX/En)-((SumX/En)*(SumX/En)))*sqrt((SumYY/En)-((SumY/En)*(SumY/En)))  
    if (Den.NE.0) Are = Num / Den
    
    print*, "  > Values of a,b,r for y=ax+b, where x= value in first .glo"
    print "(3F12.4)", Aye, Bee, Are
else
    print*, "  > Too few valid regions."
end if

end subroutine CalcLSR

!*******************************************************************************
! pos of A relative to B and C

subroutine PosABC
   
 do XReg = 1, WorkRegN    
   if (GloLoaded(UseGlo(1),XReg).NE.MissVal.AND.GloLoaded(UseGlo(2),XReg).NE.MissVal &
   					.AND.GloLoaded(UseGlo(3),XReg).NE.MissVal) then
     if (GloLoaded(UseGlo(1),XReg).LT.GloLoaded(UseGlo(2),XReg) &
     					.AND.GloLoaded(UseGlo(1),XReg).LT.GloLoaded(UseGlo(3),XReg)) then
       GloLoaded (XGloFree,XReg) = -1
     else if (GloLoaded(1,XReg).GT.GloLoaded(UseGlo(2),XReg) &
     					.AND.GloLoaded(UseGlo(1),XReg).GT.GloLoaded(UseGlo(3),XReg)) then
       GloLoaded (XGloFree,XReg) = 1
     else
       GloLoaded (XGloFree,XReg) = 0
     end if
   else
       GloLoaded (XGloFree,XReg) = MissVal
   end if
 end do

print*, "  > We have stored the result in column: ", XGloFree
if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1

end subroutine PosABC

!*******************************************************************************
! mask A using B

subroutine MaskAB

if ((XGloFree+2).LE.WorkGloN) then   
 print*, "  > Is the threshold a real (=1) or a magnitude (=2): "
 do
	read (*,*,iostat=ReadStatus), CutOffType
	if (ReadStatus.LE.0.AND.CutOffType.GE.1.AND.CutOffType.LE.2) exit
 end do

 print*, "  > Enter the threshold in B by which to separate A: "
 do
	read (*,*,iostat=ReadStatus), CutOffVal
	if (ReadStatus.LE.0) exit
 end do

 do XReg = 1, WorkRegN    
   if (GloLoaded(UseGlo(1),XReg).NE.MissVal.AND.GloLoaded(UseGlo(2),XReg).NE.MissVal) then
       if (CutOffType.EQ.2) then
       		TestVal = abs (GloLoaded(UseGlo(2),XReg))
       else
       		TestVal = GloLoaded(UseGlo(2),XReg)
       end if
       
       if      (TestVal.EQ.CutOffVal) then
       	  GloLoaded ((XGloFree+1),XReg) = GloLoaded(UseGlo(1),XReg)
       else if (TestVal.LT.CutOffVal) then
       	  GloLoaded ((XGloFree+0),XReg) = GloLoaded(UseGlo(1),XReg)
       else if (TestVal.GT.CutOffVal) then
       	  GloLoaded ((XGloFree+2),XReg) = GloLoaded(UseGlo(1),XReg)
       end if
   end if
 end do
 
 print "(a56,3i4)", "  > We have stored values <, =, > threshold in columns: ", &
 				XGloFree, (XGloFree+1), (XGloFree+2)
 if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 3
 if (XGloFree.GT.WorkGloN) XGloFree = WorkGloN
else
 print*, "  > Insufficient columns remaining to perform operation."
end if

end subroutine MaskAB

!*******************************************************************************
! compare two masked columns

subroutine CompMasked

if ((XGloFree+1).LE.WorkGloN) then   
 print*, "  > Select the min (=1) or max (=2) of A and B:"
 do
	read (*,*,iostat=ReadStatus), QMinMax
	if (ReadStatus.LE.0.AND.QMinMax.GE.1.AND.QMinMax.LE.2) exit
 end do
  
 do XReg = 1, WorkRegN    
   if (GloLoaded(UseGlo(1),XReg).EQ.MissVal) then
     if (GloLoaded(UseGlo(2),XReg).EQ.MissVal) then
       GloLoaded ((XGloFree+0),XReg) = MissVal     
       GloLoaded ((XGloFree+1),XReg) = 0     
     else
       GloLoaded ((XGloFree+0),XReg) = GloLoaded(UseGlo(2),XReg)     
       GloLoaded ((XGloFree+1),XReg) = 1     
     end if
   else
     if (GloLoaded(UseGlo(2),XReg).EQ.MissVal) then
       GloLoaded ((XGloFree+0),XReg) = GloLoaded(UseGlo(1),XReg)     
       GloLoaded ((XGloFree+1),XReg) = 2     
     else
       OpDiff = GloLoaded(UseGlo(2),XReg) - GloLoaded(UseGlo(1),XReg)		! B minus A
       
       if (OpDiff.GE.0) then
         if (QMinMax.EQ.1) then
           GloLoaded ((XGloFree+0),XReg) = GloLoaded(UseGlo(1),XReg)     
           GloLoaded ((XGloFree+1),XReg) = 3     
         else
           GloLoaded ((XGloFree+0),XReg) = GloLoaded(UseGlo(2),XReg)     
           GloLoaded ((XGloFree+1),XReg) = 4     
         end if
       else
         if (QMinMax.EQ.1) then
           GloLoaded ((XGloFree+0),XReg) = GloLoaded(UseGlo(2),XReg)     
           GloLoaded ((XGloFree+1),XReg) = 4     
         else
           GloLoaded ((XGloFree+0),XReg) = GloLoaded(UseGlo(1),XReg)     
           GloLoaded ((XGloFree+1),XReg) = 3     
         end if
       end if
     end if
   end if  
 end do
 
 print*, "  > Selected value and code go in columns:", XGloFree, (XGloFree+1)
 print*, "  > A,B=MissVal value=MissVal code=0"
 print*, "  >   A=MissVal value=B       code=1"
 print*, "  >   B=MissVal value=A       code=2"
 print*, "  >   A=chosen  value=A       code=3"
 print*, "  >   B=chosen  value=B       code=4"
  
 XGloFree = XGloFree + 2
 if (XGloFree.GT.WorkGloN) XGloFree = WorkGloN
else
 print*, "  > Insufficient columns remaining to perform operation."
end if

end subroutine CompMasked

!*******************************************************************************
! divide by constant

subroutine OperateAConstant

print*, "  > Divide=1, multiply=2, minus=3, add=4, sqroot=5, exp=6, abs=7"
do
	read (*,*,iostat=ReadStatus), Operator
	if (ReadStatus.LE.0 .AND. Operator.GE.1 .AND. Operator.LE.7) exit
end do

if (Operator.LT.5) then
  print*, "  > Enter constant (real):"
  do
	read (*,*,iostat=ReadStatus), RealConstant
	if (ReadStatus.LE.0) exit
  end do
else if (Operator.EQ.6) then
  print*, "  > Enter constant (positive integer):"
  do
	read (*,*,iostat=ReadStatus), IntConstant
	if (ReadStatus.LE.0.AND.IntConstant.GE.0) exit
  end do
end if

do XReg = 1, WorkRegN
   if (GloLoaded(UseGlo(1),XReg).NE.MissVal) then        
        if (Operator.EQ.1.AND.RealConstant.NE.0) &
        		   GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg)  /  RealConstant
        if (Operator.EQ.2) GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg)  *  RealConstant
        if (Operator.EQ.3) GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg)  -  RealConstant
        if (Operator.EQ.4) GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg)  +  RealConstant
        if (Operator.EQ.5) GloLoaded (XGloFree,XReg) = sqrt (GloLoaded(UseGlo(1),XReg))
        if (Operator.EQ.6) GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg)  ** IntConstant
        if (Operator.EQ.7) GloLoaded (XGloFree,XReg) = abs  (GloLoaded(UseGlo(1),XReg))
   end if
end do

print*, "  > We have stored the result in column: ", XGloFree
if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1

end subroutine OperateAConstant

!*******************************************************************************
! operate on A with B

subroutine OperateAB

print*, "  > A operator B"
print*, "  > Operator is divide =1, multi =2, minus =3, add =4"
do
	read (*,*,iostat=ReadStatus), Operator
	if (ReadStatus.LE.0 .AND. Operator.GE.1 .AND. Operator.LE.4) exit
end do

do XReg = 1, WorkRegN
   if (GloLoaded(UseGlo(1),XReg).NE.MissVal.AND.GloLoaded(UseGlo(2),XReg).NE.MissVal) then        
        if (Operator.EQ.1.AND.GloLoaded(UseGlo(2),XReg).NE.0) &
        		   GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg) / GloLoaded(UseGlo(2),XReg)
        if (Operator.EQ.2) GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg) * GloLoaded(UseGlo(2),XReg)
        if (Operator.EQ.3) GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg) - GloLoaded(UseGlo(2),XReg)
        if (Operator.EQ.4) GloLoaded (XGloFree,XReg) = GloLoaded(UseGlo(1),XReg) + GloLoaded(UseGlo(2),XReg)
   end if
end do

print*, "  > We have stored the result in column: ", XGloFree
if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1

end subroutine OperateAB

!*******************************************************************************
! get class specifications

subroutine GetClassSpec

print*, "  > Define equal-interval classes (=0) or manually (=1) ?"
do
	read (*,*,iostat=ReadStatus), QClassMethod
	if (ReadStatus.LE.0.AND.QClassMethod.GE.0.AND.QClassMethod.LE.1) exit
end do

if (QClassMethod.EQ.0) then
  call EqualClass
else
  call ManualClass
end if

end subroutine GetClassSpec

!*******************************************************************************
! get class specifications

subroutine EqualClass

print*, "  > Enter the min,max,interval: "
do
	read (*,*,iostat=ReadStatus), RangeMin,RangeMax,RangeGap
	if (ReadStatus.LE.0 .AND. RangeMax.GT.RangeMin .AND. &
		mod((RangeMax-RangeMin),RangeGap).EQ.0.0) exit
end do

ClassN = nint((RangeMax-RangeMin) / RangeGap)
allocate (UseClass(ClassN,3), stat=AllocStat)		! lower & upper bounds, class value
if (AllocStat.NE.0) print*, "  > ##### ERROR: EqualClass: Allocation failure #####"

do XClass=1,ClassN
  UseClass (XClass,1) = RangeMin + ((XClass-1)*RangeGap)
  UseClass (XClass,2) = RangeMin + ((XClass+0)*RangeGap)
  UseClass (XClass,3) = (UseClass (XClass,1) + UseClass (XClass,2)) / 2.0
end do

end subroutine EqualClass

!*******************************************************************************
! get class specifications

subroutine ManualClass

print*, "  > Enter the number of classes:"
do
	read (*,*,iostat=ReadStatus), ClassN
	if (ReadStatus.LE.0 .AND. ClassN.GE.1) exit
end do

allocate (UseClass(ClassN,3), stat=AllocStat)		! lower & upper bounds, class value
if (AllocStat.NE.0) print*, "  > ##### ERROR: ManualClass: Allocation failure #####"
UseClass = MissVal

print*, "  > Enter each class boundary, in ascending order, totalling:", (ClassN-1)
LastBound = 10 - huge(ClassBound)
do XBound = 1, (ClassN-1)
  do
  	read (*,*,iostat=ReadStatus), ClassBound
	if (ReadStatus.NE.0 .OR.  ClassBound.LE.LastBound) print*, "  > Boundary not accepted. Re-enter it."
	if (ReadStatus.LE.0 .AND. ClassBound.GT.LastBound) exit
  end do
  
  UseClass (XBound,2)     = ClassBound
  UseClass ((XBound+1),1) = ClassBound
end do

print*, "  > Enter each class value, in class order, totalling:", ClassN
do XClass = 1, ClassN
  do
  	read (*,*,iostat=ReadStatus), ClassValue
	if (ReadStatus.NE.0) print*, "  > Value not accepted. Re-enter it."
	if (ReadStatus.LE.0) exit
  end do
  
  UseClass (XClass,3) = ClassValue
end do

end subroutine ManualClass

!*******************************************************************************
! classify A and assign a class value to each box

subroutine ClassifyABox

call GetClassSpec

do XGlo = 1, SelectGloN
 do XReg = 1, WorkRegN
  if (GloLoaded(UseGlo(XGlo),XReg).NE.MissVal) then        
    XClass  = 0    
    ExitNow = 0
    
    do
      XClass = XClass + 1

      if (GloLoaded(UseGlo(XGlo),XReg).LE.UseClass(XClass,2).OR.XClass.EQ.ClassN) then
        	GloLoaded(XGloFree,XReg) = UseClass(XClass,3)
        	ExitNow = 1
      end if

      if (ExitNow.EQ.1) exit
    end do
  end if
 end do

 print*, "  > Column classified to column: ", UseGlo(XGlo), XGloFree
 if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1
end do

deallocate (UseClass,UseGlo, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ClassifyABox: Deallocation failure #####"

end subroutine ClassifyABox

!*******************************************************************************
! classify A and count frequency in each group

subroutine ClassifyAGroup

call GetClassSpec

allocate (FreqClass(ClassN), stat=AllocStat)		! no. of values in class
if (AllocStat.NE.0) print*, "  > ##### ERROR: ClassifyAGroup: Allocation failure #####"

print "(a22)", " col class value  freq"

do XGlo = 1, SelectGloN
 FreqClass = 0
 
 do XReg = 1, WorkRegN
  if (GloLoaded(UseGlo(XGlo),XReg).NE.MissVal) then        
    XClass  = 0    
    ExitNow = 0
    
    do
      XClass = XClass + 1

      if (GloLoaded(UseGlo(XGlo),XReg).LE.UseClass(XClass,2).OR.XClass.EQ.ClassN) then
        	FreqClass(XClass) = FreqClass(XClass) + 1
        	ExitNow = 1
      end if

      if (ExitNow.EQ.1) exit
    end do
  end if
 end do

 print*
 do XClass = 1, ClassN
  print "(i4,f12.4,i6)", UseGlo(XGlo), UseClass(XClass,3), FreqClass(XClass)
 end do
end do

deallocate (UseClass,FreqClass, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: ClassifyABox: Deallocation failure #####"

end subroutine ClassifyAGroup

!*******************************************************************************
! compare A and B

subroutine CompareAB

ClassN = 3
allocate (UseClass(ClassN,1), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CompareAB: Allocation failure #####"
UseClass = MissVal

print*, "  > We classify by: A<B, A=B, A>B. Enter a value to assign to each."
do XClass = 1, ClassN
  do
  	read (*,*,iostat=ReadStatus), ClassValue
	if (ReadStatus.NE.0) print*, "  > Value not accepted. Re-enter it."
	if (ReadStatus.LE.0) exit
  end do
  
  UseClass (XClass,1) = ClassValue
end do

do XReg = 1, WorkRegN
  if (GloLoaded(UseGlo(1),XReg).NE.MissVal.AND.GloLoaded(UseGlo(2),XReg).NE.MissVal) then        
      if      (GloLoaded(UseGlo(1),XReg).LT.GloLoaded(UseGlo(2),XReg)) then
        	GloLoaded(XGloFree,XReg) = UseClass(1,1)
      else if (GloLoaded(UseGlo(1),XReg).EQ.GloLoaded(UseGlo(2),XReg)) then
        	GloLoaded(XGloFree,XReg) = UseClass(2,1)
      else
        	GloLoaded(XGloFree,XReg) = UseClass(3,1)
      end if
  else
                GloLoaded(XGloFree,XReg) = MissVal
  end if
end do

print*, "  > We have stored the result in column: ", XGloFree
if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1

deallocate (UseClass, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: CompareAB: Deallocation failure #####"

end subroutine CompareAB

!*******************************************************************************
! regionalise A using new .ref
! no missing values are permitted, but this could be altered fairly simply

subroutine RegionaliseA

call RegSelect (WorkGrid, WorkLongN, WorkLatN, WorkDataN, NewMapIDLReg, &
			NewRegSizes, NewRegNames, NewRegTitle, NewRegN)

print*, "  > Convert as totals (=1) or means (=2) ? "
do
	read (*,*,iostat=ReadStatus), QConvert
	if (ReadStatus.LE.0.AND.QConvert.GE.1.AND.QConvert.LE.2) exit
end do

print*, "  > Save to .glo (=1), print to screen (=2), or both (=3) ? "
do
	read (*,*,iostat=ReadStatus), QOutput
	if (ReadStatus.LE.0.AND.QOutput.GE.1.AND.QOutput.LE.3) exit
end do

allocate (NewTot (NewRegN), &
	  NewEn  (NewRegN), &
	  ThisGlo(NewRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: RegionaliseA: Allocation failure #####"
NewTot = 0.0
NewEn  = 0.0
ThisGlo= MissVal

do XLong = 1, WorkLongN
  do XLat = 1, WorkLatN
    if (NewMapIDLReg (XLong,XLat).NE.MissVal) then		! grid box is wanted for new region set    
      if (WorkMapIDLReg (XLong,XLat).NE.MissVal) then		! grid box is specified in old region set    	
    	if (GloLoaded ( UseGlo(1), WorkMapIDLReg (XLong,XLat) ).NE.MissVal) then	! value is valid    	
    	  NewTot ( NewMapIDLReg (XLong,XLat) ) = NewTot ( NewMapIDLReg (XLong,XLat) ) + &
    	  						GloLoaded ( UseGlo(1), WorkMapIDLReg (XLong,XLat) )
          NewEn  ( NewMapIDLReg (XLong,XLat) ) = NewEn  ( NewMapIDLReg (XLong,XLat) ) + 1.0
        end if
      end if    
    end if    
  end do
end do

do XReg = 1, NewRegN
  if (NewEn(XReg).EQ.NewRegSizes(XReg)) then
    if (QConvert.EQ.1) then
      ThisGlo (XReg) = NewTot (XReg)
    else
      ThisGlo (XReg) = NewTot (XReg) / NewEn (XReg)
    end if
  end if
end do

if (QOutput.NE.2) call SaveGlo (WorkLongN,WorkLatN,NewRegN,WorkGridFilePath,Blank,Blank,ThisGlo,NewMapIDLReg)
if (QOutput.NE.1) then
  print '(a6,a1,a20,a1,a6,a12)', ' Index', ' ', '              Region', ' ', &
    	' Boxes', '       Value'
  do XReg = 1, NewRegN
    print '(i6,a1,a20,a1,i6,f20.4)', XReg, ' ', NewRegNames(XReg), ' ', &
    	NewRegSizes(XReg), ThisGlo (XReg)
  end do
end if

write (99,'(a6,a1,a20,a1,a6,a12)'), ' Index', ' ', '              Region', ' ', &
    	' Boxes', '       Value'
do XReg = 1, NewRegN
    write (99,'(i6,a1,a20,a1,i6,f20.4)'), XReg, ' ', NewRegNames(XReg), ' ', &
    	NewRegSizes(XReg), ThisGlo (XReg)
end do

deallocate (NewMapIDLReg, NewRegSizes, NewRegNames, NewTot, NewEn, ThisGlo, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: RegionaliseA: Deallocation failure #####"

end subroutine RegionaliseA

!*******************************************************************************
! regionalise A using new .ref for percentage column
! no missing values are permitted, but this could be altered fairly simply

subroutine RegPerCentA

print*, "  > Load the base values against which percentages were calculated."
call LoadGloVec (Blank, WorkMapIDLReg, BaseGlo, Silent=1)
  
call RegSelect (WorkGrid, WorkLongN, WorkLatN, WorkDataN, NewMapIDLReg, &
			NewRegSizes, NewRegNames, NewRegTitle, NewRegN)

print*, "  > Save to .glo (=1), print to screen (=2), or both (=3) ? "
do
	read (*,*,iostat=ReadStatus), QOutput
	if (ReadStatus.LE.0.AND.QOutput.GE.1.AND.QOutput.LE.3) exit
end do

allocate (PertRaw(WorkRegN),&
	  PertTot(NewRegN), &
	  BaseTot(NewRegN), &
	  NewEn  (NewRegN), &
	  ThisGlo(NewRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: RegPerCentA: Allocation failure #####"
PertTot = 0.0
BaseTot = 0.0
NewEn   = 0.0
ThisGlo = MissVal
PertRaw = MissVal

do XReg = 1, WorkRegN
  if (GloLoaded (UseGlo(1),XReg).NE.MissVal.AND.BaseGlo(XReg).NE.MissVal) &
  	PertRaw(XReg) = ((GloLoaded(UseGlo(1),XReg) * BaseGlo(XReg)) / 100.0) + BaseGlo(XReg)
end do

do XLong = 1, WorkLongN
  do XLat = 1, WorkLatN
    if (NewMapIDLReg (XLong,XLat).NE.MissVal) then		! grid box is wanted for new region set    
      if (WorkMapIDLReg (XLong,XLat).NE.MissVal) then		! grid box is specified in old region set    	
       if (PertRaw (WorkMapIDLReg(XLong,XLat)).NE.MissVal) then		! raw  is valid    	
        if (BaseGlo (WorkMapIDLReg(XLong,XLat)).NE.MissVal) then	! base is valid
    	  PertTot (NewMapIDLReg(XLong,XLat)) = PertTot (NewMapIDLReg(XLong,XLat)) + &
    	  						PertRaw (WorkMapIDLReg(XLong,XLat))
    	  BaseTot (NewMapIDLReg(XLong,XLat)) = BaseTot (NewMapIDLReg(XLong,XLat)) + &
    	  						BaseGlo (WorkMapIDLReg(XLong,XLat))
          NewEn   (NewMapIDLReg(XLong,XLat)) = NewEn   (NewMapIDLReg(XLong,XLat)) + 1.0
        end if
       end if
      end if    
    end if    
  end do
end do

do XReg = 1, NewRegN
  if (NewEn(XReg).EQ.NewRegSizes(XReg)) then
      ThisGlo (XReg) = 100.0 * (PertTot(XReg)-BaseTot(XReg)) / BaseTot(XReg)
  end if
end do

if (QOutput.NE.2) call SaveGlo (WorkLongN,WorkLatN,NewRegN,WorkGridFilePath,Blank,Blank,ThisGlo,NewMapIDLReg)
if (QOutput.NE.1) then
  print '(a6,a1,a20,a1,a6,a12)', ' Index', ' ', '              Region', ' ', &
    	' Boxes', '       Value'

  do XReg = 1, NewRegN
    print '(i6,a1,a20,a1,i6,f12.4)', XReg, ' ', NewRegNames(XReg), ' ', &
    	NewRegSizes(XReg), ThisGlo (XReg)
  end do
end if

deallocate (NewMapIDLReg, NewRegSizes, NewRegNames, NewEn, ThisGlo, BaseGlo, &
			PertRaw, PertTot, BaseTot, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: RegPerCentA: Deallocation failure #####"

end subroutine RegPerCentA

!*******************************************************************************
! transform one value into another

subroutine TransformVal
 
print*, "  > Enter the value to transform wherever it occurs: "
do
	read (*,*,iostat=ReadStatus), ValueToChange
	if (ReadStatus.LE.0) exit
end do

print*, "  > Enter the new value: "
do
	read (*,*,iostat=ReadStatus), ValueToLeave
	if (ReadStatus.LE.0) exit
end do

ValidN = 0    
do XReg = 1, WorkRegN
  if (GloLoaded(UseGlo(1),XReg).EQ.ValueToChange) then
  	GloLoaded(XGloFree,XReg) = ValueToLeave
  	ValidN = ValidN + 1
  else
  	GloLoaded(XGloFree,XReg) = GloLoaded(UseGlo(1),XReg)
  end if
end do    
print*, "  > Values changed total: ", ValidN

print*, "  > We have stored the result in column: ", XGloFree
if (XGloFree.LT.WorkGloN) XGloFree = XGloFree + 1

end subroutine TransformVal

!*******************************************************************************
! save .glo file

subroutine SaveToGlo

print*, "  > Enter the column to save as a .glo: "
do
	read (*,*,iostat=ReadStatus), XGloSave
	if (ReadStatus.LE.0.AND.XGloSave.GE.1.AND.XGloSave.LE.WorkGloN) exit
end do

allocate (GloSaved(WorkRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveToGlo: Allocation failure #####"
GloSaved = MissVal

do XReg = 1, WorkRegN
  GloSaved(XReg) = GloLoaded(XGloSave,XReg)
end do

call SaveGlo (WorkLongN,WorkLatN,WorkRegN,WorkGridFilePath,Blank,Blank,GloSaved,WorkMapIDLReg)

deallocate (GloSaved, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: SaveToGlo: Deallocation failure #####"

end subroutine SaveToGlo

!*******************************************************************************
! print details to screen

subroutine PrintToScreen

print "(a4,a6,5a12)", "File", "   Reg", "         Sum", "  Sum Squares", "        Mean", "       StDev", &
					"       AvMag" 

do XGlo = 1, (XGloFree-1)
  OpTot   = 0.0
  OpSqTot = 0.0
  OpAbsTot= 0.0
  OpEn    = 0.0
  OpMean  = MissVal
  OpStDev = MissVal
  AvMag   = MissVal
  
  do XReg = 1, WorkRegN
    if (GloLoaded(XGlo,XReg).NE.MissVal) then
      OpTot    = OpTot    +  GloLoaded(XGlo,XReg)
      OpSqTot  = OpSqTot  + (GloLoaded(XGlo,XReg) ** 2)
      OpAbsTot = OpAbsTot + abs(GloLoaded(XGlo,XReg))
      OpEn     = OpEn     +  1.0
    end if
  end do
  
  RegKept = nint (OpEn)
  
  if (OpEn.GE.1) then
  	OpMean  =  OpTot / OpEn
  	OpStDev = (OpSqTot / OpEn) - (OpMean ** 2)
  	OpStDev = sqrt (OpStDev)
        AvMag   = OpAbsTot / OpEn
  end if
  
  print "(i4,i6,5e12.4)", XGlo, RegKept, OpTot, OpSqTot, OpMean, OpStDev, AvMag
end do

end subroutine PrintToScreen

!*******************************************************************************
! print extra details to screen

subroutine PrintToScreenExtra

allocate (ThisGlo(WorkRegN), &
	  Indices(WorkRegN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PrintToScreenExtra: Allocation failure #####"

print "(a4,a6,4a12)", "File", "   Reg", "      MinVal", "      MaxVal", "      Median", "        Mode"

do XGlo = 1, (XGloFree-1)
  MinVal  =  100000.0
  MaxVal  = -100000.0
  Median  =  MissVal
  Mode    =  MissVal
  OpEn    =  0.0
  OpTot   =  0.0
  
  do XReg = 1, WorkRegN
    ThisGlo(XReg) = GloLoaded(XGlo,XReg)
    
    if (GloLoaded(XGlo,XReg).NE.MissVal) then
      if (GloLoaded(XGlo,XReg).GT.MaxVal) MaxVal = GloLoaded(XGlo,XReg)
      if (GloLoaded(XGlo,XReg).LT.MinVal) MinVal = GloLoaded(XGlo,XReg)
      
      OpEn     = OpEn     +  1.0 
      OpTot    = OpTot    +  GloLoaded(XGlo,XReg)
    end if
  end do
  
  Median = MedianValue (WorkRegN,ThisGlo)
  Mode   = (OpTot / OpEn) - (3 * ( (OpTot / OpEn) - Median ) )
  
  RegKept = nint (OpEn)
  
  print "(i4,i6,4e12.4)", XGlo, RegKept, MinVal, MaxVal, Median, Mode
end do

deallocate (ThisGlo,Indices, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: PrintToScreenExtra: Deallocation failure #####"

end subroutine PrintToScreenExtra

!*******************************************************************************
! write region details to log file

subroutine WriteRegLog

do XReg = 1, WorkRegN
    write (99,'(i6,a1,a20,a1,i6,e12.5)'), XReg, ' ', WorkRegNames(XReg), ' ', &
    	WorkRegSizes(XReg), GloLoaded (UseGlo(1),XReg)
end do

end subroutine WriteRegLog

!*******************************************************************************
! write region details to the screen

subroutine WriteRegScreen

print '(a6,a1,a20,a1,a6,a12)', ' Index', ' ', '              Region', ' ', &
    	' Boxes', '       Value'

do XReg = 1, WorkRegN
    print '(i6,a1,a20,a1,i6,f12.4)', XReg, ' ', WorkRegNames(XReg), ' ', &
    	WorkRegSizes(XReg), GloLoaded (UseGlo(1),XReg)
end do

end subroutine WriteRegScreen

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

end program CompareGlo2