source: trunk/tools/transfo/grads_read.pro

;+
;
; Reading ECHAM grads files
;
; @uses
; <pro>common</pro>
;
; @history
; EG 24/2/99
;
; @version
; $Id$
;
;-
FUNCTION grads_read, file_name, var_name $
         , NCDF_DB=ncdf_db $
         , TIME_1=time_1 $
         , TIME_2=time_2 $
         , ALL_DATA=all_data $
         , NO_KEY_SHIFT=no_key_shift $
         , LEVEL=level
;
  compile_opt idl2, strictarrsubs
;
@common
;
; inits
   IF NOT keyword_set(NCDF_DB) THEN ncdf_db = iodir
   IF NOT keyword_set(TIME_1) THEN time_1 = 1
   IF NOT keyword_set(TIME_2) THEN time_2 = time_1
   IF NOT keyword_set(LEVEL) THEN level = 1

   numgrd = 20
   z = 0.
;
; decide which subdomain of data to read
;
   IF keyword_set(ALL_DATA) THEN BEGIN
      premierx = 0
      premiery = 0
      premierz = 0
      nx = jpi
      ny = jpj
      nz = jpk
      dernierx = jpi-1
      derniery = jpj-1
      dernierz = jpk-1
   ENDIF ELSE BEGIN
      grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz
      mask = 1
   ENDELSE 

; meta data file <file_name>.ctl
   m_file = ncdf_db+file_name+'.ctl'


   jpia = LONG(grep('grep -i XDEF '+m_file, ' ', 1))
   jpja = LONG(grep('grep -i YDEF '+m_file, ' ', 1))
   
   IF jpia LE 0 OR jpja LE 0 THEN BEGIN
      print, 'error in ', m_file
      return, z
   ENDIF 
;
; find record of field
;
   var_name = STRUPCASE(var_name)

; find number of levels

   spawn, 'grep ZDEF '+m_file+' | awk '' {print $2}''', line
   nlevels = long(line[0])

; find number of time interval

   spawn, 'grep TDEF '+m_file+' | awk '' {print $2}''', line
   ntime = long(line[0])

; find number of fields

   spawn, 'grep VARS '+m_file+' | awk '' {print $2}''', line
   nvars = long(line[0])
;
; nrecf = line of <var_name> minus line of VARS
;
   nlvars = LONG(grep('grep -in "VARS " '+m_file, ':', 0))
   nlfild = LONG(grep('grep -in "'+var_name+'    " '+m_file, ':', 0))

   nrecf = nlfild - nlvars - 1

; nrecfl = index of field

   delta_t = time_2 - time_1

   nrecfl = nvars*nlevels*delta_t + nrecf*nlevels + level - 1


; open data file <file_name>.grd
   s_file = ncdf_db+file_name+'.grd'
   openr, numgrd, s_file, /get_lun
;
; compute offset
;
   tab_size = jpia*jpja*4 + 8
   offset = nrecfl*tab_size + 4
;
; read data
;
   z = assoc(numgrd, fltarr(jpia, jpja, /nozero), offset)
   IF NOT keyword_set(no_key_shift) THEN $
    z = reverse(shift(z[0], key_shift, 0), 2) ELSE $
    z = reverse(shift(z[0], 0, 0), 2)
;
   IF nlevels EQ 1 THEN BEGIN 
      print, '    Read rec. ', string(nrecfl, format = '(I3)'), ':', string(var_name, format = '(A8)'), ' dims : ', string(jpia, format = '(I4)'), string(jpja, format = '(I4)'), '  (', string(delta_t, format = '(I3)'),') - min/max : ', min(z), max(z)
   ENDIF ELSE BEGIN    
      print, '    Read rec. ', string(nrecfl, format = '(I3)'), ':', string(var_name, format = '(A8)'), ' (level ',string(level, format = '(I2.2)'),') dims : ', string(jpia, format = '(I4)'), string(jpja, format = '(I4)'), string(nlevels, format = '(I4)'), '  (', string(delta_t, format = '(I3)'), ') - min/max : ', min(z), max(z)
   ENDELSE 

   close, numgrd
   free_lun, numgrd

; Field Attributes

   field = {name: '', data: z, legend: '', units: '', origin: '', dim: 0}
   field.name = var_name
   field.origin = ncdf_db+file_name
   field.legend = grepsed('grep -i "'+var_name+'    " '+m_file, 'sed ''s/.*99\(.*\)\[.*/\1/''')
   field.units = grepsed('grep -i "'+var_name+'    " '+m_file, 'sed ''s/.*\[\(.*\)\].*/\1/''')
   field.dim = 2

;
   return, field
END