;+
; NAME:
;	grd_read
; PURPOSE:
;	Very provisional reader for .grd files
; CATEGORY:
;	gen/idl/util
; CALLING SEQUENCE:
	FUNCTION grd_read, InFile, data , $
		dimension	= dimension 	, $
		norow		= norow 		, $
		errormessage= errormessage	, $
		silent		= silent
; INPUTS:
;	file		scalar; type: string
;					fully-qualified file name
; OPTIONAL INPUT PARAMETERS:
;	/dimension	is set then only the array dimension from the header
;					is returned
; OUTPUTS:
;	result		scalar; type: integer
;					always 1
;	data			array[n,m]; type: float
;					data array
;				array[2]: type: integer
;					(only if /dimension is set)
;					Size of array (i.e. values of n and m)
; OPTIONAL OUTPUT PARAMETERS:
;	error_message=error_message
;				scalar; type: string
;					always null-string
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, flt_string
; SEE ALSO:
; PROCEDURE:
; MODIFICATION HISTORY:
;	JUL-2003, Paul Hick (UCSD/CASS)
;	SEP-2003, Paul Hick (UCSD/CASS)
;		Added /dimension keyword
;	SEP-2003, Paul Hick (UCSD/CASS)
;		Added /norow keyword
;	APR-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added code to handle binary grd files in addition
;		to ascii.
;-

InitVar, dimension	, /key
InitVar, silent 	, 0
InitVar, norow		, /key

fi = (file_search(InFile))[0]

; If file found then unzip if necessary.
; If file not found try adding .gz extension.

gzipped = gunzip_file(InFile+(['','.gz'])[fi eq ''], rawfile, check=fi eq '',isgz=isgz)
CASE gzipped OF
0: IF isgz THEN fi = ''
1: fi = rawfile
ENDCASE

status = 0
nrec = -1L

IF fi EQ '' THEN BEGIN
	errormessage = 'file not found or unzip error: '+hide_env(InFile)
	goto, EARLY_RETURN
ENDIF

openr, /get_lun, iu, fi, error=i
IF i NE 0 THEN BEGIN
	errormessage = !error_state.msg	; Open error
	goto, EARLY_RETURN
ENDIF

grd_type = bytarr(4)
readu, iu, grd_type
grd_type = string(grd_type)

IF grd_type NE 'DSAA' AND grd_type NE 'DSBB' THEN BEGIN
	error_message = 'unrecognized grd type'
	goto, EARLY_RETURN
ENDIF

free_lun, iu
openr, /get_lun, iu, fi

on_ioerror, IOERROR					; Establish error handler

CASE grd_type OF
'DSAA'	: BEGIN 					; Ascii file

	rec  = ''
	rec0 = ''

	; Read 5-line header. The second line contains the image dimensions

	readf, iu, rec
	readf, iu, rec

	CASE dimension OF

	0: BEGIN

		nr = round(flt_string(rec)) 	; Pick up image dimensions

		readf, iu, rec
		readf, iu, rec
		readf, iu, rec

		readf, iu, rec0

		nx = nr[0]
		ny = nr[1]

		nb = flt_string(rec0)			; # numbers/record
		nb   = n_elements(nb)

		CASE norow OF
		0: BEGIN
			;test for grid_read dlm
			command = "data = grid_read('" + InFile + "')"
			;idl hack to run command and see if it worked without erroring out on a failure
			status = execute(command,1,1)
			CASE status OF
			; here 1 = true and 0 = false for command execution
			1: BEGIN
				; grid_read should have worked and returned an array data
				; close the file unit
				free_lun, iu
			END				
			0: BEGIN
											; # records per row of the image
				nb = ceil(float(nx)/nb)

				rec = strarr(nb*ny-1)
				readf, iu, rec				; Read remainder of image
				free_lun, iu

				; Prefix the first record. rec now contains the whole file.

				rec = [rec0,rec]

				; Process the first full row of nb records to establish
				; a format for reading the remaining rows.

				ib = 0L
				ie = nb-1
				nr = flt_string(strjoin(rec[ib:ie]), fmt=fmt)

				data = fltarr(n_elements(nr),ny, /nozero)
				data[*,0] = nr

				; Use established format for 1st row to process blocks of
				; nb records.

				FOR irec=1L,ny-1 DO BEGIN

					ib = ib+nb
					ie = ie+nb

					reads, strjoin(rec[ib:ie]), format=fmt, nr
					data[*,irec] = nr

				ENDFOR
			END
			ENDCASE
		END
		1: BEGIN

			nrec = ceil(float(nx)*ny/nb)	; # records in file

			rec = strarr(nrec-1)
			readf, iu, rec					; Read remainder of image
			free_lun, iu

			; Prefix the first record. rec now contains the whole file.

			rec = [rec0,rec]
			nr = flt_string(rec[0], fmt=fmt); 1st record establishes format

			data = fltarr(nx*ny,/nozero)
			ib = -1L
			FOR irec=0L,nrec-2 DO BEGIN
				reads, rec[irec], format=fmt, nr
				data[ib+1:ib+nb] = nr
				ib = ib+nb
			ENDFOR

			nr = flt_string(rec[nrec-1])
			data[ib+1:nx*ny-1] = nr

			data = reform(data, nx,ny, /overwrite)
		END
		ENDCASE

	END

	1: BEGIN

		free_lun, iu

		data = round(flt_string(rec)) 	; Pick up image dimensions

	END

	ENDCASE

END

'DSBB'	: BEGIN 					; Binary grd file

	dsbb = bytarr(4)
	dims = intarr(2)
	readu, iu, dsbb, dims

	CASE dimension OF

	0: BEGIN
		scale = dblarr(6)
		data  = fltarr(dims)
		readu, iu, scale, data
	END

	1: data = long(dims)

	ENDCASE

	free_lun, iu

END

ENDCASE

errormessage = ''
status = 1

on_ioerror, NULL					; Cancel error handler

EARLY_RETURN:						; Return before ioerror handler is established

IF IsType(iu, /defined) THEN free_lun, iu

CASE status OF
0: BEGIN
	IF silent LE 1 THEN message, /info, errormessage
	destroyvar, data
END
1: IF silent LE 0 THEN message, /info, hide_env(fi)
ENDCASE

IF gzipped THEN i = do_file(/delete, rawfile, silent=silent)

RETURN, status

IOERROR:							; I/O error label

free_lun, iu
on_ioerror, NULL

status =  0
destroyvar, data

errormessage = !error_state.msg
IF silent LE 1 THEN message, /info, errormessage

IF gzipped THEN i = do_file(/delete, rawfile, silent=silent)

RETURN, status  &  END