************************************************************************
*
*  To use:  "grads -l"
*
*  Script makes a mask for a limited area of domain...
*
*  see GrADS Listserv email from 12/11/2007 by Sachin D Ghude
*  <sachinghude@TROPMET.RES.IN>
*
************************************************************************
function ta_mask

'reinit'

*-----------------------------------------------------------------------
*  GrADS metafile output?
*-----------------------------------------------------------------------

do_print = 'yes'

if( do_print = 'yes' )
  'enable print mask.gx'
endif

*-----------------------------------------------------------------------
*  Specify GrADS ctl file
*-----------------------------------------------------------------------

ctl_file = '/archive/cjs/zetac/old/KWAJEX/old_history/SKEWT.06:30z24Jul1999_2.ctl'
'open 'ctl_file
*
*  define GrADS dataset file information...
*
'q file'
rec5 = sublin(result,5)
nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9) ; nt = subwrd(rec5,12)
*
*  define limited area for mask...
*
x1 = nx/4 ; x2 = 3*nx/4
y1 = ny/4 ; y2 = 3*ny/4
say 'x1, x2, nx = 'x1', 'x2', 'nx
say 'y1, y2, ny = 'y1', 'y2', 'ny
say 'nt = 'nt
nxy1 = nx*ny
nxy2 = (x2-x1+1)*(y2-y1+1)
nxya = nxy1-nxy2
say 'nxy1, nxy2 = 'nxy1', 'nxy2
say 'nxya = 'nxya

*-----------------------------------------------------------------------
say '...make "ta" mask'
*-----------------------------------------------------------------------
*
*  define a "mask" which has "-1" in area (x1:x2,y1:y2) and "1"
*  everywhere else (ref: http://grads.iges.org/grads/gadoc/gradfuncconst.html)
*  (use "ta" variable to start...)
*
'set x 1 'nx
'set y 1 'ny
'set z 1'
'set t 1'
'mask = const(ta,1,-a)'
  yy = y1
while ( yy <= y2 )
  xx = x1
while ( xx <= x2 )
  'set defval mask 'xx' 'yy' -1'
  xx=xx+1
endwhile
  yy=yy+1
endwhile
*
*  write to GrADS fwrite file...
*
'!rm mask.ieee'
'set gxout fwrite'
'set fwrite -le ./mask.ieee'
'd mask'
'disable fwrite'
*
*  construct a local GrADS ctl file for "mask", using model GrADS ctl file
*  "ctl_file";  references: skewt_aave2.gs,
*  http://grads.iges.org/grads/gadoc/script.html#intro (read, write)
*  /home/cjs/zetac_fms.ana/processor/nc_extract (sed)
*
'q ctlinfo'
say result
rec = sublin(result,3) ; undef_rec = rec ; rec_word1 = subwrd(rec,1) ; undef_value = subwrd(undef_rec,2)
xdef_rec1 = ''
ydef_rec1 = ''
zdef_rec1 = ''
tdef_rec1 = ''
nrec = 4
while ( rec_word1 != vars )
  rec = sublin(result,nrec) ; rec_word1 = subwrd(rec,1)
  if( rec_word1 = xdef )
    say
    nrec1_xdef = nrec
    xdef_rec1 = rec
    say '"xdef" record (nrec='nrec') is:'
    say xdef_rec1
    say
  endif
  if( rec_word1 = ydef )
    say
    nrec1_ydef = nrec
    ydef_rec1 = rec
    say '"ydef" record (nrec='nrec') is:'
    say ydef_rec1
    say
  endif
  if( rec_word1 = zdef )
    say
    nrec1_zdef = nrec
    zdef_rec1 = rec
    say '"zdef" record (nrec='nrec') is:'
    say zdef_rec1
    rec = sublin(result,nrec1_zdef) ; z1 = subwrd(rec,4) ; say 'z1 = 'z1
    say
  endif
  if( rec_word1 = tdef )
    say
    tdef_rec = rec
    say '"tdef" record (nrec='nrec') is:'
    say tdef_rec
    t0 = subwrd(tdef_rec,4)
    say 'GrADS "t0" is: 't0
    dt = subwrd(tdef_rec,5)
    say 'GrADS "dt" is: 'dt
    say
  endif
  nrec = nrec+1
endwhile
if( xdef_rec1 = '' )
  say
  say 'did not find a record word = xdef...'
  say '...exiting script'
  say
  exit
endif
if( ydef_rec1 = '' )
  say
  say 'did not find a record word = ydef...'
  say '...exiting script'
  say
  exit
endif
if( zdef_rec1 = '' )
  say
  say 'did not find a record word = zdef...'
  say '...exiting script'
  say
  exit
endif
if( tdef_rec = '' )
  say
  say 'did not find a record word = tdef...'
  say '...exiting script'
  say
  exit
endif

'! grep -n xdef 'ctl_file' | cut -f1 -d: > nl_xdef'
res = read ('nl_xdef') ; nl_xdef = sublin(res,2) ; say 'nl_xdef = 'nl_xdef
'! grep -n ydef 'ctl_file' | cut -f1 -d: > nl_ydef'
res = read ('nl_ydef') ; nl_ydef = sublin(res,2) ; say 'nl_ydef = 'nl_ydef
'! grep -n zdef 'ctl_file' | cut -f1 -d: > nl_zdef'
res = read ('nl_zdef') ; nl_zdef = sublin(res,2) ; say 'nl_zdef = 'nl_zdef

res = write ('mask.ctl','dset ^mask.ieee')
res = write ('mask.ctl','title mask')
res = write ('mask.ctl','options little_endian')
res = write ('mask.ctl',''undef_rec'')
rec = close ('mask.ctl')
'! sed -n "'nl_xdef','nl_zdef-1'"p 'ctl_file' >> mask.ctl'
res = write ('mask.ctl','zdef 1 levels 'z1'',append)
res = write ('mask.ctl','tdef 1 linear 't0' 'dt'')
res = write ('mask.ctl','vars 1')
res = write ('mask.ctl',' mask 0 99 dimensionless')
res = write ('mask.ctl','endvars')
res = close ('mask.ctl')
status = sublin(res,1)
if( status != 0 )
  say
  say '"mask.ctl" status number = 'status
  say 'exiting script...'
  say
  exit
endif
say

*-----------------------------------------------------------------------
*  plot mask...
*-----------------------------------------------------------------------

'open mask.ctl'
'set dfile 2'
'set x 1 'nx
'set y 1 'ny
'set z 1'
'set t 1'
'set gxout shaded'
'c'
'set clevs 0'
'd mask'
'draw title "mask"'
'run cbar 0.7'
if( do_print = 'yes' )
  'print'
  'disable print'
endif

'!run_gxpscw mask.gx'
'!run_ps2pdf mask.ps'
'!gzip mask.ps'

say
say '************************************************'
say ' Finished with this script.'
say '************************************************'
say

return
