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

'reinit'

*-----------------------------------------------------------------------
*  Allow external GrADS functions
*-----------------------------------------------------------------------

rc = gsfallow("on")

*-----------------------------------------------------------------------
*  Batch mode?
*-----------------------------------------------------------------------

batch_mode = 'no'

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

do_print = 'yes'

if( do_print = 'yes' )
  'enable print ta_maskout.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)
*
*  calculate area-average of "ta", excluding area (x1:x2,y1:y2)...
*
*  define limited area for averaging...
*
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
pause()

*-----------------------------------------------------------------------
*  method 1: modify technique sent to GrADS Listserv on 12/6/2007 by
*            Jayakrishnan PR <prjayakrishnan@GMAIL.COM>
*-----------------------------------------------------------------------

'set x 1'
'set y 1'
'set z 1'
'set t 1 'nt
'tas1 = asum(ta,x=1,x='nx',y=1,y='ny')'
'tas2 = asum(ta,x='x1',x='x2',y='y1',y='y2')'
'taa1 = (tas1-tas2)/'nxya

'set gxout print'
'd taa1'
say 'taa1 = (tas1-tas2)/'nxya
say result

'set gxout line'
'c'
'd taa1'
'draw title "taa1"'
pause()

*-----------------------------------------------------------------------
*  method 2: define and use "tanew" which has missing values in area (x1:x2,y1:y2)...
*-----------------------------------------------------------------------

'set x 1 'nx
'set y 1 'ny
'set z 1'
'set t 1'
'tatmp = ta'
  yy = y1
while ( yy <= y2 )
  xx = x1
while ( xx <= x2 )
  'set defval tatmp 'xx' 'yy' -1'
  xx=xx+1
endwhile
  yy=yy+1
endwhile
*
*  new temperature field, using "maskout" to set negatives to "undef" in
*  limited area...
*
'set x 1 'nx
'set y 1 'ny
'set z 1'
'set t 1 'nt
'tanew = maskout(ta,tatmp)'

'set gxout shaded'
'set x 1 'nx
'set y 1 'ny
'set z 1'
'set t 1'
'c'
'd ta(t=1)'
'draw title "ta(t=1)"'
'run cbar 0.7'
if( do_print = 'yes' )
  'print'
endif
pause()
'c'
'd tanew(t=1)'
'draw title "tanew(t=1)"'
'run cbar 0.7'
if( do_print = 'yes' )
  'print'
endif
pause()
  tt = 1
while ( tt <= nt )
  say '"tt" = 'tt
  'c'
  'd ta(t='tt')-tanew(t='tt')'
  'draw title "ta(t='tt')-tanew(t='tt')"'
  'run cbar 0.7'
**pause()
  tt=tt+1
endwhile
pause()
  tt = 1
while ( tt <= nt )
  say '"tt" = 'tt
  'c'
  'd tanew(t='tt')-tanew(t=1)'
  'draw title "tanew(t='tt')-tanew(t=1)"'
  'run cbar 0.7'
**pause()
  tt=tt+1
endwhile
pause()
'c'
'd tanew(t='nt')-tanew(t=2)'
'draw title "tanew(t='nt')-tanew(t=2)"'
'run cbar 0.7'
if( do_print = 'yes' )
  'print'
endif
pause()
say
*
*  calculate area-average of "ta" by excluding area (x1:x2,y1:y2), and for
*  total domain (two ways)...
*
'set x 1'
'set y 1'
'set z 1'
'set t 1 'nt
'taa2 = aave(tanew,x=1,x='nx',y=1,y='ny')'
'taa = aave(ta,x=1,x='nx',y=1,y='ny')'
'tas = asum(ta,x=1,x='nx',y=1,y='ny')/('nx'*'ny')'

'set gxout line'
'c'
'd taa'
'draw title "taa"'
if( do_print = 'yes' )
  'print'
endif
pause()
'c'
'd taa1'
'd taa2'
'draw title "taa1", "taa2"'
if( do_print = 'yes' )
  'print'
endif
pause()
'c'
'd (taa2-taa1)/taa'
'draw title "(taa2-taa1)/taa"'
if( do_print = 'yes' )
  'print'
endif
pause()
'c'
'd (taa-taa1)/taa'
'draw title "(taa-taa1)/taa"'
if( do_print = 'yes' )
  'print'
endif
pause()
'c'
'd (taa-tas)/taa'
'draw title "(taa-tas)/taa"'
if( do_print = 'yes' )
  'print'
endif
pause()
'c'
'd (taa-tas)/taa'
'd (taa2-taa1)/taa'
'draw title "(taa-tas)/taa", "(taa2-taa1)/taa"'
if( do_print = 'yes' )
  'print'
endif
pause()
*
*  verify for a few individual time levels... do time-level by time-level...
*
'set x 1 'nx
'set y 1 'ny
'set z 1'
'set t 1'
  tt = 1
while ( tt <= 3 )
  say '"tt" = 'tt
  'tatmp = ta(t='tt')'
    yy = y1
  while ( yy <= y2 )
    xx = x1
  while ( xx <= x2 )
    'set defval tatmp 'xx' 'yy' -1'
    xx=xx+1
  endwhile
    yy=yy+1
  endwhile
*
*  new temperature field, using "maskout" to set negatives to "undef" in
*  limited area...
*
  'set x 1 'nx
  'set y 1 'ny
  'set z 1'
  'set t 'tt
  'tanew = maskout(ta,tatmp)'
  'c'
  'd tanew(t='tt')-tatmp(t='tt')'
  'draw title "tanew(t='tt')-tatmp(t='tt')"'
  'run cbar 0.7'
pause()
  'c'
  'd ta(t='tt')-tanew(t='tt')'
  'draw title "ta(t='tt')-tanew(t='tt')"'
  'run cbar 0.7'
pause()
  'undefine tanew'
  tt=tt+1
endwhile

*-----------------------------------------------------------------------
*  method 3: use external "mask" field...
*-----------------------------------------------------------------------

'open ta_mask.ctl'
'set dfile 1'
'set x 1 'nx
'set y 1 'ny
'set z 1'
'set t 1'
'mask = mask.2'
'undefine tanew'
'set t 1 'nt
'tanew = maskout(ta,mask)'
'set x 1'
'set y 1'
'set z 1'
'set t 1 'nt
'taa3 = aave(tanew,x=1,x='nx',y=1,y='ny')'

'set gxout print'
'd taa3'
say 'taa3 = aave(tanew,x=1,x='nx',y=1,y='ny')'
say result

'set gxout line'
'c'
'd taa1'
'd taa2'
'd taa3'
'draw title "taa1", "taa2", "taa3"'
if( do_print = 'yes' )
  'print'
endif
pause()
'c'
'd (taa2-taa1)/taa'
'd (taa3-taa1)/taa'
'draw title "(taa2-taa1)/taa", "(taa3-taa1)/taa"'
if( do_print = 'yes' )
  'print'
endif
pause()
'c'
'd (taa3-taa2)/taa'
'draw title "(taa3-taa2)/taa"'
if( do_print = 'yes' )
  'print'
endif
pause()
'c'
'd (taa3-taa)/taa'
'draw title "(taa3-taa)/taa"'
if( do_print = 'yes' )
  'print'
endif
pause()

if( do_print = 'yes' )
  'disable print'
endif
say
say '************************************************'
say ' Finished with this script.'
say '************************************************'
say

if( batch_mode = 'yes' )
  quit
endif

return
