***********
* Name:    print.asciigrid
* Args:    <expression>
* Purpose: prints out the given expression given by the argument to an
*          ASCII grid in ESRI format that can be imported by ArcView/ArcInfo.
*          It will write a file called ddMMMyyyy.asc for each time step.
*          ESRI's ASCII grids are defined only on quadratic grid cells; the
*          result will be distorted if the lat- and lon-resolution of the GrADS
*          expression are not equal and if they are not linear.
*
*          This script is modified from Daniel de Castro Victoria's posting
*          on the gradsusr list:
*          http://list.cineca.it/cgi-bin/wa?A2=ind04&L=gradsusr&P=R127145&X=162D1877E38271538F
* Author:  Stefan Fronzek, e-mail: Stefan(dot)Fronzek(at)ymparisto(dot)fi
* Version: 1.0 (Thu Jun 10 10:29:46 FLEDT 2004)
***********

function main(args)

var=subwrd(args,1)
if(var = '')
say "print.asciigrid - exporting grads grids to  ASCII grid (ESRI)"
say "                  output file names: ddMMMyyyy.asc"
say "Usage: run print.asciigrid <expr>"
return
endif

say "print.asciigrid - exporting ("var") with the following dimensions:"
say
'q dims'
say result

* get the X, Y, and T dimensions
'q dims'
lx=sublin(result,2)
ly=sublin(result,3)
lt=sublin(result,5)

minlon=subwrd(lx,6)
maxlon=subwrd(lx,8)
minx=subwrd(lx,11)
maxx=subwrd(lx,13)
* don't allow fixed X or Y dimensions as the resolution would be undefined
if(subwrd(lx,3) = 'fixed')
say 'X dimension may not be fixed'
return
endif
minlat=subwrd(ly,6)
maxlat=subwrd(ly,8)
miny=subwrd(ly,11)
maxy=subwrd(ly,13)
if(subwrd(ly,3) = 'fixed')
say 'Y dimension may not be fixed'
return
endif
initime=subwrd(lt,11)
fimtime=subwrd(lt,13)
if(subwrd(lt,3) = 'fixed')
initime=subwrd(lt,9)
fimtime=initime
endif

* calculate the resolution
res=(maxlat-minlat)/(maxy-miny)

** time-step loop
ti=initime
while (ti <= fimtime)
   'set t 'ti
*  set the filename to the time: <ddMMMyyyy>.asc
   'q time'
   nomel=subwrd(result,3)
   nome=substr(nomel,4,9)
   nome=nome'.asc'

   say "Exporting "nome" ..."

*  write the header information
   lixo=write(nome,"ncols       "maxx-minx+1)
   lixo=write(nome,"nrows       "maxy-miny+1)
   lixo=write(nome,"xllcenter   "minlon)
   lixo=write(nome,"yllcenter   "minlat)
   lixo=write(nome,"cellsize    "res)
   lixo=write(nome,"NODATA_value        -9999")

*  loop Y - line or lat (variable j)
   j=maxy
   while (j>=miny)
*      *isolate a line
      'set y 'j
*      *loop X - column or lon (variable i)
      i=minx
      linha=""
      while (i<=maxx)
*         *isolate a columns
         'set x 'i
*         *get the date and write to the end of a string
         'd 'var
         temp=subwrd(result,4)
         linha=linha''temp' '
         i=i+1
      endwhile
*      *write one line of data
      lixo=write(nome,linha)
      j=j-1
   endwhile

   lixo=close(nome)
   ti=ti+1
endwhile
* reset the dimensions
'set x 'minx' 'maxx
'set y 'miny' 'maxy
'set t 'initime' 'fimtime
