* write out a time series of an averaged none-rectangular area in ASCII.
* time and space interval has to be set in grads first
*
* USAGE: run wrseries.quer filename fld (x1/y1) (x2/y2) (x3/y3) (x4/y4) [-ovy]
* where 
*   fld:   grads expressions
*
*       x4/y4***
*        *       **x3/y3
*       *            *
*      *              *
*     *                 *
*    *                   *
*  x2/y2*****              *
*             *****         *
*                   *****     *
*                        ****x1/y1 
*
* POINT 1 is southernmost point POINT 2(west),3(east) are contected to POINT 1 
*
*
*    -ov:   option for output of the data to screen.
*    -oy:   scale time step by 1/12
*    -ovy:   both -ov and -oy
*
* NOTE: The space units must be set to an exact data point!
*
* Example: run wrseries.quer datafile.txt zeta 10 1 2 2 8 10 4 15
*-> write out a time series of an averaged area in ASCII. the area is:
*
*       4/15***
*        *      **8/10
*       *            *
*      *              *
*     *                 *
*    *                   *
*  2/2  *****              *
*             *****         *
*                   *****     *
*                        ****10/1
*
*
*
*AUTHOR: dietmar dommenget
*
*
*
*
function wrseries(arg)
file=subwrd(arg,1)
fld=subwrd(arg,2)
x1=subwrd(arg,3)
x2=subwrd(arg,5)
x3=subwrd(arg,7)
x4=subwrd(arg,9)
y1=subwrd(arg,4)
y2=subwrd(arg,6)
y3=subwrd(arg,8)
y4=subwrd(arg,10)
if(y4= '')
  say 'NOT ENOUGH KOORDINATES !!!'
  say 'USAGE: run wrseries filename fld x1 y1 x2 y2 x3 y3 x4 y4 [-ovy] '
  say 'where  '
  say '  fld:   grads expressions '
  say '* '
  say '      x4/y4*** '
  say '       *       **x3/y3 '
  say '      *            * '
  say '     *              * '
  say '    *                 * '
  say '   *                   * '
  say ' x2/y2*****              * '
  say '            *****         * '
  say '                  *****     * '
  say '                       ****x1/y1  '
  say '* '
  say '* '
  say 'POINT 1 is southernmost point POINT 2(west),3(east) are contected to POINT 1 '
  say '* '
  say '* '
  say '   -ov:   option for output of the data to screen. '
  say '   -oy:   scale time step by 1/12 '
  say '   -ovy:   both -ov and -oy '
  return
endif
vrbos=subwrd(arg,11)
yrbos=subwrd(arg,11)
vopt=''
yopt=''
if(vrbos= '-ov')
   vopt='v'
endif
if(vrbos= '-oy')
   yopt='y'
endif
if(vrbos= '-ovy')
   yopt='y'
   vopt='v'
endif
if(vrbos= '-oyv')
   yopt='y'
   vopt='v'
endif

gxstate()
say '     writing 'fld' for '_ts'<=t<='_te' to 'file'.'
say ' dim:   x1/y1: 'x1'/'y1'  x2/y2: 'x2'/'y2
say '        x3/y3: 'x3'/'y3'  x4/y4: 'x4'/'y4
*
*
yy=y1
xs1=x1
xe1=x1
if(y3 >= y4 & y3 >= y2 );ynorth=y3;endif
if(y2 >= y3 & y2 >= y4 );ynorth=y2;endif
if(y4 >= y3 & y4 >= y2 );ynorth=y4;endif
if(y2-y1 !=0)
  slopxs1=(x2-x1)/(y2-y1)
else
  xs1=x2
endif
if(y3-y1 !=0)
  slopxe1=(x3-x1)/(y3-y1)
else
  xe1=x3
endif
say '****'
'!rm -f 'file
tt=_ts
while (tt <= _te & yy <= ynorth )
*
* berechne slop fuer xe, xs
*
   if( yy = y1)
     xstart=xs1
     xende=xe1
     slopxs=slopxs1
     slopxe=slopxe1
   else
     xstart= slopxs + xs
     xende= slopxe + xe
   endif
   diff=yy-y3
   if( (diff>=0 & diff<=0.5) | (diff<0 & diff>=-0.5) )
     if( y4-y3 !=0 );slopxe=(x4-x3)/(y4-y3);endif
   endif
   diff=yy-y2
   if( (diff>=0 & diff<=0.5) | (diff<0 & diff>=-0.5) )
     if( y4-y2 !=0 );slopxs=(x4-x2)/(y4-y2);endif
   endif
   xs=xstart
   xe=xende
   'set t 'tt
   time=tt
   if(yopt='y')
     time=tt/12
   endif
   if(tt = _ts);say 'xstart='xs'  xende='xe'  y='yy;endif
*
* addiere lat-linien
*
   xline=xe-xs+1
   if( yy = y1)
     xsum=xline
     'define line='xline'*aave('fld',x='xs',x='xe',y='yy',y='yy')'
   else
     xsum=xsum+xline
     'define line=line+'xline'*aave('fld',x='xs',x='xe',y='yy',y='yy')'     
   endif
*
* t - schleife
*
   yy=yy + 1
   if(yy > ynorth)
     'd line/'xsum
     ydum=subwrd(result,4)
     if(vopt='v');say '  'fld'('time')='ydum;endif
     good=write(file,time' 'ydum,append)
     tt=tt + 1
     yy=y1
   endif
endwhile
'set t '_ts' '_te
return

******************************************************
******************************************************
function gxstate()
*
*  Get dimension information as global var _xs,_xe,...
*
'query dim'
dinf = result
lx = sublin(dinf,2)
ly = sublin(dinf,3)
lz = sublin(dinf,4)
lt = sublin(dinf,5)
if ( subwrd(lx,7) = 'to')
  _lons = subwrd(lx,6)
  _lone = subwrd(lx,8)
  _xs = subwrd(lx,11)
  _xe = subwrd(lx,13)
else
  _lons = subwrd(lx,6)
  _lone = subwrd(lx,6)
  _xs = subwrd(lx,9)
  _xe = subwrd(lx,9)
endif
if ( subwrd(ly,7) = 'to')
  _lats = subwrd(ly,6)
  _late = subwrd(ly,8)
  _ys = subwrd(ly,11)
  _ye = subwrd(ly,13)
else
  _lats = subwrd(ly,6)
  _late = subwrd(ly,6)
  _ys = subwrd(ly,9)
  _ye = subwrd(ly,9)
endif
if ( subwrd(lz,7) = 'to')
  _levs = subwrd(lz,6)
  _leve = subwrd(lz,8)
  _zs = subwrd(lz,11)
  _ze = subwrd(lz,13)
else
  _levs = subwrd(lz,6)
  _leve = subwrd(lz,6)
  _zs = subwrd(lz,9)
  _ze = subwrd(lz,9)
endif
if ( subwrd(lt,7) = 'to')
  _tims = subwrd(lt,6)
  _time = subwrd(lt,8)
  _ts = subwrd(lt,11)
  _te = subwrd(lt,13)
else
  _tims = subwrd(lt,6)
  _time = subwrd(lt,6)
  _ts = subwrd(lt,9)
  _te = subwrd(lt,9)
endif
return
******************************************************
*****************************************************