*
*  Script to allow user to interactively create a 
*  mask data set for any resolution.  Given a GrADS
*  descriptor file name, this script will open that file, 
*  and display a grid at the resolution of the file, along
*  with three buttons.  When the mode is "zoom", the "zoom"
*  button will be highlighted.  If the user clicks and 
*  drags across the grid, a 'rubber-band' rectangle will 
*  appear, defining the zoom area.  If the "rectangle" button
*  is clicked, the user may rubber-band a rectangle to toggle
*  all those grid points to the on state.  If the toggle
*  button is hit, the user may click on individual grid points
*  to toggle them on and off.  When quit is selected the script
*  ends, and writes the mask grid to an output file, which can
*  then be set up by the user to be described by a GrADS
*  descriptor file, opened, and used as a mask grid in maskout
*  operations.
*

'reinit'
'set rgb 90 100 100 100 '
'set rgb 91  50 50 50 '
'set rgb 92 200 200 200'
_fnum = 0

*  put input file name into the ofile argument

rc = ofile('d:\model.ctl')

*  Select gridpoints for the mask

'set mpdset mres'
rc = newar(1)

'set gxout fwrite'
 
*  output file name as the fwrite operand

'set fwrite d:\maskarea.dat'
  'set x 0.5 '%(_xnum.fn+0.5)
  'set y 0.5 '%(_ynum.fn+0.5)
  'set z 1'
  'set t 1'
'd maskarea'
say result

function newar(fn)
  'set dfile 'fn
  'set x 0.5 '%(_xnum.fn+0.5)
  'set y 0.5 '%(_ynum.fn+0.5)
  'set z 1'
  'set t 1'
  varnm = subwrd(_vars.fn,1)
  'define maskarea = const('varnm',-1)'
  arnam = 'maskarea'  
  'set button 1 90 91 92 7 90 92 91 6'
  'set parea 0.8 10.5 1.5 8.0'
  _mode = 1
  'q dims'
  rec = sublin(result,2)
  _lonlo = subwrd(rec,6)
  _lonhi = subwrd(rec,8)
  rec = sublin(result,3)
  _latlo = subwrd(rec,6)
  _lathi = subwrd(rec,8)

*  loop to plot the field, and wait for a 
*  point and click, then modify the field.


   'set dbuff on'
   ret = dodisp(arnam)

   while (1)
     'swap'
     'q pos'
     x = subwrd(result,3)
     y = subwrd(result,4)
     mb = subwrd(result,5)
     cl = subwrd(result,6)
     if (cl=1) 
       bn = subwrd(result,7)
     endif
     if (cl=2) 
       rn = subwrd(result,7)
       x2 = subwrd(result,8)
       y2 = subwrd(result,9)
     endif

     if (cl=1) 
       if (bn=1); return; endif;
       if (bn=10); _mode = 1; endif; 
       if (bn=11); _mode = 2; endif; 
       if (bn=12); _mode = 3; endif; 
       if (bn=2)
         'set x 0.5 '%(_xnum.fn+0.5)
         'set y 0.5 '%(_ynum.fn+0.5)
         'q dims'
          rec = sublin(result,2)
         _lonlo = subwrd(rec,6)
         _lonhi = subwrd(rec,8)
         rec = sublin(result,3)
         _latlo = subwrd(rec,6)
         _lathi = subwrd(rec,8)
       endif
     endif
     if (cl=2) 
       if (x>x2)
         x1 = x2
         x2 = x
       else
         x1 = x
       endif
       if (y>y2)
         y1 = y2
         y2 = y
       else
         y1 = y
       endif
       if (x2-x1>0.1 & y2-y1>0.1)
         if (_mode = 1)
           lon1 = _lonlo + (x1-_xlo)*(_lonhi-_lonlo)/(_xhi-_xlo)
           lat1 = _latlo + (y1-_ylo)*(_lathi-_latlo)/(_yhi-_ylo)
           lon2 = _lonlo + (x2-_xlo)*(_lonhi-_lonlo)/(_xhi-_xlo)
           lat2 = _latlo + (y2-_ylo)*(_lathi-_latlo)/(_yhi-_ylo)
           _lonlo = lon1; _lonhi = lon2
           _latlo = lat1; _lathi = lat2
         endif
         if (_mode = 2)
           i1 = _ilo + (x1-_xlo)*(_ihi-_ilo)/(_xhi-_xlo)
           j1 = _jlo + (y1-_ylo)*(_jhi-_jlo)/(_yhi-_ylo)
           i2 = _ilo + (x2-_xlo)*(_ihi-_ilo)/(_xhi-_xlo)
           j2 = _jlo + (y2-_ylo)*(_jhi-_jlo)/(_yhi-_ylo)
           i1 = intrnd(i1)
           j = intrnd(j1)
           i2 = intrnd(i2) + 1
           j2 = intrnd(j2) + 1
           while (j<j2)
             i = i1
             while (i<i2)
               'set defval 'arnam' 'i' 'j' 1'
               i = i + 1
             endwhile
             j = j + 1
           endwhile
         endif
       endif
     endif
     if (cl=0)
       if (x<_xlo | x>_xhi | y<_ylo | y>_yhi)
         say 'Out of Range'
       else
         i = _ilo + (x-_xlo)*(_ihi-_ilo)/(_xhi-_xlo)
         j = _jlo + (y-_ylo)*(_jhi-_jlo)/(_yhi-_ylo)
         i = intrnd(i)
         j = intrnd(j)
         'q defval 'arnam' 'i' 'j
         val = subwrd(result,3)
         if (val<0) 
           'set defval 'arnam' 'i' 'j' 1'
           if (i<_xmin.arnum); _xmin.arnum = i; endif;
           if (i>_xmax.arnum); _xmax.arnum = i; endif;
           if (j<_ymin.arnum); _ymin.arnum = j; endif;
           if (j>_ymax.arnum); _ymax.arnum = j; endif;
         else 
           'set defval 'arnam' 'i' 'j' -1'
         endif
       endif
     endif
     ret = dodisp(arnam)
   endwhile

   'set dbuff off'
     
      
function dodisp (arnam)
   'set line 90'
   'draw recf 0 0 11 8.5'
   'set lon '_lonlo' '_lonhi
   'set lat '_latlo' '_lathi
   'set grid off'
   'set map 15 1 6'
   'set gxout fgrid'
   'set fgvals -1 0 1 8'
   'set grads off'
   'd 'arnam
   'set gxout grid'
   'set grads off'
   'set ccolor 1'
   'd maskout('arnam','arnam'-100)'
   'draw button 1 0.5 0.3 1 0.5 Quit'
   'draw button 2 3.5 0.3 1.3 0.5 UnZoom'
   'draw button 10 6.4 0.3 1.8 0.5 Zoom'
   'draw button 11 8.2 0.3 1.8 0.5 Rectangle'
   'draw button 12 10.0 0.3 1.8 0.5 Toggle'
   if (_mode=1) 
     'redraw button 10 1'
   endif
   if (_mode=2) 
     'redraw button 11 1'
   endif
   if (_mode=3) 
     'redraw button 12 1'
   endif
   'q gxinfo'
   rec = sublin(result,3)
   _xlo = subwrd(rec,4)
   _xhi = subwrd(rec,6)
   rec = sublin(result,4)
   _ylo = subwrd(rec,4)
   _yhi = subwrd(rec,6)
   'q xy2gr '_xlo' '_ylo
   _ilo = subwrd(result,3)
   _jlo = subwrd(result,6)
   'q xy2gr '_xhi' '_yhi
   _ihi = subwrd(result,3)
   _jhi = subwrd(result,6)
   if (_mode=1 | _mode = 2) 
     'set rband 1 box '_xlo' '_ylo' '_xhi' '_yhi
   endif
return

* int function

function intrnd(stuff)

  if (stuff<0) 
    nflag = 1
    stuff = -1 * stuff
  else
    nflag = 0
  endif
  stuff = stuff + 0.5
  res = ''
  i = 1
  c = substr(stuff,i,1)
  while (c!='' & ('x'%c)!='x.') 
    res = res%c
    i = i + 1
    c = substr(stuff,i,1)
  endwhile
  if (nflag = 1 & res > 0) 
    res = -1 * res
  endif
  return res

*  Function to open a file
*  Variables filled in:  
*     _fnum  -  number of open files
*     _znum.fn  - number of levels in file
*     _tnum.fn  - number of times in file
*     _vnum.fn  - number of variables in file
*     _levs.fn  - levels (character); one long string
*     _vars.fn  - vars (character); one long string
*     _vname.fn.vn  - long var description

function ofile (fname) 

  'open 'fname
  rec = sublin(result,2)
  if (subwrd(rec,1)="Open") 
    say 'Open error on file name: 'fname
    return
  endif
  fn = subwrd(rec,8)
  _fnum = _fnum + 1
  if (fn!=_fnum) 
    say "Internal logic error on open"
    _fnum = _fnum - 1
    return
  endif
 
* Get info on file 

  'set dfile 'fn
  'q file 'fn
  res = result
  dum = sublin(res,5)
  _xnum.fn = subwrd(dum,3)
  _ynum.fn = subwrd(dum,6)
  _znum.fn = subwrd(dum,9)
  _tnum.fn = subwrd(dum,12)
  xnum.fn = subwrd(dum,3)
  ynum.fn = subwrd(dum,6)
  dum = sublin(res,6)
  _vnum.fn = subwrd(dum,5)

* Get world values for each level 

  i = 1
  _levs.fn = ''
  while (i<=_znum.fn) 
    'set z 'i
    _levs.fn = _levs.fn % ' ' % subwrd(result,4)
    i = i + 1
  endwhile

* Get each variable abbreviation

  i = 0
  _vars.fn = ''
  while (i<_vnum.fn) 
    dum = sublin(res,i+7)
    _vars.fn = _vars.fn % ' ' % subwrd(dum,1)
    j = i + 1
    pos = wrdpos(dum,4)
    _vname.fn.j = substr(dum,pos,100)
    i = i + 1
  endwhile

return