; MakeRegions.pro
; written by Tim Mitchell on 3.8.00
; last modified on 16.8.00

; ******************************************************************************
; initialise

initx

MissVal  = -999.0
ImageLatMin  = 0.0 & ImageLatMax  = 0.0 & ImageLongMin = 0.0 & ImageLongMax = 0.0
Black = 0 & White = 255 & Grey = 254
DatumCha = "" & DatumInt = 0L & DatumFlt = 0.0D
View0 = 0 & View1 = 0 & View2 = 0 & View3 = 0
Outline0 = 0 & Outline1 = 0 & Outline2 = 0
ALat = 0 & BLat = 0 & CLat = 0 & ALong = 0 & BLong = 0 & CLong = 0
MinLong = 0L & MinLat = 0L & MaxLong = 0L & MaxLat = 0L
TotPlus = 0L & TotMinus = 0L
SelInProg = 0

ModelChosen = -999
SelectModel, ModelChosen, ModelFilePath, ModelLongN, ModelLatN
ModelBoxN = 0L
ModelBoxN = ModelLongN * ModelLatN
  
GetImageBounds, ModelChosen, ImageLatMin, ImageLatMax, ImageLongMin, ImageLongMax
GetViewBounds, ViewN, ViewName, ViewKeyWords, ViewBounds
GetScales, ScaleN, ScaleName, ScaleColTab, ScaleSeg, ScaleLimits

SelArray        = lonarr (ModelLongN,ModelLatN)		; contains region number
SelArray [*,*]  = long(MissVal)
InfoArray       = dblarr (ModelLongN,ModelLatN)		; contains .glo file loaded
InfoArray [*,*] = MissVal

RegName       = strarr (ModelBoxN)			; name up to 20 characters
RegName [*]   = ""
RegSize       = lonarr (ModelBoxN)			; number of grid boxes
RegSize [*]   = 0
RegColour     = intarr (ModelBoxN)			; colour index (ct=27)
RegColour [*] = White

RegBounds       = dblarr (ModelBoxN,4)			; lonmin, lonmax, latmin, latmax
RegBounds [*,*] = MissVal
Spread          = 5.0

DisplayChosen = 2			; 0=Sel-View, 1=Sel-Reg, 2=Info-View, 3=Info-Reg
ViewChosen    = 0
ScaleChosen   = 0
Outlines      = [1,0,0]			; coasts, countries, rivers; 0=no, 1=yes

ExtraLines      = fltarr (1,4)		; start(long,lat) - end(long,lat)
ExtraLines[*,*] = MissVal
ExtraN = 1

CurrentReg    = 0
FirstFreeReg  = 0

TogglePlots   = 4
P0Long        = 0.0

LastPlotBounds= [-180.0,180.0,-90.0,90.0]
Button        = 0

PlotDisplay   = intarr (4)		; top-left, top-right, bott-left, bott-right
PlotDisplay   = [1,3,0,2]

DisplayBounds = fltarr (4,4)		; second dim: 0=lonmin, 1=lonmax, 2=latmin, 3=latmax
for XDisplay = 0, 3 do begin
	DisplayBounds [XDisplay,*] = ViewBounds [0,*]
endfor

InfoColour       = intarr (ModelLongN,ModelLatN)
InfoColour [*,*] = Grey
SelColour        = intarr (ModelLongN,ModelLatN)
SelColour  [*,*] = Grey

;******************************************************************************
; main loop

Choice = -1

while (Choice LT 99) do begin
	
case 1 of
	
;******************************************************************************
; load existing region specifications

(Choice EQ-1): begin	

; don't print menu

end

;******************************************************************************
; load existing region specifications

(Choice EQ 0): begin	

LoadRef, ModelChosen, ModelLongN, ModelLatN, SelArray, RefSize, RefName, RefColour, RefRegN

RegName [*]     = ""
RegSize [*]     = 0
RegColour [*]   = Grey
RegBounds [*,*] = MissVal

for XReg = 0L, (RefRegN-1) do begin
  RegName [XReg]   = RefName [XReg]
  RegSize [XReg]   = RefSize [XReg]
  RegColour [XReg] = RefColour [XReg]
endfor

CurrentReg   = 0
FirstFreeReg = RefRegN

print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'
end

;******************************************************************************
; load .glo file

(Choice EQ 1): begin	

GloFilePath = ""
LoadGlo, ModelChosen, InfoArray, GloFilePath, GloFileInfo
									; info array coloured
InfoColour = DataToColours (InfoArray, reform(ScaleLimits(ScaleChosen,*)), ModelLongN, ModelLatN, Grey)

end

;******************************************************************************
; change order of plots

(Choice EQ 2): begin	

print, "  > Select display to plot last (-1=list): "
DisplayChosen = -2
while (DisplayChosen LT 0 OR DisplayChosen GE 4) do begin
  if (DisplayChosen EQ -1) then begin
      print, "  > 0=Sel-View, 1=Sel-Reg, 2=Info-View, 3=Info-Reg"
  endif

  read, DisplayChosen
endwhile

if (DisplayChosen EQ 0) then PlotDisplay = [3,1,2,0]
if (DisplayChosen EQ 1) then PlotDisplay = [2,0,3,1]
if (DisplayChosen EQ 2) then PlotDisplay = [1,3,0,2]
if (DisplayChosen EQ 3) then PlotDisplay = [0,2,1,3]

end

;******************************************************************************
; change view plotted

(Choice EQ 3): begin	

print, "  > Select view (-1=list): "
ViewChosen = -2
while (ViewChosen LT 0 OR ViewChosen GE ViewN) do begin
  if (ViewChosen EQ -1) then begin
    for XView = 0, (ViewN-1) do begin
      print, XView, " : ", ViewName[XView], format='(i4,a3,a10)'
    endfor  
  endif

  read, ViewChosen
endwhile

DisplayBounds [0,*] = ViewBounds [ViewChosen,*]
DisplayBounds [2,*] = ViewBounds [ViewChosen,*]

end

;******************************************************************************
; change scale plotted

(Choice EQ 4): begin	

print, "  > Select scale (-1=list): "
ScaleChosen = -2
while (ScaleChosen LT 0 OR ScaleChosen GE ScaleN) do begin		; scale selected
  if (ScaleChosen EQ -1) then begin					; scales listed
    for XScale = 0, (ScaleN-1) do begin
      print, XScale, " : ", ScaleName[XScale], ScaleLimits[XScale,0], ScaleLimits[XScale,1], $
  		 format='(i4,a3,a10,2f10.2)'
    endfor  
  endif

  read, ScaleChosen
endwhile
									; info array coloured
InfoColour = DataToColours (InfoArray, reform(ScaleLimits(ScaleChosen,*)), ModelLongN, ModelLatN, Grey)
   
end

;******************************************************************************
; reverse all scales

(Choice EQ 5): begin	

for XScale = 0, (ScaleN-1) do begin
  DatumFlt = ScaleLimits[XScale,2]
  ScaleLimits[XScale,2] = ScaleLimits[XScale,3]
  ScaleLimits[XScale,3] = DatumFlt
endfor
									; info array coloured
InfoColour = DataToColours (InfoArray, reform(ScaleLimits(ScaleChosen,*)), ModelLongN, ModelLatN, Grey)

end

;******************************************************************************
; select outlines

(Choice EQ 6): begin	

print, "  > Select outlines for coasts (0=N,1=Y), countries (0=N,1=Y), user(0=N): "
read, Outline0, Outline1, Outline2
Outlines = [Outline0, Outline1, Outline2]
if (Outlines[2] NE 0) then SelectPoly, Outlines[2]

end

;******************************************************************************
; load lines to plot

(Choice EQ 7): begin	

DatTitle = '' & DatFilePath = ''
ExtraN = 0
StartX = 0.0 & StartY = 0.0 & EndX = 0.0 & EndY = 0.0

print, "  > Enter the filepath of the .dat file specifying lines: "

OpenStatus = -1
while (OpenStatus NE 0) do begin
  read, DatFilePath
  DatFilePath = strtrim(DatFilePath,2)
  openr, lunDat, DatFilePath, /get_lun, error=OpenStatus
  if (OpenStatus NE 0) then print, "  > File cannot be opened. Re-enter filepath."
endwhile  

readf, lunDat, DatTitle, format='(a80)'
readf, lunDat, ExtraN,   format='(i6)'

ExtraLines = fltarr (ExtraN,4)

for XDat = 0, (ExtraN-1) do begin
  readf, lunDat, StartX, StartY, EndX, EndY, format='(4f10.2)'
  ExtraLines [XDat,*] = [StartX, StartY, EndX, EndY]
endfor

free_lun, lunDat

end

;******************************************************************************
; toggle between 1 and 4 plots

(Choice EQ 8): begin	

if (TogglePlots EQ 1) then begin
  TogglePlots = 4
endif else begin
  TogglePlots = 1
endelse

end

;******************************************************************************
; alter focus on region

(Choice EQ 9): begin	

print, "  > Enter the number of degrees to plot around the region: "
read, Spread

print, "  > Enter the longitude on which to centre the projection: "
read, P0Long

end

;******************************************************************************
; refresh screen or print to .eps

(Choice EQ 10 OR Choice EQ 40 OR Choice EQ 41): begin	

print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

for XLong = 0, (ModelLongN-1) do begin				; refresh SelColour
  for XLat = 0, (ModelLatN-1) do begin
    if (SelArray(XLong,XLat) NE MissVal) then begin
      SelColour (XLong,XLat) = RegColour (SelArray(XLong,XLat))
    endif else begin
      SelColour (XLong,XLat) = Grey
    endelse    
  endfor
endfor

if (Choice EQ 10) then begin					; refresh RegBounds
  GetRegBounds, SelArray, CurrentReg, ModelLongN, ModelLatN, ModelChosen, Spread, CurrentRegBounds
  RegBounds (CurrentReg,*) = CurrentRegBounds (*)
  ThisDisplay = PlotDisplay
  ThisView    = [ViewChosen,ViewChosen,ViewChosen,ViewChosen]
endif else begin				; get views to plot
  print, "  > Enter the four views to plot: "
  read, View0, View1, View2, View3
  ThisDisplay = [0,0,0,0]
  ThisView    = [View0, View1, View2, View3]
endelse

if (Choice EQ 10) then begin
  set_plot, 'X'	
endif else begin
  PSPath = ""
  print, "  > Enter the filepath of the .eps file: "
  read, PSPath

  Set_Plot, 'ps', /copy
  device, filename=PSPath, bits_per_pixel=8, /Color, /encapsulated
endelse

if (TogglePlots EQ 4) then begin
  !P.Multi = [0, 2, 2, 0, 0]					; initialise screen
endif else begin
  !P.Multi = [0, 1, 1, 0, 0]
endelse


for XPlot = (4-TogglePlots), 3 do begin
  if (TogglePlots EQ 4) then begin
    !P.Multi = [(4-XPlot), 2, 2, 0, 0]
  endif else begin
    !P.Multi = [1, 1, 1, 0, 0]
  endelse

  if (Choice EQ 10) then begin
      Set_Plot, 'X'
  endif else begin
      Set_Plot, 'ps', /copy
  endelse
  							; load colour table
  if (ThisDisplay(XPlot) EQ 2 OR ThisDisplay(XPlot) EQ 3) then begin 
    LoadAnyCT, ScaleColTab (ScaleChosen)
  endif else begin
    loadct, 38, /silent
    tvlct,   0,   0,   0,   0		; black
    tvlct, 200, 200, 200, 254		;  grey
    tvlct, 255, 255, 255, 255		; white
  endelse
  							; specify bounds
  if (ThisDisplay(XPlot) EQ 0 OR ThisDisplay(XPlot) EQ 2 OR RegBounds[CurrentReg,0] EQ MissVal) then begin
    LongMin = ViewBounds[ThisView(XPlot),0]
    LongMax = ViewBounds[ThisView(XPlot),1]
    LatMin  = ViewBounds[ThisView(XPlot),2]
    LatMax  = ViewBounds[ThisView(XPlot),3]
  endif else begin
    LongMin = RegBounds [CurrentReg,0]
    LongMax = RegBounds [CurrentReg,1]
    LatMin  = RegBounds [CurrentReg,2]
    LatMax  = RegBounds [CurrentReg,3]
  endelse  
  							; apply projection details
  if (ThisDisplay(XPlot) EQ 0 OR ThisDisplay(XPlot) EQ 2) then begin 
    if            (ViewKeyWords[ThisView(XPlot),0] EQ 'Cyli') then begin
      Map_set, 0, P0Long, /cylindrical, color=White, /noerase, /advance, limit=[LatMin,LongMin,LatMax,LongMax] 
    endif else if (ViewKeyWords[ThisView(XPlot),0] EQ 'Lamb') then begin
      Map_set, 0, P0Long, /lambert,     color=White, /noerase, /advance, limit=[LatMin,LongMin,LatMax,LongMax]
    endif else if (ViewKeyWords[ThisView(XPlot),0] EQ 'Hamm') then begin
      Map_set, 0, P0Long, /hammer,      color=White, /noerase, /advance, limit=[LatMin,LongMin,LatMax,LongMax]
    endif
  endif else begin
      Map_set, 0, P0Long, /cylindrical, color=White, /noerase, /advance, limit=[LatMin,LongMin,LatMax,LongMax] 
  endelse
    							; map colour array
  if (ThisDisplay(XPlot) EQ 0 OR ThisDisplay(XPlot) EQ 1) then begin
    Map = map_image (SelColour,  startx, starty, xs, ys, $
	lonmin=ImageLongMin, lonmax=ImageLongMax, $
	latmin=ImageLatMin,  latmax=ImageLatMax,  $
	compress=1, min_value=(MissVal+1), missing=Grey)
  endif else begin
    Map = map_image (InfoColour, startx, starty, xs, ys, $
	lonmin=ImageLongMin, lonmax=ImageLongMax, $
	latmin=ImageLatMin,  latmax=ImageLatMax,  $
	compress=1, min_value=(MissVal+1), missing=Grey)
  endelse
    	
  TV, Map, startx, starty, xsize=xs, ysize=ys		; televise map
							; overlay continents and countries
  if (Outlines[0] EQ 1) then Map_continents, /coasts,    color=Black, /noerase
  if (Outlines[1] EQ 1) then Map_continents, /countries, color=Black, /noerase
  if (Outlines[2] GE 1) then PlotPoly, Outlines[2], linestyle=0, color=Black
  
  for XExtra = 0, (ExtraN-1) do begin
    if (ExtraLines[XExtra,0] NE MissVal) then begin
       DrawTrueLine, ExtraLines[XExtra,*], LineStyle=0, Color=Black

;      plots, [ExtraLines[XExtra,0],ExtraLines[XExtra,2]], [ExtraLines[XExtra,1],ExtraLines[XExtra,3]], $
;      		linestyle=0, color=Black
    endif
  endfor
  
endfor

LastPlotBounds = [LongMin,LongMax,LatMin,LatMax]

if (Choice EQ 40) then 	device, /close

set_plot, 'X'	

if (SelInProg EQ 30) then Choice = 30		; return to mouse operations
if (SelInProg EQ 31) then Choice = 31
if (SelInProg EQ 32) then Choice = 32
if (SelInProg EQ 34) then Choice = 34

end

;******************************************************************************
; specify new region

(Choice EQ 11): begin	

RegToFill    = FirstFreeReg
NewRegion, RegToFill, RegToFillNext, RegToMakeCurrent, RegName, RegColour, RegSize, RegBounds
FirstFreeReg = RegToFillNext
CurrentReg   = RegToMakeCurrent

print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

end

;******************************************************************************
; select region

(Choice EQ 12): begin	

CurrentReg = SelectReg (RegName,ModelBoxN)

print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

end

;******************************************************************************
; rename region

(Choice EQ 13): begin	

CurrentReg = SelectReg (RegName,ModelBoxN)

print, "  > Enter name of region: "			; choose name
DatumCha = ""
while (DatumCha EQ "") do begin
  read, DatumCha
endwhile
RegName (CurrentReg) = DatumCha

print, "  > Enter colour of region: "			; choose colour
DatumInt = -1
while (DatumInt LT 0 OR DatumInt GT 255) do begin
  read, DatumInt
endwhile
RegColour (CurrentReg) = DatumInt

print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

end

;******************************************************************************
; delete region

(Choice EQ 14): begin	

ThisReg = SelectReg (RegName,ModelBoxN)

print, "  > deleted: ", RegName(ThisReg), " =", ThisReg, ", size = ", RegSize(ThisReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

RegName (ThisReg)     = ""					; initialise region arrays
RegSize (ThisReg)     = 0
RegColour (ThisReg)   = Grey
RegBounds (ThisReg,*) = MissVal
        
for XLong = 0, (ModelLongN-1) do begin				; set boxes in region to MissVal
  for XLat = 0, (ModelLatN-1) do begin
    if (SelArray(XLong,XLat) EQ ThisReg) then SelArray(XLong,XLat) = MissVal
  endfor
endfor
        
if (FirstFreeReg GT ThisReg) then FirstFreeReg = ThisReg	; reset FirstFreeReg
if (CurrentReg   EQ ThisReg) then begin				; reset CurrentReg
          DatumCha = ""
          XReg     = -1L
          while (DatumCha EQ "") do begin
            XReg = XReg + 1
            DatumCha = RegName(XReg)
            if (XReg EQ (ModelBoxN-1)) then DatumCha = "allblank"
          endwhile
          
          if (DatumCha EQ "allblank") then CurrentReg = 0
          if (DatumCha NE "allblank") then CurrentReg = XReg
endif
        
print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

end

;******************************************************************************
; weed out blanks

(Choice EQ 15): begin	

WeedBlanks, ModelLongN, ModelLatN, Grey, SelArray, RegSize, RegName, RegColour, FirstFreeReg
RegBounds [*,*] = MissVal
CurrentReg = FirstFreeReg - 1
print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

end

;******************************************************************************
; list regions

(Choice EQ 16): begin	

LogFilePath = '/cru/u2/f709762/data/scratch/log-mak.dat'
openw, lunLog, LogFilePath, /get_lun, error=OpenStatus

for XReg = 0L, (ModelBoxN-1) do begin
  if (RegName(XReg) NE "") then begin
    	print, "  > ", RegName(XReg), " =", XReg, ", size = ", RegSize(XReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'
    	printf, lunLog, RegName(XReg), XReg, RegSize(XReg), RegColour(XReg), format='(a20,i4,i4,i3)'
  endif
endfor

free_lun, lunLog

end

;******************************************************************************
; find unfilled

(Choice EQ 17): begin	

TotUnfilled = 0
for XLat = 0, (ModelLatN-1) do begin
  for XLong = 0, (ModelLongN-1) do begin
    if (SelArray(XLong,XLat) EQ MissVal AND InfoArray(XLong,XLat) NE MissVal) then begin
    	print, "  > ", XLong, XLat, format='(a4,i4,i4)'
    	TotUnfilled = TotUnfilled + 1
    endif
  endfor
endfor

print, "  > Search for .ref=missing & .glo=valid: total unfilled: ", TotUnfilled, format='(a58,i8)'

end

;******************************************************************************
; mask .ref by .glo

(Choice EQ 18): begin	

for XLong = 0, (ModelLongN-1) do begin
  for XLat = 0, (ModelLatN-1) do begin
    if (InfoArray(XLong,XLat) EQ MissVal AND SelArray(XLong,XLat) NE MissVal) then begin
    	RegSize (SelArray(XLong,XLat)) = RegSize (SelArray(XLong,XLat)) - 1

    	if (RegSize (SelArray(XLong,XLat)) EQ 0) then begin
          RegName   (SelArray(XLong,XLat)) = ""
          RegColour (SelArray(XLong,XLat)) = Grey
        endif

    	SelArray(XLong,XLat) = MissVal
    endif
  endfor
endfor

WeedBlanks, ModelLongN, ModelLatN, Grey, SelArray, RegSize, RegName, RegColour, FirstFreeReg
RegBounds [*,*] = MissVal
CurrentReg = FirstFreeReg - 1
print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

end

;******************************************************************************
; make land-sea
; combines all valid/non-valid boxes into two regions: land and sea

(Choice EQ 19): begin	

print, "  > valid/non-valid --> land/sea (=1); valid/non-valid --> sea/land (=2)"
read, QFlip

ValidReg = FirstFreeReg					; state valid region

XReg     = FirstFreeReg
DatumCha = "gotcha"
while (DatumCha NE "") do begin				; find non-valid region
  XReg     = XReg + 1
  DatumCha = RegName (XReg)
endwhile
NoValReg = XReg

if (QFlip EQ 1) then begin
    ValidName = "land"
    NoValName = "sea"
endif else begin
    ValidName = "sea"
    NoValName = "land"
endelse

RegName   (ValidReg) = ValidName
RegName   (NoValReg) = NoValName

RegColour (ValidReg) = 2
RegColour (NoValReg) = 242

RegSize   (ValidReg) = 0L
RegSize   (NoValReg) = 0L

RegBounds [ValidReg,*] = MissVal
RegBounds [NoValReg,*] = MissVal

for XLong = 0, (ModelLongN-1) do begin
  for XLat = 0, (ModelLatN-1) do begin
    if (SelArray(XLong,XLat) NE MissVal) then begin
      RegSize(SelArray(XLong,XLat)) = RegSize(SelArray(XLong,XLat)) - 1
      SelArray(XLong,XLat)          = ValidReg
      RegSize(ValidReg)             = RegSize(ValidReg) + 1
    endif else begin
      SelArray(XLong,XLat)          = NoValReg
      RegSize(NoValReg)             = RegSize(NoValReg) + 1
    endelse
  endfor
endfor

WeedBlanks, ModelLongN, ModelLatN, Grey, SelArray, RegSize, RegName, RegColour, FirstFreeReg
RegBounds [*,*] = MissVal
CurrentReg = FirstFreeReg - 1
print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

end

;******************************************************************************
; save .ref file

(Choice EQ 20 OR Choice EQ 21): begin	

if (Choice EQ 21) then begin
  Weightsfile=""
  print, "  > Enter the rpb .glo file with the appropriate weights (""=1.0): "
  read, WeightsFile
  Weights=fltarr(ModelLatN)
  
  if (WeightsFile EQ "") then begin
    Weights(*)=1.0  
  endif else begin
    LoadGlo, ModelChosen,GridWeights,WeightsFile,WeightsInfo
    Weights(*)=GridWeights(0,*)
  endelse
endif

WeedBlanks, ModelLongN, ModelLatN, Grey, SelArray, RegSize, RegName, RegColour, FirstFreeReg

RefRegN = FirstFreeReg
RefSize = intarr (RefRegN)
RefName = strarr (RefRegN)

for XReg = 0L, (RefRegN-1) do begin
  RefSize (XReg) = RegSize (XReg)
  RefName (XReg) = RegName (XReg)
endfor

if (Choice EQ 20) then SaveRef, RefRegN,ModelLongN,ModelLatN,SelArray,RefSize,RefName
if (Choice EQ 21) then SaveSpec,RefRegN,ModelLongN,ModelLatN,SelArray,RegName,Weights

CurrentReg = FirstFreeReg - 1
print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'
		
end

;******************************************************************************
; save .ref file as .glo

(Choice EQ 22): begin	

WeedBlanks, ModelLongN, ModelLatN, Grey, SelArray, RegSize, RegName, RegColour, FirstFreeReg

RefRegN = FirstFreeReg
RefCol  = intarr (RefRegN)

for XReg = 0L, (RefRegN-1) do begin
  RefCol  (XReg) = RegColour (XReg)
endfor

SaveGlo,ModelChosen,SelArray
; SaveRefGlo, RefRegN, ModelLongN, ModelLatN, SelArray, RefCol

CurrentReg = FirstFreeReg - 1
print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'
		
end

;******************************************************************************
; select boxes

(Choice EQ 30): begin	

SelInProg = 30
Button = 1
while (Button EQ 1) do begin
  UseMouse, LastPlotBounds, ModelChosen, ThisLong, ThisLat, Button, /Loud
  
  if (Button EQ 1) then begin
   if            (SelArray(ThisLong,ThisLat) EQ MissVal) then begin
    SelArray(ThisLong,ThisLat) = CurrentReg
    RegSize (CurrentReg)       = RegSize (CurrentReg) + 1
;    print, "  > added:   ", ThisLong,ThisLat, format='(a13,i4,i4)'
   endif else if (SelArray(ThisLong,ThisLat) EQ CurrentReg) then begin
    SelArray(ThisLong,ThisLat) = MissVal
    RegSize (CurrentReg)       = RegSize (CurrentReg) - 1
;    print, "  > deleted: ", ThisLong,ThisLat, format='(a13,i4,i4)'
   endif
  endif
  
  if (Button EQ 2) then begin			; returns to main menu
    Button = 0
    print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'
    SelInProg = 0
  endif
  
  if (Button EQ 4) then Choice = 10		; refresh screen
endwhile

end

;******************************************************************************
; select rectangles
; select by left-clicking on two opposite corners of the rectangle in turn
; confirm by right-clicking (any other choice makes selection void)

(Choice EQ 31): begin	

SelInProg = 31
Button = 1
while (Button EQ 1) do begin
  print, "  > Left = select corner ; Right = refresh ; Centre = exit"
  UseMouse, LastPlotBounds, ModelChosen, ALong, ALat, Button, /Loud
  
  if (Button EQ 1) then begin
    print, "  > Left = select opposite corner ; other = restart"
    UseMouse, LastPlotBounds, ModelChosen, BLong, BLat, Button, /Loud
    
    if (Button EQ 1) then begin
      print, "  > Left = confirm ; other = restart"
      UseMouse, LastPlotBounds, ModelChosen, CLong, CLat, Button
      
      if (Button EQ 1) then begin        
        MinLong = min ([ALong,BLong])
        MaxLong = max ([ALong,BLong])
        MinLat  = min ([ALat,BLat])
        MaxLat  = max ([ALat,BLat])
        TotPlus = 0
        TotMinus= 0
        
        for XLong = MinLong, MaxLong do begin
          for XLat = MinLat, MaxLat do begin
            if            (SelArray(XLong,XLat) EQ MissVal) then begin
    		SelArray(XLong,XLat) = CurrentReg
    		TotPlus              = TotPlus + 1
            endif else if (SelArray(XLong,XLat) EQ CurrentReg) then begin
    		SelArray(XLong,XLat) = MissVal
    		TotMinus             = TotMinus + 1
            endif
          endfor
        endfor
    	
    	RegSize (CurrentReg) = RegSize (CurrentReg) + TotPlus - TotMinus
        print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'
      endif else begin
        Button = 1
      endelse
    endif else begin
      Button = 1
    endelse
  endif
  
  if (Button EQ 2) then begin			; returns to main menu
    Button = 0
    print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'
    SelInProg = 0
  endif
  
  if (Button EQ 4) then Choice = 10		; refresh screen
endwhile

end

;******************************************************************************
; deselect all but rectangles
; select by left-clicking on two opposite corners of the rectangle in turn
; confirm by right-clicking (any other choice makes selection void)

(Choice EQ 32): begin	

SelInProg = 32
Button = 1
while (Button EQ 1) do begin
  print, "  > Left = select corner ; Right = refresh ; Centre = exit"
  UseMouse, LastPlotBounds, ModelChosen, ALong, ALat, Button
  
  if (Button EQ 1) then begin
    print, "  > Left = select opposite corner ; other = restart"
    UseMouse, LastPlotBounds, ModelChosen, BLong, BLat, Button
    
    if (Button EQ 1) then begin
      print, "  > Left = confirm ; other = restart"
      UseMouse, LastPlotBounds, ModelChosen, CLong, CLat, Button
      
      if (Button EQ 1) then begin        
        MinLong = min ([ALong,BLong])
        MaxLong = max ([ALong,BLong])
        MinLat  = min ([ALat,BLat])
        MaxLat  = max ([ALat,BLat])
        TotMinus= 0
        
        for XLong = 0, (ModelLongN-1) do begin
          for XLat = 0, (ModelLatN-1) do begin
            Wanted = 0
            if (XLong GE MinLong AND XLong LE MaxLong) then Wanted = Wanted + 1
            if (XLat GE MinLat AND XLat LE MaxLat)     then Wanted = Wanted + 1
            if (Wanted LT 2) then begin
    		if (SelArray(XLong,XLat) NE MissVal) then $
    				RegSize(SelArray(XLong,XLat)) = RegSize(SelArray(XLong,XLat)) - 1
    		SelArray(XLong,XLat) = MissVal
    		TotMinus             = TotMinus + 1
            endif
          endfor
        endfor
          
        print, "  > deleted: ", TotMinus, format='(a13,i6)'
      endif else begin
        Button = 1
      endelse
    endif else begin
      Button = 1
    endelse
  endif
  
  if (Button EQ 2) then begin			; returns to main menu
    Button = 0
    print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'
    SelInProg = 0
  endif
  
  if (Button EQ 4) then Choice = 10		; refresh screen
endwhile

end

;******************************************************************************
; make specified long/lat a region

(Choice EQ 33): begin	

print, "  > Specify box (long,lat = integers) to make into a region: "
read, ALong, ALat

RegToFill    = FirstFreeReg
NewRegion, RegToFill, RegToFillNext, RegToMakeCurrent, RegName, RegColour, RegSize, RegBounds
SelArray (ALong,ALat) = RegToFill
RegSize  (RegToFill)  = 1

FirstFreeReg = RegToFillNext
CurrentReg   = RegToMakeCurrent

print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

end

;******************************************************************************
; deselect all but rectangles
; select by left-clicking on two opposite corners of the rectangle in turn
; confirm by right-clicking (any other choice makes selection void)

(Choice EQ 34): begin	

SelInProg = 34
OpChoice  = 1
Button    = 0
while (OpChoice EQ 1) do begin
  OpChoice = 0
  print, "  > Select rectangle (=1), refresh screen (=2), or exit (=3): "
  while (OpChoice LT 1 OR OpChoice GT 3) do begin
    read, OpChoice
  endwhile

  if            (OpChoice EQ 1) then begin
    ManualMouse, LastPlotBounds, ModelChosen, ALong, ALat
    print,"  > x,y ref: ",Along,ALat
    
    OpChoice = 0
    print, "  > Continue selection (=1), or reselect (=4): "
    while (OpChoice NE 1 AND OpChoice NE 4) do begin
      read, OpChoice
    endwhile

    if (OpChoice EQ 1) then begin
      ManualMouse, LastPlotBounds, ModelChosen, BLong, BLat
      print,"  > x,y ref: ",Blong,BLat
        
      OpChoice = 0      
      print, "  > Confirm (=1), or reselect (=4): "
      while (OpChoice NE 1 AND OpChoice NE 4) do begin
        read, OpChoice
      endwhile
      
      if (OpChoice EQ 1) then begin
        MinLong = min ([ALong,BLong])
        MaxLong = max ([ALong,BLong])
        MinLat  = min ([ALat,BLat])
        MaxLat  = max ([ALat,BLat])
        TotPlus = 0
        TotMinus= 0
        
        for XLong = MinLong, MaxLong do begin
          for XLat = MinLat, MaxLat do begin
            if            (SelArray(XLong,XLat) EQ MissVal) then begin
    		SelArray(XLong,XLat) = CurrentReg
    		TotPlus              = TotPlus + 1
            endif else if (SelArray(XLong,XLat) EQ CurrentReg) then begin
    		SelArray(XLong,XLat) = MissVal
    		TotMinus             = TotMinus + 1
            endif
          endfor
        endfor
    	
    	RegSize (CurrentReg) = RegSize (CurrentReg) + TotPlus - TotMinus
        print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'
      endif else begin
        OpChoice = 1
      endelse
    endif else begin
      OpChoice = 1
    endelse  
  endif else if (OpChoice EQ 2) then begin
    Choice = 10
    Button = 4
  endif else if (OpChoice EQ 3) then begin
    Button = 0
    print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'
    SelInProg = 0
  endif
endwhile

end

;******************************************************************************
; convert region's constituent boxes into a new region each

(Choice EQ 35): begin	

ThisReg = SelectReg (RegName,ModelBoxN)

print, "  > converting..."

for XLong = 0, (ModelLongN-1) do begin				; transform
  for XLat = 0, (ModelLatN-1) do begin
    if (SelArray(XLong,XLat) EQ ThisReg) then begin
      NewReg               = FirstFreeReg
      RegName (NewReg)     = "auto"
      RegSize (NewReg)     = 1
      RegColour (NewReg)   = 2 + (16 * (NewReg MOD 16L))
      SelArray(XLong,XLat) = NewReg
      
      while (RegName(FirstFreeReg) NE "") do begin		; update FirstFreeReg
        FirstFreeReg = FirstFreeReg + 1
      endwhile      
    endif
  endfor
endfor

print, "  > convert: ", RegName(ThisReg), " =", ThisReg, ", size = ", RegSize(ThisReg), $
		", col = ", RegColour(ThisReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

RegName   (ThisReg)   = ""					; remove chosen reg
RegSize   (ThisReg)   = 0
RegColour (ThisReg)   = Grey
RegBounds (ThisReg,*) = MissVal
        
if (RegName(CurrentReg) EQ "") then begin			; update CurrentReg
  CurrentReg = 0
  while (RegName(CurrentReg) EQ "") do begin
        CurrentReg = CurrentReg + 1
  endwhile      
endif

print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'

end

;******************************************************************************
; name single-boxes by ref

(Choice EQ 36): begin	

for XLong = 0, (ModelLongN-1) do begin				; transform
  for XLat = 0, (ModelLatN-1) do begin
    XReg=SelArray(XLong,XLat)
    if (XReg NE MissVal) then begin
      RegName (XReg) = string((XLong+1),format="(i4)") + '_' + string((XLat+1),format="(i4)")
    endif
  endfor
endfor

end

;******************************************************************************
; convert .glo to .ref
; this works by converting all .glo values to positive integers
; the old .ref is effectively eliminated
; the .glo values are used as region indices
; then all the blanks are pruned
; 7.10.03

(Choice EQ 37): begin	

print, "  > Delete existing regions ..."
for ThisReg = 0L, (RefRegN-1) do begin
 RegName (ThisReg)     = ""					; delete existing regions
 RegSize (ThisReg)     = 0
 RegColour (ThisReg)   = Grey
 RegBounds (ThisReg,*) = MissVal
endfor

print, "  > Delete allocation of boxes to regions ..."
for XLong = 0, (ModelLongN-1) do begin				; de-assign boxes
  for XLat = 0, (ModelLatN-1) do begin
    SelArray(XLong,XLat) = MissVal
  endfor
endfor

print, "  > Reassign boxes using .glo ..."
for XLong = 0, (ModelLongN-1) do begin				; re-assign boxes
  for XLat = 0, (ModelLatN-1) do begin
    if (InfoArray(XLong,XLat) NE MissVal) then begin
      ThisReg = abs(fix(InfoArray(XLong,XLat)))			; identify region index
      
      if (ThisReg LT RefRegN) then begin
       if (RegSize(ThisReg) EQ 0) then begin			; if new region
        RegName (ThisReg)     = string(ThisReg,format="(i10)")
        RegColour (ThisReg)   = Grey
       endif
      
       RegSize(ThisReg)=RegSize(ThisReg)+1			; augment region size
;       GetRegBounds,SelArray,ThisReg,ModelLongN,ModelLatN, $	; reidentify bounds
;      		ModelChosen,Spread,CurrentRegBounds
;       RegBounds (ThisReg,*) = CurrentRegBounds (*)
      
       SelArray(XLong,XLat) = ThisReg				; assign box to region
       CurrentReg=ThisReg					; make reg=current
      endif
    endif
  endfor
endfor

print, "  > Weed blanks ..."
WeedBlanks, ModelLongN, ModelLatN, Grey, SelArray, RegSize, RegName, RegColour, FirstFreeReg
RegBounds [*,*] = MissVal
CurrentReg = FirstFreeReg - 1

end

;******************************************************************************
; merge one region with another

(Choice EQ 38): begin	

DropReg=0

while (DropReg NE MissVal) do begin
  print, "  > Merge regions: A-B=reg to drop, C=reg to grow ; -999=finish"
  DropReg = SelectReg (RegName,ModelBoxN)
  ClotReg = SelectReg (RegName,ModelBoxN)
  GrowReg = SelectReg (RegName,ModelBoxN)
  
  if (DropReg NE -999 AND ClotReg NE -999 AND GrowReg NE -999) then begin
   for XReg=DropReg,ClotReg do begin
    RegName   (XReg)   = ""				; initialise region arrays
    RegSize   (XReg)   = 0
    RegColour (XReg)   = Grey
    RegBounds (XReg,*) = MissVal
        
    for XLong = 0, (ModelLongN-1) do begin		; set boxes in region to new reg
     for XLat = 0, (ModelLatN-1) do begin
      if (SelArray(XLong,XLat) EQ XReg) then begin
       SelArray(XLong,XLat) = GrowReg
       RegSize (GrowReg)    = RegSize (GrowReg) + 1
      endif
     endfor
    endfor
   endfor
  endif
endwhile

print, "  > Weed blanks ..."
WeedBlanks, ModelLongN, ModelLatN, Grey, SelArray, RegSize, RegName, RegColour, FirstFreeReg
RegBounds [*,*] = MissVal
CurrentReg = FirstFreeReg - 1

end

; *********************************************************************
; end of main loop

else: begin
	print, ""
	print, "  >  0 : load .ref file"
	print, "  >  1 : load .glo file"
	print, "  >  2 : change order of plots"
	print, "  >  3 : change view"
	print, "  >  4 : change scale"
	print, "  >  5 : reverse all scales"
	print, "  >  6 : select outlines"
	print, "  >  7 : specify lines"
	print, "  >  8 : toggle between 1,4 plots"
	print, "  >  9 : alter region plot"
	print, "  > 10 : refresh screen"
	print, "  > 11 : specify new region"
	print, "  > 12 : select region"
	print, "  > 13 : rename region"
	print, "  > 14 : delete region"
	print, "  > 15 : weed out blanks"
	print, "  > 16 : list regions"
	print, "  > 17 : find unfilled"
	print, "  > 18 : mask .ref by .glo"
	print, "  > 19 : make land-sea"	 
	print, "  > 20 : save .ref file"
	print, "  > 21 : save .spec file"
	print, "  > 22 : save .ref as .glo"
	print, "  > 30 : select boxes"
	print, "  > 31 : select rectangles"
	print, "  > 32 : delete all but rectangles"
	print, "  > 33 : make long/lat box a region"
	print, "  > 34 : select rectangles manually"
	print, "  > 35 : convert region into rpb"
	print, "  > 36 : name single-boxes by ref"
	print, "  > 37 : convert .glo to .ref"
	print, "  > 38 : merge reg with another"
	print, "  > 40 : print .ref to .eps [4]"
	print, "  > 41 : print .ref to .eps [1]"
	print, "  > 99 : end"	
end

endcase

   if (Button EQ 0) then begin						; only if not on mouse operations        
     print, ""
     print, "  > Main menu: Enter a number:"
     read, Choice
   endif
									; ensures current region is named
   if (Choice GE 10 AND RegName(CurrentReg) EQ "" AND Choice NE 99) then begin
     print, "  > Current region has no name or colour."
     RegToFill    = CurrentReg
     NewRegion, RegToFill, RegToFillNext, RegToMakeCurrent, RegName, RegColour, RegSize, RegBounds
     FirstFreeReg = RegToFillNext
     CurrentReg   = RegToMakeCurrent

     print, "  > current: ", RegName(CurrentReg), " =", CurrentReg, ", size = ", RegSize(CurrentReg), $
		", col = ", RegColour(CurrentReg), format='(a13,a20,a2,i6,a9,i6,a8,i3)'
   endif

endwhile

print, ""
!P.Multi = [0, 1, 1, 0, 0]

end