* write_ascii_data_table.gs
* by Matthias Fripp, 6 November 2009

* this script creates an ascii file with a grid of values from the current data file
* each column of the ascii file corresponds to a particular longitude (shown at the top)
* each row of the ascii file corresponds to a particular latitude (shown at the top)
* before calling this script, open a data file and use "set" commands to choose a single level and time,
* and an appropriate range of x and y values.
* you will be prompted for the name of the variable to export and the name of the file to create

*find the dimensions that are currently active
*you may need to adjust this depending on your version of grads
'query dims'
xdims=sublin(result, 2)
ydims=sublin(result, 3)
zdims=sublin(result,4)
tdims=sublin(result,5)

xmin = subwrd(xdims, 11)
xmax = subwrd(xdims, 13)
ymin = subwrd(ydims, 11)
ymax = subwrd(ydims, 13)

latmin = subwrd(ydims, 6)
latmax = subwrd(ydims, 8)
lonmin = subwrd(xdims, 6)
lonmax = subwrd(xdims, 8)

if (substr(zdims, 1, 10) != "Z is fixed" | substr(tdims, 1, 10) != "T is fixed")
  say zdims
  say tdims
  say "You must use 'set' to select a single level and time before calling this script."
  return
endif

say "Dimensions:"
say result
say "xmin: " xmin " xmax: " xmax
say "ymin: " ymin " ymax: " ymax
say "Latitude is " latmin " - " latmax
say "Longitude is " lonmin " - " lonmax
say ""

prompt  'Please give the name of the variable to export, or press "return" to quit: '
pull varname
if (varname = ""); say "quitting, no data has been written"; return; endif

prompt  'Please give the name of the file to write, or press "return" to quit: '
pull fname
if (fname = ""); say "quitting, no data has been written"; return; endif

* define a variable for the relevant data
* this forces grads to read it all in at once, which is faster
'define myvar = ' varname
'q defval myvar -1 -1'

say ""
say "Creating file " fname
say "Exporting data for variable " varname
say "Missing-data value is " subwrd(result, 3)
say ""

* make a list of all the active longitudes
outstr=""
xstep=xmin
while xstep <= xmax
  outstr = outstr " " lonmin + (lonmax-lonmin) * (xstep-xmin)/(xmax-xmin)
  xstep = xstep + 1
endwhile

* write the longitudes to the file as a header
rc = write(fname,outstr)
rc = subwrd(rc, 1)
if (rc != 0)
  say "Error writing to file " fname "; result code: " rc
  return
endif

*loop through all latitudes, writing out data
ystep = ymin
while ystep <= ymax
* find the current latitude
  outstr = latmin + (latmax-latmin) * (ystep-ymin)/(ymax-ymin)

* get a list of values for each longitude
  xstep = xmin
  while xstep <= xmax
    'q defval myvar ' xstep ' ' ystep
    outstr = outstr " " subwrd(result,3)
    xstep = xstep + 1
  endwhile

* write the values to the output file
  rc = write(fname,outstr,append)
  rc = subwrd(rc, 1)
  if (rc != 0)
    say "Error writing to file " fname
    break
  endif

  ystep = ystep + 1
endwhile

if (rc = 0)
  numpoints = (xmax-xmin+1) * (ymax-ymin+1)
  say "Successfully wrote " numpoints " values of " varname " to file " fname
endif

rc = close(fname)