function fwrite(arg)
*
* this GrADS script writes out a "map" and also generates a control file
* v1.1 7/2004
*
if (arg = '')
   say '  run fwrite.gs file expression [title]'
   exit
else
   cmdline=arg
endif

file=subwrd(cmdline,1)
data=file '.dat'
ctl=file '.ctl'
exp=subwrd(cmdline,2)
title=subwrd(cmdline,3)

if (exp = '')
   say '  run fwrite.gs file expression [title]'
   exit
endif

if (title = '')
   title=exp
endif

* change dimensions to be on grid
'query dim'
qdim = result  
xline=sublin(qdim,2)
yline=sublin(qdim,3)
zline=sublin(qdim,4)
tline=sublin(qdim,5)
xvary=subwrd(xline,3)
yvary=subwrd(yline,3)
zvary=subwrd(zline,3)
tvary=subwrd(tline,3)

* get filename
i=1
j=1
while(i<200)
  a=substr(data, i, 1)
  if ( a = ''); break; endif
  if ( a = '/')
    j=i+1
  endif
  i=i+1
endwhile
data1=substr(data,j,i-j)
err=write(ctl,'dset ^' data1)


if (xvary != 'fixed')
   'query dim'
   qdim = result
   line=sublin(qdim,2)
   lon0=subwrd(line,6)
   lon1=subwrd(line,8)
   i0=subwrd(line,11)
   i1=subwrd(line,13)
   n=math_nint(i1-i0+1)
   i0=subwrd(xline,11)
   i1=subwrd(xline,13)
   i0=math_nint(i0)
   i1=math_nint(i1)
   'set x ' i0 ' ' i1
   'query dim'
   qdim = result
   line=sublin(qdim,2)
   lon0=subwrd(line,6)
   lon1=subwrd(line,8)
   i0=subwrd(line,11)
   i1=subwrd(line,13)
   n=math_nint(i1-i0+1)
   dx=(lon1-lon0)/(n-1)
   say 'xdef ' n ' linear ' lon0 ' ' dx
   err=write(ctl, 'xdef ' n ' linear ' lon0 ' ' dx, append)
else
   'query dim'
   qdim = result
   line=sublin(qdim,2)
   lon0=subwrd(line,6)
   say 'xdef ' 1 ' linear ' lon0 ' ' 1
   err=write(ctl, 'xdef ' 1 ' linear ' lon0 ' ' 1, append)
endif

if (yvary != 'fixed')
   'query dim'
   qdim = result  
   line=sublin(qdim,3)
   lat0=subwrd(line,6)
   lat1=subwrd(line,8)
   i0=subwrd(yline,11)
   i1=subwrd(yline,13)
   i0=math_nint(i0)
   i1=math_nint(i1)
   'set y ' i0 ' ' i1
   'query dim'
   qdim = result  
   line=sublin(qdim,3)
   lat0=subwrd(line,6)
   lat1=subwrd(line,8)
   i0=subwrd(yline,11)
   i1=subwrd(yline,13)
   i0=math_nint(i0)
   i1=subwrd(line,13)
   n=math_nint(i1-i0+1)
   dx=(lat1-lat0)/(n-1)
   say 'ydef ' n ' linear ' lat0 ' ' dx
   err=write(ctl, 'ydef ' n ' linear ' lat0 ' ' dx, append)
else
   'query dim'
   qdim = result
   line=sublin(qdim,3)
   lat0=subwrd(line,6)
   err=write(ctl, 'ydef ' 1 ' linear ' lat0 ' ' 1, append)
endif

if (zvary != 'fixed')
   say 'Z can not be varying'
   exit
else
   'query dim'
   qdim = result
   line=sublin(qdim,4)
   z=subwrd(line,6)
   say 'zdef ' z ' linear ' z ' 1'
   err=write(ctl, 'zdef ' 1 ' linear ' z ' 1', append)
endif

if (tvary != 'fixed')
   'query dim'
   qdim = result
   line=sublin(qdim,5)
   i0=subwrd(line,11)
   i1=subwrd(line,13)
   i0=math_nint(i0)
   i1=math_nint(i1)
   'set t ' i0 ' ' i1
   'query dim'
   qdim = result
   line=sublin(qdim,5)
   t0=subwrd(line,6)
   i0=subwrd(line,11)
   i1=subwrd(line,13)
   n=i1-i0+1

   'query ctlinfo'
   qctlinfo = result
   m=7
   while (m<120)
       line=sublin(qctlinfo, m)
       tdf=subwrd(line,1)
       if (tdf = 'tdef'); break; endif
       if (tdf = 'TDEF'); break; endif
       m=m+1
   endwhile
   tintv=subwrd(line,5)
   say tintv
   err=write(ctl, 'tdef ' n ' linear ' t0 ' ' tintv, append)
else
   'query dim'
   qdim = result
   line=sublin(qdim,5)
   t=subwrd(line,6)
   say 'tdef ' 1 ' linear ' t ' 1mo'
   err=write(ctl, 'tdef ' 1 ' linear ' t ' 1mo', append)
endif

err=write(ctl, 'undef 9.999e20',append)
err=write(ctl, 'options little_endian',append)
err=write(ctl, 'title ' title, append)
err=write(ctl, 'vars 1', append)
err=write(ctl, 'x 0 0 ** ' title, append)
err=write(ctl, 'endvars', append)

* save 
'q gxinfo'
line=sublin(result,1)
gxout=subwrd(line,4)
if (gxout = 'Clear')
  gxout='shaded'
endif
if (gxout = 'Line')
  gxout='contour'
endif

* write data file
'set fwrite -le ' data
'set gxout fwrite'
'display ' exp
'disable fwrite'
'set gxout ' gxout
say 'set gxout ' gxout