;+
; NAME:
;	read_netCDF.pro
;
; PURPOSE:
;	Read netCDF file into structure variable
;
; CATEGORY:
;	All levels of processing
;
; CALLING SEQUENCE:  
;	read_netCDF, filename, data, attributes, status, $
;		pointers_for_arrays=pointers_for_arrays, $
;		single_pointers=single_pointers, fields=fields, $
;		indices=indices, dim_to_index=dim_to_index, $
;		structure=structure, tempdir=tempdir, $
;		/halt_on_error
;
; INPUTS:
;	filename = filename for existing netCDF file
; OPTIONAL INPUTS:
;	pointers_for_arrays: if set to 1, make all fields of type array in the structure
;		that use the unlimited dimension along with additional dimensions
;		(i.e. the resulting field would normally be a multi-dimensional array)
;		into an pointer array of size unlimited.
;		If set to a string or array of strings, then make each field that
;		matches the list of strings (in the same way as the field param) into an
;		array of pointers of the same length as the last dimension in the netcdf
;		variable, even single-dimension unlimited fields.
;		If a 'fillValue' attribute exists for a variable, then make the field into
;		an array of pointers regardless of this parameter's setting.
;	single_pointers: (only affects cases where the resulting structure would have
;		1 element due to not all fields using the unlimited dimension)
;		If set to 1, then make each var that uses the unlimited dimension into a
;		single pointer to the data that the field would otherwise use.
;		E.g., if the field was going to be an array of pointers as per the
;		array_pointers param or having a fillValue, then the field
;		will be a single pointer to that array of pointers.
;		If set to a string or array of strings, then apply the change
;		to each field that matches the list of strings.
; 	fields: a string or array of strings to match a subset of the variables in the file;
;		only matching fields will be extracted from the file.  Substructures of a
;		matching field will also match (e.g. if the netcdf file has vars 'foo.bar'
;		and 'foo.baz', the string 'foo' will match both vars.)
; 	indices: index(es) into the unlimited dimension of structure elements to retrieve
;	dim_to_index: name of a dimension to treat as unlimited (for indices, etc.)
;	structure: a (non-array) structure, which can contain nested
;		structures, that matches the fields of the netcdf file.
;		If datatypes are different, a netcdf field will
;		automatically be converted to the datatype of the 
;		structure field.  It will even account for pointers in
;		a structure field.
;	tempdir: if set, should be the name of an existing directory;
;		read_netCDF will copy the file to tempdir, read the
;		copy, and then delete it.
;	halt_on_error: if this flag is set, when read_netCDF is unable to populate
;		a field (due to type conversion or array sizing conflicts)
;		it will halt to let the user check the data structure and value.
;
; OUTPUTS:  
;	data = structure variable for data read from netCDF file
;	attributes = array of strings of the attributes from the netCDF file
;	status = result status: 0 = OK_STATUS, -1 = BAD_PARAMS, -2 = BAD_FILE,
;			-3 = BAD_FILE_DATA, -4 = FILE_ALREADY_OPENED,
;			-5 = TEMPDIR_FAILURE
;
; COMMON BLOCKS:
;	None
;
; PROCEDURE:
;	Check for valid input parameters
;	Open the netCDF file
;	Read the attributes into a string array
;	Create structures based on the netCDF definitions
;	Once structures are defined, then read the netCDF variables into the structure's data
;	Close the netCDF file
;
;	NetCDF IDL Procedures / Process:
;	1.	NCDF_OPEN: Open an existing netCDF file.
;	2.	NCDF_INQUIRE: Call this function to find the format of the netCDF file.
;	3.	NCDF_DIMINQ: Retrieve the names and sizes of dimensions in the file.
;	4.	NCDF_VARINQ: Retrieve the names, types, and sizes of variables in the file.
;	5.	NCDF_ATTINQ: Optionally, retrieve the types and lengths of attributes.
;	6.	NCDF_ATTNAME: Optionally, retrieve attribute names.
;	7.	NCDF_ATTGET: Optionally, retrieve the attributes.
;	8.	NCDF_VARGET: Read the data from the variables.
;	9.	NCDF_CLOSE: Close the file.
;
; USAGE EXAMPLES:
;	read_netcdf, 'afile.nc', data, attr, status, fields=['date', 'item.values']
;		- this returns a structure array containing tags named DATE and ITEM.VALUES,
;		including any subfields under those tags (eg. ITEM.VALUES.COLOR would
;		exist if found in the file.)  If the fields param is set but no matching
;		variables are found, it aborts with BAD_PARAMS status.
;		
;	read_netcdf, 'afile.nc', data, attr, status, indices=[3,6]
;		- this returns a structure array containing all fields from 'afile.nc'
;		but only indices 3,6 of the unlimited dimension from the file.
;		If there are just 5 elements in the file, only index 3 is returned;
;		if there are no corresponding indices in the file, it aborts
;		with BAD_PARAMS status.
;
;	read_netcdf, 'afile.nc', tmp, attr, status, fields='date'
;	read_netcdf, 'afile.nc', data, attr, status, indices=where(tmp.date gt 100 and tmp.date lt 250)
;		- this first retrieves all the dates from the file, and then
;		reads the full structure for all the indices where DATE meets
;		a given criteria
;
;	read_netcdf, 'afile.nc', tmp, attr, status, structure=mydata[0]
;		- builds a structure array using the structure of mydata,
;		and uses it to build the returned structure array.  If no
;		variables are found that match fields of mydata, then
;		abort with BAD_PARAMS status.
;
;	read_netcdf, 'afile.nc', tmp, attr, status, structure={s: ptr_new() }
;		- builds a structure array using the given structure.
;		Here tmp will be a single element structure, and the variable S will
;		be rendered as it would were there no given structure, and then
;		placed inside a single pointer.
;
;	read_netcdf, 'afile.nc', tmp, attr, status, structure={s: ptr_new( ptrarr(1) ) }
;		- builds a structure array using the given structure.
;		Here tmp will be a single element structure, and the variable S will
;		be converted into an array of pointers placed inside a single pointer.
;
; MODIFICATION HISTORY:
;	9/20/1999		Tom Woods		Original release of code, Version 1.00
;	12/3/1999		Tom Woods		Removed BYTE array conversion to STRING
;   05/23/2004  Don Woodraska  Prevents IDL reserved words by appending an
;                              underscore to tag names if necessary.
; 	1/26/2006	Karie Shipley	Support for variable-sized arrays
;		using pointers, extracting a subset of structure indices
;		and/or structure tags
;	3/13/2006	Karie Shipley	Support passing a structure to use
;		for the return structure, to avoid overhead from building
;		the return structure ad hoc
;	3/22/2006	Karie Shipley	Support 6 levels of nested structures
;
;	5/08/2006	Karie Shipley	Make indices apply to the unlimited
;		dimension, even if not all variables use it.
;
;	8/11/2006	Karie Shipley	No more nesting limits; support for
;		complex numbers, temporary directory
;
;-

pro	read_netCDF, filename, data, attributes, status, $
		pointers_for_arrays=pointers_for_arrays, $
		single_pointers=single_pointers, fields=fields, $
		indices=indices, dim_to_index=dim_to_index, $
		structure=structure, tempdir=tempdir
;
;	Supported IDL types (and their codes) for fields:
;		BYTE (1), INT (2), LONG (3), FLOAT (4), DOUBLE (5),
;		COMPLEX (6), STRING (7), DCOMPLEX(9)
;
;	netCDF does not currently support either 64-bit or unsigned integers,
;	but read_netCDF can cast a netCDF double variable to LONG64 (14),
;	with some resulting loss of precision for very large values 
;	(more than 15 digits.)  Note that this only works correctly
;	if the original value was a valid LONG64 value; otherwise
;	it will overflow.
;	Unsigned integers are stored as twos-complement signed integers,
;	and read_netCDF can cast INT to UINT (12) or LONG to ULONG (13)
;	with no loss of data; double can be cast to ULONG64 (15) as well,
;	with similar loss in precision as a cast to LONG64.
;


;
;	Generic "status" values
;
OK_STATUS = 0
BAD_PARAMS = -1
BAD_FILE = -2
BAD_FILE_DATA = -3
FILE_ALREADY_OPENED = -4
TEMPDIR_FAILURE = -5

debug_mode = 0			; set to 1 if want to debug this procedure

;
;	check for valid parameters
;
status = BAD_PARAMS
if (n_params(0) lt 1) then begin
	print, 'USAGE: read_netCDF, filename, data, attributes, status'
	return
endif
if (n_params(0) lt 2) then begin
	filename = ''
	read, 'Enter filename for the existing netCDF file : ', filename
	if (strlen(filename) lt 1) then return
endif

full_file = filename

; check optional parameters
fields_regexp = '.*'
if (size(fields, /type) ne 0) then begin
	fields = [ strupcase(fields) ]
	fields_regexp = '^(' + strjoin(fields + '[.].*','|') + ')'
endif else $
	fields = ''

ptr_fields_regexp = ' '
if size(pointers_for_arrays, /type) eq 7 then begin
	if pointers_for_arrays[0] eq '*' then begin
		ptr_fields_regexp = '.*'
		ptr_fields = ''
	endif else begin
		ptr_fields = [ strupcase(pointers_for_arrays) ]
		ptr_fields_regexp = '^(' + strjoin(ptr_fields + $
			'[.].*','|') + ')'
	endelse
endif else $
	ptr_fields = ''

single_ptr_fields_regexp = ' '
if size(single_pointers, /type) eq 7 then begin
	if single_pointers[0] eq '*' then begin
		single_ptr_fields_regexp = '.*'
		single_ptr_fields = ''
	endif else begin
		single_ptr_fields = [ strupcase(single_pointers) ]
		single_ptr_fields_regexp = '^(' + $
			strjoin(single_ptr_fields + '[.].*','|') + ')'
	endelse
endif else $
	single_ptr_fields = ''

if (size(indices, /type) ne 0) then $
	indices = [ indices ] $
else $
	indices = 0
ind_ndim = size(indices, /n_dimensions)
ind_nelem = n_elements(indices)

if size(dim_to_index, /type) eq 7 then $
	index_dim_name = dim_to_index[0] $
else $
	index_dim_name = ''
index_dim = -1
common_dim = -1

have_struct = 0
if keyword_set(structure) and size(structure, /type) eq 8 then $
	have_struct = 1

status = BAD_FILE
if not file_test(full_file, /read) then begin
	print, 'ERROR: file not found or file not readable by this process'
	return
endif

status = TEMPDIR_FAILURE
use_temp_file = 0
if keyword_set(tempdir) and (size(tempdir, /type) eq 7) then begin
	if not file_test(tempdir, /directory) then begin
		file_mkdir, tempdir
		if not file_test(tempdir, /directory) then begin
			print, 'ERROR: tempdir is not a directory'
			return
		endif
		print, "created temp directory '", tempdir, $
			"'; this directory will not be removed when", $
			' read_netCDF finishes.'
	endif else if not file_test(tempdir, /write) then begin
		print, 'ERROR: tempdir is not writeable by this process'
		return
	endif
	path_sep = '/'
	file_copy, full_file, tempdir, /overwrite, /require_directory
	file_path = strsplit(full_file, path_sep, /extract)
	file_basename = file_path[n_elements(file_path)-1]
	full_file = strjoin([ tempdir, file_basename ], path_sep)
	if not file_test(full_file, /read) then begin
		print, 'ERROR: temp file not found or unreadable by this process'
		return
	endif else $
		print, "created temp file '", full_file, $
			"', which will be removed when", $
			' read_netCDF finishes.'

	use_temp_file = 1
endif

null_ptr = ptr_new()

reserved_words=['AND','BEGIN','BREAK','CASE','COMMON','COMPILE_OPT', $
                'CONTINUE','DO','ELSE','END','ENDCASE','ENDELSE','ENDFOR', $
                'ENDIF','ENDREP','ENDSWITCH','ENDWHILE','EQ','FOR', $
                'FORWARD_FUNCTION','FUNCTION','GE','GOTO','GT','IF', $
                'INHERITS','LE','LT','MOD','NE','NOT','OF','ON_IOERROR', $
                'OR','PRO','REPEAT','SWITCH','THEN','UNTIL','WHILE','XOR']

;
;	Open the netCDF file
;	1.	NCDF_OPEN: Open an existing netCDF file.
;
status = BAD_FILE

if (debug_mode gt 0) then print, 'Opening ', filename, ' ...'
fid = NCDF_OPEN( full_file, /NOWRITE )

status = OK_STATUS

;
;	Create structures based on the netCDF definitions
;	2.	NCDF_INQUIRE: Call this function to find the format of the netCDF file.
;	3.	NCDF_DIMINQ: Retrieve the names and sizes of dimensions in the file.
;	4.	NCDF_VARINQ: Retrieve the names, types, and sizes of variables in the file.
;
error_status = 0
catch, error_status
if error_status ne 0 then begin
	catch, /cancel
	print, "ERROR: Couldn't process netCDF header:"
	print, "    ", !error_state.msg
	print, 'Aborting...'
	NCDF_CONTROL, fid, /ABORT
	status = BAD_FILE_DATA
	goto, cleanup
endif

finq = NCDF_INQUIRE( fid )		; finq /str = ndims, nvars, ngatts, recdim

;
;	get dimension definitions first
;	get unlimited dimension (finq.recdim)
;	
dim_unlimited = finq.recdim		; = -1 if undefined, otherwise index into dim array
complex_dim = -1
if ( finq.ndims gt 0 ) then begin
	dimstr = ' '
	dimsize = 0L
	dim_name = strarr( finq.ndims )
	dim_size = lonarr( finq.ndims )
	for k=0,finq.ndims-1 do begin
		NCDF_DIMINQ, fid, k, dimstr, dimsize
		dim_name[k] = dimstr
		dim_size[k] = dimsize
		if dimstr eq 'complex_number' then $
			complex_dim = k $
		else if dimstr eq index_dim_name then begin
			index_dim = k
			dim_unlimited = k
		endif
	endfor
endif

if (dim_unlimited ge 0) then $
	unlim_size = dim_size[dim_unlimited] $
else $
	unlim_size = 0

;
;	get variable definitions next
;	also determine nested structure levels, max. dimension, and command dimension value
;
;	LIMITATION: 6 dimensions allowed per variable
;	netCDF does not really define unsigned variable types
;
;	Have internal structure definition for tracking variables / structures
;		name = name from netCDF file
;		var_name = name from structure definition (last word after last '.')
;		type = data type value (same values as used by size())
;		natts = number of attributes for this variable
;		ndims = number of dimensions in "dim"
;		dim = dimension index into dim_size[]
;		nest_level = nest level of structures (number of '.' in name)
;		nesting_ptr = pointer to array of nesting data structures:
;			nest_name = structure name (nested)
;			nest_id = index to first case of structure name (nested)
;			nest_cnt = index of variable within a single structure (nested)
;			str_ptr = structure pointer (if first case of new structure)
;		unlim_dim = -1 if doesn't use unlimited dimension; otherwise the dimension #
;		last_dim_size_1 = if the last dim for the var in the file has size 1
;		ptr = data variable pointer
;		struct_dims = if the final variable will be an array of pointers,
;			this holds the dimensions for that pointer array;
;			otherwise just the dimensions of the var in the given struct
;		atts = pointer to attributes for this variable
;		make_ptr = whether this variable should be an array of pointers
;		make_single_ptr = whether this variable should be a single pointer
;		dont_contract = whether or not to contract the var into an array of pointers
;		
nest_item = { nest_name : '', nest_id : 0L, nest_cnt : 0L, str_ptr : null_ptr }
var_inq1 = { name : " ", var_name : " ", type : 0, natts : 0L, ndims : 0L, $
	ndim_arr : 0L, dim: lonarr(8), nest_level : 0, nesting_ptr : null_ptr, $
	unlim_dim : -1, last_dim_size_1 : 0, ptr : null_ptr, $
	struct_dims : lonarr(8), n_struct_dims : 0, atts : null_ptr, $
	make_ptr : 0, make_single_ptr : 0, dont_contract : 0}

var_inq = replicate( var_inq1, finq.nvars )
max_level = 0			; track max structure nest level while getting variable definitions
max_base_dim = 1		; track max base structure dimension required
has_common_dim = 1		; assume TRUE to start out, any conflict makes it FALSE
				; means all vars have a common last dim
var_without_unlimited_dim = 0	; assume FALSE to start out, any conflict makes it TRUE
				; means some var doesn't use the unlimited dim

;
;	sort out first the dimensions and attribute numbers
;	check for max. dim needed for base structure
;	and if should have base structure array (if all the same last dim)
;
first_var_last_dim = -1
for k=0, finq.nvars-1 do begin
	var_def = NCDF_VARINQ( fid, k )
	var_inq[k].name = strupcase(var_def.name)
	var_inq[k].ndims = var_def.ndims
	if (var_def.ndims gt 0) then $
		var_inq[k].dim[0:var_def.ndims-1] = var_def.dim
	; only set type if in the list of fields to retrieve (or no field list given)
	if ((fields_regexp ne '.*') and ((where(var_inq[k].name eq fields))[0] eq -1) and $
		((where((stregex(var_inq[k].name, fields_regexp))[0] ne -1))[0] $
		eq -1)) then $
		continue

	var_inq[k].natts = var_def.natts
	if (var_def.ndims gt 0) then begin
		if first_var_last_dim lt 0 then $
			first_var_last_dim = var_inq[k].dim[var_def.ndims-1]
		unlim_ind = (where(var_def.dim eq dim_unlimited))[0]
		if (unlim_ind ge 0) then $
			var_inq[k].unlim_dim = unlim_ind $
		else $
			var_without_unlimited_dim = 1
		lastdim = var_def.dim[var_def.ndims-1]
		lastdim_size = dim_size[lastdim]
		if lastdim_size eq 1 then var_inq[k].last_dim_size_1 = 1
		if (lastdim_size gt max_base_dim) then max_base_dim = lastdim_size
		if (lastdim ne first_var_last_dim) then $
			has_common_dim = 0
	endif else begin
		var_without_unlimited_dim = 1
		has_common_dim = 0
	endelse
	
	; get & initialize IDL datatype from the netCDF datatype
	case strupcase(var_def.datatype) of
		'BYTE': begin
			theType = 1	; use size() definitions for data type numbers
			end
		'CHAR': begin
			theType = 7	; expect STRING type
			if (debug_mode gt 0) then print, $
				'STRING type for ', var_def.name
			end
		'SHORT': theType = 2
		'INT': theType = 2
		'LONG': theType = 3
		'LONG64': theType = 14
		'DOUBLE': begin
			if var_def.dim[0] eq complex_dim then begin
				; dcomplex type is stored as DOUBLE
				; with an extra complex dimension
				theType = 9
				if (debug_mode gt 0) then print, $
					'DCOMPLEX type for ', var_def.name
			endif else $
				theType = 5
			end
		'FLOAT': begin
			if var_def.dim[0] eq complex_dim then begin
				; complex type is stored as FLOAT
				; with an extra complex dimension
				theType = 6
				if (debug_mode gt 0) then print, $
					'COMPLEX type for ', var_def.name
			endif else $
				theType = 4
			end
		'UNKNOWN': begin
			; UNKNOWN datatype indicates a truncated file; abort
			print, 'ERROR: found an UNKNOWN type for variable ', var_def.name
			theType = -1
			end
		;
		;	NOTE: Once netCDF supports unsigned types, this code
		;	block should be updated to recognize them.  Refer to
		;	the IDL size type codes and the netcdf datatype names.
		;
		else: begin
			print, "WARNING: read_netCDF doesn't recognize netCDF datatype '", $
				var_def.datatype, "', skipping ", var_def.name
			theType = 0
			end
	endcase
	;
	;	initialize structure variable definitions
	;
	var_inq[k].type = theType
	if theType ne 0 then begin
		var_inq[k].nesting_ptr = ptr_new( [ nest_item ] )
		(*var_inq[k].nesting_ptr)[0].nest_id = 0
	endif
endfor

; check for multiple vars with UNKNOWN type - abort if so
dum = where(var_inq.type eq -1, cnt)
if (cnt gt 1) and (var_inq[finq.nvars-1].type eq -1) then begin
	print,'multiple UNKNOWN variables: indicates a truncated file'
	print, 'Aborting...'
	NCDF_CONTROL, fid, /ABORT
	status = BAD_FILE
	goto, cleanup
endif

if (debug_mode gt 0) then begin
	print, ' '
	if (has_common_dim) then print, 'Array dimension for base structure = ', strtrim(max_base_dim, 2) $
	else print, 'Single structure element will be defined - max last dim. seen though is ', strtrim(max_base_dim, 2)
endif

if (has_common_dim eq 0) then begin
	max_base_dim = 1		;  make single-element structure only
endif else begin
	common_dim = first_var_last_dim
	if dim_unlimited lt 0 then begin
		; if no unlimited dimension for the fields,
		; use the common dim as "unlimited"
		dim_unlimited = common_dim
		unlim_size = dim_size[dim_unlimited]
	endif
endelse

;	make sure the indices fit in the range of 0..unlim_size
num_elements = max_base_dim
if (dim_unlimited ge 0) then begin
	if (size(indices, /n_dimensions) eq 1) then begin
		ind = where(indices lt unlim_size and indices ge 0, cnt)
		if cnt lt n_elements(indices) then $
			print, 'WARNING: index dimension has just ', unlim_size, $
				' elements; only the indices within the range [0:', $
				unlim_size-1, '] will be retrieved.'
		if (cnt ne 0) then $
			indices = temporary(indices[ind]) $
		else $
			indices = 0
		num_elements = (size(indices, /dimensions))[0]

		; set the unlimited dim size to the number of existing indices
		if (num_elements gt 0) then begin
			if unlim_size eq 1 then $
				indices = 0
			dim_size[dim_unlimited] = num_elements
			unlim_size = num_elements
		endif
	endif else indices = 0
	if (has_common_dim eq 0) then $
		num_elements = max_base_dim
endif else indices = 0
ind_ndim = size(indices, /n_dimensions)
ind_nelem = n_elements(indices)

if (num_elements eq 0) then begin
	print, 'ERROR: no matching indices'
	print, 'Aborting...'
	NCDF_CONTROL, fid, /ABORT
	status = BAD_PARAMS
	goto, cleanup
endif

;
;	now define "attributes" as string array and read attributes from the netCDF file
;	5.	NCDF_ATTINQ: If atts exist, retrieve the types and lengths of attributes.
;	6.	NCDF_ATTNAME: If atts exist, retrieve attribute names.
;	7.	NCDF_ATTGET: If atts exist, retrieve the attributes.
;
;	LIMITATION: limit attributes with more than 1 parameter are compressed into single string
;
CRLF = string( [ 13B, 10B ] )
num_att = 0L
;	finq.ngatts	= number of GLOBAL attributes from NCDF_INQUIRE earlier
if (finq.ngatts gt 0) then num_att = finq.ngatts + 1
for k=0, finq.nvars-1 do if (var_inq[k].natts gt 0) then num_att = num_att + var_inq[k].natts + 1

name_val = { name: '', val: '' }

if ( num_att gt 0 ) then begin
	attributes = strarr( num_att )
	acnt = 0L
	;
	;	do global attributes first
	;
	if ( finq.ngatts gt 0) then begin
		attributes[acnt] = 'GLOBAL:' ;	+ CRLF
		acnt = temporary(acnt) + 1L
		for jj=0,finq.ngatts-1 do begin
			att_name = NCDF_ATTNAME( fid, /GLOBAL, jj )
			NCDF_ATTGET, fid, /GLOBAL, att_name, att_value
			att_str = string( att_value )
			n_str = n_elements(att_str)
			if (n_str gt 1) then $
				att_str = strjoin(strtrim(temporary(att_str),2), ' ')
			attributes[acnt] = '    ' + att_name + ' = ' + att_str ; + CRLF
			acnt = temporary(acnt) + 1L
		endfor
	endif
	;
	;	now get variable attributes
	;
	for k=0, finq.nvars-1 do begin
		if (var_inq[k].type gt 0 and var_inq[k].natts gt 0) then begin
			attributes[acnt] = var_inq[k].name + ':' ;  + CRLF
			acnt = temporary(acnt) + 1L
			tmp_att = replicate(name_val, var_inq[k].natts)
			for jj=0,var_inq[k].natts-1 do begin
				att_name = NCDF_ATTNAME( fid, k, jj )
				NCDF_ATTGET, fid, k, att_name, att_value
				att_str = string( att_value )
				n_str = n_elements(att_str)
				if (n_str gt 1) then $
					att_str = strjoin(strtrim(temporary(att_str),2), ' ')
				tmp_att[jj].name = att_name
				tmp_att[jj].val = att_str
				attributes[acnt] = '    ' + att_name + ' = ' + att_str ; + CRLF
				acnt = temporary(acnt) + 1L
			endfor
			var_inq[k].atts = ptr_new(tmp_att)
		endif
	endfor
endif else begin
	attributes = "NONE"
endelse

; done catching errors in reading the netCDF header
catch, /cancel

if (debug_mode gt 0) then begin
	stop, 'Check out var_inq and dim_name, dim_size...'
endif


;
;	now define variables
;

data = 0
prev_var_ind = -1
match_cnt = 0
used_var_ind = where(var_inq.type ne 0, used_var_cnt)
if have_struct ne 0 then begin
	;
	; user provided a structure; use its data types when reading netcdf data
	;
	n_struct = n_elements(structure)
	first_var_last_dim = -1
	var_without_unlimited_dim = 0
	has_common_dim = 1
	has_single_item = 0

	; find the structure fields (w/nesting) that match the netcdf data
	if used_var_cnt gt 0 then for i=0, used_var_cnt-1 do begin
		k = used_var_ind[i]

		theStruct = structure[0]
		dotpos = 0
		lastpos = 0
		tmp_name = var_inq[k].name
		nesting_ptr = var_inq[k].nesting_ptr

		while (dotpos ge 0 and size(theStruct, /type) eq 8) do begin
			lastpos = dotpos
			dotpos = strpos( tmp_name, '.', lastpos )
			if (dotpos ge 0) then $
				newname = strmid(tmp_name, lastpos, dotpos-lastpos) $
			else $
				newname = strmid(tmp_name, lastpos)
			; avoid IDL reserved word or illegal name character conflicts
			if !version.release ge 6.0 then begin
				; need to use call_function to keep pre-IDL 6.0
				; versions from failing at compile time
				newname = call_function('idl_validname', $
					newname, /convert_all)
			endif else begin
				; earlier IDL version; need to fix names directly
				dum = where(newname eq reserved_words, test_result)
				dum = 0
				if test_result ne 0 then $
					newname = '_' + newname $
				else begin
					res_name = ''
					offset = 0
					if strmid(newname, 0, 1) eq '!' then begin
						res_name = '!'
						offset = 1
					endif
					pos = stregex(strmid(newname, offset), '[^A-Za-z0-9_\$]')
					while pos ge 0 do begin
						res_name = res_name + strmid(newname, offset, pos) + '_'
						offset = offset + pos + 1
						pos = stregex(strmid(newname, offset), '[^A-Za-z0-9_\$]')
					endwhile
					res_name = res_name + strmid(newname, offset)
					if stregex(res_name, '[0-9\$]') eq 0 then $
						res_name = '_' + res_name $
					else if stregex(res_name, '[^A-Za-z_\!]') eq 0 then $
						res_name = '_' + strmid(res_name, 1)
					newname = res_name
				endelse
			endelse
			if (dotpos ge 0) then begin
				(*nesting_ptr) = [ (*nesting_ptr), nest_item ]
				var_inq[k].nest_level = var_inq[k].nest_level + 1
				nn = var_inq[k].nest_level
				if (nn gt max_level) then max_level = nn
				(*nesting_ptr)[nn].nest_name = newname
				if (prev_var_ind lt 0) then k1=0 else k1 = prev_var_ind
				k1_nesting = *var_inq[k1].nesting_ptr
				if var_inq[k1].nest_level ge nn then $
					k1_nest_name = k1_nesting[nn].nest_name $
				else $
					k1_nest_name = ''
				if (prev_var_ind ge 0) and (k1_nest_name eq newname) then begin
					ncnt = k1_nesting[nn-1].nest_cnt
					(*nesting_ptr)[nn-1].nest_cnt = ncnt
					theStruct = temporary(theStruct.(ncnt))
				endif else begin
					spos = (where(newname eq tag_names(theStruct)))[0]
					if spos ge 0 then begin
						(*nesting_ptr)[nn-1].nest_cnt = spos
						theStruct = temporary(theStruct.(spos))
					endif else $
						theStruct = 0
				endelse
				dotpos = dotpos + 1
				k1_nesting = 0
			endif
		endwhile

		if (size(theStruct, /type) ne 8) then begin
			var_inq[k].type = 0
			continue
		endif

		tmp_name = newname
		var_inq[k].var_name = tmp_name

		spos = (where(tmp_name eq tag_names(theStruct)))[0]
		if spos ge 0 then begin
			; found a match with a structure field
			(*nesting_ptr)[var_inq[k].nest_level].nest_cnt = spos
			theStruct = temporary(theStruct.(spos))
			atype = size(theStruct, /type)
			struct_dims = size(theStruct, /dimensions)
			if n_struct gt 1 then begin
				if struct_dims[0] eq 0 then $
					struct_dims = [ n_struct ] $
				else $
					struct_dims = [ struct_dims, n_struct ]
			endif
			var_inq[k].struct_dims = struct_dims
			var_inq[k].n_struct_dims = n_elements(struct_dims)
			if (atype eq 8) then begin
				var_inq[k].type = 0
				continue
			endif else if (atype eq 6) or (atype eq 9) then begin
				; force complex type onto numeric var type
				if (where(var_inq[k].type eq [0,7,8,10,11]))[0] lt 0 then $
					var_inq[k].type = atype
			endif else if (atype eq 10) then begin
				; need to account for pointer type of given structure field
				ptr_dims = struct_dims
				if ptr_dims[0] eq 0 then begin
					var_inq[k].make_single_ptr = 1
					if theStruct ne null_ptr then begin
						; see if theStruct is a pointer to an array of pointers
						; if so, make the data into a pointer to a ptrarr
						if (size(*theStruct, /type) eq 10) and $
							(size(*theStruct, /n_dimensions) gt 0) then $
							var_inq[k].make_ptr = 1
					endif
				endif else $
					var_inq[k].make_ptr = 1
				var_inq[k].struct_dims = ptr_dims
				var_inq[k].n_struct_dims = n_elements(ptr_dims)
				if (keyword_set(pointers_for_arrays)) then begin
					if ((ptr_fields_regexp eq '.*') or $
						((where(var_inq[k].name eq ptr_fields))[0] eq -1) or $
						((where((stregex(var_inq[k].name, ptr_fields_regexp))[0] $
						ne -1))[0] eq -1)) then begin
						var_inq[k].make_ptr = 1
						var_inq[k].dont_contract = 1
					endif
				endif
				if (var_inq[k].natts gt 0) then begin
					atts = *var_inq[k].atts
					if ((where(atts.name eq 'fillValue'))[0] ge 0) then begin
						var_inq[k].make_ptr = 1
						var_inq[k].dont_contract = 0
					endif
				endif
				if (var_inq[k].ndims gt 0) then begin
					dims = var_inq[k].dim
					if var_inq[k].make_ptr and $
						(var_inq[k].struct_dims[0] eq 0) then begin
						var_inq[k].struct_dims[0] = $
							dim_size[dims[ var_inq[k].ndims-1 ] ]
						var_inq[k].n_struct_dims = 1
					endif
				endif


			endif
			match_cnt = match_cnt + 1
		endif else begin
			var_inq[k].type = 0
			continue
		endelse
		prev_var_ind = k
		theStruct = 0
		n_dims = var_inq[k].ndims
		if var_inq[k].n_struct_dims ne 0 then begin
			last_struct_dim = $
				var_inq[k].struct_dims[var_inq[k].n_struct_dims-1]
			if (var_inq[k].make_single_ptr ne 0) or $
				(n_dims eq 0) then $
				has_single_item = 1 $
			else if last_struct_dim eq $
				dim_size[ var_inq[k].dim[n_dims-1] ] then $
				has_single_item = 1
		endif
		if (n_dims gt 0) then begin
			lastdim = var_inq[k].dim[n_dims-1]
			if first_var_last_dim lt 0 then $
				first_var_last_dim = lastdim
			if (lastdim ne first_var_last_dim) then $
				has_common_dim = 0
			if (lastdim eq dim_unlimited) then $
				var_without_unlimited_dim = 1
			lastdim_size = dim_size[lastdim]
		endif else begin
			var_without_unlimited_dim = 1
			has_common_dim = 0
			has_single_item = 1
		endelse
	endfor

endif else if used_var_cnt gt 0 then for i=0, used_var_cnt-1 do begin
	k = used_var_ind[i]
	;
	;	no structure provided; determine the structure nesting
	;
    
	nesting_ptr = var_inq[k].nesting_ptr
	if (prev_var_ind lt 0) then (*nesting_ptr)[0].nest_cnt = 0 $
	else (*nesting_ptr)[0].nest_cnt = $
		(*var_inq[prev_var_ind].nesting_ptr)[0].nest_cnt + 1

	;  assume nest level 0 before looking for '.'
	;  increase nest_level for each '.' found and set new nest_name, nest_id, nest_cnt
	;
	dotpos = 0
	lastpos = 0
	tmp_name = var_inq[k].name
	while (dotpos ge 0) do begin
		lastpos = dotpos
		dotpos = strpos( tmp_name, '.', lastpos )
		if (dotpos ge 0) then $
			newname = strmid(tmp_name, lastpos, dotpos-lastpos) $
		else $
			newname = strmid(tmp_name, lastpos)
		; avoid IDL reserved word or illegal var name character conflict
		if !version.release ge 6.0 then begin
			newname = call_function('idl_validname', $
				newname, /convert_all)
		endif else begin
			dum = where(newname eq reserved_words, test_result)
			dum = 0
			if test_result ne 0 then $
				newname = '_' + newname $
			else begin
				res_name = ''
				offset = 0
				if strmid(newname, 0, 1) eq '!' then begin
					res_name = '!'
					offset = 1
				endif
				pos = stregex(strmid(newname, offset), '[^A-Za-z0-9_\$]')
				while pos ge 0 do begin
					res_name = res_name + strmid(newname, offset, pos) + '_'
					offset = offset + pos + 1
					pos = stregex(strmid(newname, offset), '[^A-Za-z0-9_\$]')
				endwhile
				res_name = res_name + strmid(newname, offset)
				if stregex(res_name, '[0-9\$]') eq 0 then $
					res_name = '_' + res_name $
				else if stregex(res_name, '[^A-Za-z_\!]') eq 0 then $
					res_name = '_' + strmid(res_name, 1)
				newname = res_name
			endelse
		endelse
		if (dotpos ge 0) then begin
			(*nesting_ptr) = [ (*nesting_ptr), nest_item ]
			var_inq[k].nest_level = var_inq[k].nest_level + 1
			nn = var_inq[k].nest_level
			if (nn gt max_level) then max_level = nn
			newname = strmid(tmp_name, lastpos, dotpos-lastpos)
			(*nesting_ptr)[nn].nest_name = newname
			if (prev_var_ind lt 0) then k1=0 else k1 = prev_var_ind
			is_new_nest = 1
			if var_inq[k1].type ne 0 then begin
				k1_nesting = *var_inq[k1].nesting_ptr
				if (prev_var_ind ge 0) and ( var_inq[k1].nest_level ge nn ) $
					and (n_elements(k1_nesting) gt nn) then begin
					if (k1_nesting[nn].nest_name eq newname) then begin
						is_new_nest = 0
						(*nesting_ptr)[nn-1].nest_cnt = $
							(*nesting_ptr)[nn-1].nest_cnt - 1
						(*nesting_ptr)[nn].nest_id = k1_nesting[nn].nest_id
						(*nesting_ptr)[nn].nest_cnt = $
							k1_nesting[nn].nest_cnt + 1
					endif
				endif
			endif
			if is_new_nest eq 1 then begin
				(*nesting_ptr)[nn].nest_id = k
				(*nesting_ptr)[nn].nest_cnt = 0
			endif
			dotpos = dotpos + 1
			k1_nesting = 0
		endif
	endwhile
	var_inq[k].var_name = newname

	match_cnt = match_cnt + 1
	prev_var_ind = k

	;	now set variable parameters
	;	uses dumb dimension rules : 
	;	 ndim_var = ndim_total - 1	for base structure being an array
	;	 ndim_var = ndim_var - 1	for string, complex definitions
	;
	ndim_array = var_inq[k].ndims
	atype = var_inq[k].type
	if (has_common_dim) then ndim_array = ndim_array - 1
	if (atype eq 6) or (atype eq 7) or (atype eq 9) then $
		ndim_array = ndim_array - 1
	if (ndim_array lt 0) then ndim_array = 0
	if (ndim_array eq 0) then begin
		if (var_inq[k].natts gt 0) then begin
			atts = *var_inq[k].atts
			if ((where(atts.name eq 'fillValue'))[0] ge 0) then begin
				var_inq[k].make_ptr = 1
				var_inq[k].dont_contract = 0
			endif
		endif
	endif else begin
		if (keyword_set(pointers_for_arrays)) then begin
			if ((not ((not has_common_dim) and (ndim_array eq 1) $
				and (dim_size[ var_inq[k].dim[0] ] eq 1)) $
				and (var_inq[k].unlim_dim ge 0) and $
				(ptr_fields_regexp eq ' '))) or $
				(ptr_fields_regexp eq '.*') or $
				((where(var_inq[k].name eq ptr_fields))[0] ge 0) or $
				((where((stregex(var_inq[k].name, $
				ptr_fields_regexp))[0] ne -1))[0] ge 0) then begin
				var_inq[k].make_ptr = 1
				var_inq[k].dont_contract = 1
			endif
		endif
		if (ndim_array le 7 and var_inq[k].natts gt 0) then begin
			atts = *var_inq[k].atts
			if ((where(atts.name eq 'fillValue'))[0] ge 0) then begin
				var_inq[k].make_ptr = 1
				var_inq[k].dont_contract = 0
			endif
		endif
		if keyword_set(single_pointers) and (not has_common_dim) then begin
			if ((var_inq[k].unlim_dim ge 0) and (single_ptr_fields_regexp eq ' ')) $
				or (single_ptr_fields_regexp eq '.*') or $
				((where(var_inq[k].name eq single_ptr_fields))[0] ge 0) or $
				((where((stregex(var_inq[k].name, single_ptr_fields_regexp))[0] $
				ne -1))[0] ge 0) then begin
				var_inq[k].make_single_ptr = 1
			endif
		endif
		if (var_inq[k].ndims gt 0) then begin
			dims = var_inq[k].dim
			if var_inq[k].make_ptr then $
				var_inq[k].struct_dims[0] = dim_size[dims[ var_inq[k].ndims-1 ] ]
		endif else begin
			ndim_array = 0
			var_inq[k].dim = 0
		endelse
		if ndim_array ge 8 then begin
			print, 'ERROR: read_netCDF can only handle 7 dimensions for arrays'
			print, 'Aborting...'
			NCDF_CONTROL, fid, /ABORT
			status = BAD_FILE_DATA
			goto, cleanup
		endif 
	endelse
	var_inq[k].ndim_arr = ndim_array
endfor

;
;	if no matching field variables were found, abort
;
if (match_cnt eq 0) then begin
	print, 'ERROR: no matching field names'
	print, 'Aborting...'
	NCDF_CONTROL, fid, /ABORT
	data = 0
	status = BAD_PARAMS
	goto, cleanup
endif

if have_struct ne 0 then begin
	;
	;	structure provided
	;	check the set of fields matching the structure fields,
	;	set up the final data structure
	;
	if has_common_dim and (has_single_item eq 0) then begin
		dim_unlimited = lastdim
		num_elements = lastdim_size
	endif
	if (has_single_item eq 1) then $
		data = structure[0] $
	else if n_struct ge num_elements then begin
		if n_struct gt num_elements then begin
			print, 'WARNING: the provided structure contains ', $
				strtrim(n_struct,2), ' elements, but the netCDF file ', $
				' would produce a structure of just ', $
				strtrim(num_elements,2), ' elements.', $
				'  read_netCDF will use the full given ', $
				'structure and populate its first ', $
				strtrim(num_elements,2), ' with the data.'
		endif
		data = structure
	endif else if (has_common_dim) and (num_elements gt 1) then begin
		print, 'WARNING: using the given structure as a basis ', $
			'for the new ', strtrim(num_elements,2), '-element ', $
			'data structure.'
		if n_struct eq 1 then begin
			print, 'The elements of the new structure will come ', $
				'from the first element of the given structure.'
		endif else begin
			print, 'The first ', strtrim(n_struct,2), $
			'elements of the new structure will come from the given ', $
			'structure, the rest will come from the first element ', $
			'of the given structure.'
		endelse
		data = replicate( structure[0], num_elements )
		if n_struct gt 1 then $
			data[0:n_struct-1] = structure
	endif

endif

if (debug_mode gt 0) then begin
	print, ' '
	nvar = n_elements( var_inq )
	print, 'Indx Lvl -- 0  1  2 ID 3   4   5--< 0  1  2 CT 3   4   5 >  NAME'
	for jj=0,nvar-1 do print, jj, var_inq[jj].nest_level, var_inq[jj].nest_id[0:5], var_inq[jj].nest_cnt[0:5], $
		var_inq[jj].name, form="(14I4,'    ',A)"
	stop, 'Check out var_inq and dim_name, dim_size...'
endif

used_var_ind = where(var_inq.type ne 0, used_var_cnt)

;
;	define structures based on var and dim definitions from netCDF file
;	using anonymous structure name with CREATE_STRUCT()
;
;	start with largest nest level and work down to zero level
;	store higher level structures as PTR (in var_inq[XX].str_ptr)
;
;	search backwards in variables for structure definitions
;   assume structure variables are grouped together
;
if (size(structure, /type) ne 8) then for nn=max_level,0,-1 do begin
	prev_var_ind = -1
	if used_var_cnt gt 0 then for i=0, used_var_cnt-1 do begin
		k = used_var_ind[i]
		;
		;	check if new structure found (same nest level as "nn" and cnt = 0)
		;	if new, then ss = CREATE_STRUCT( tag, value ) for first parameter and
		;	then ss = CREATE_STRUCT( ss, tag, value ) for other parameters
		;
		nesting = *var_inq[k].nesting_ptr
		n_nestings = n_elements(nesting)
		last_dim = -1
		ndims = var_inq[k].ndims

		if ndims gt 0 then $
			last_dim = var_inq[k].dim[ ndims ]
		if (prev_var_ind lt 0) then begin
			firstzero = 1
			if n_nestings gt nn then begin
				firstzero = nesting[nn].nest_cnt eq 0
			endif
		endif else begin
			prev_var_nesting = *var_inq[prev_var_ind].nesting_ptr
			firstzero_k = 1
			k_nest_id = 0
			if n_nestings gt nn then begin
				firstzero_k = nesting[nn].nest_cnt eq 0
				k_nest_id = nesting[nn].nest_id
			endif
			prev_nest_id = 0
			firstzero_prev = 0
			if firstzero_k eq 1 then begin
				if n_elements(prev_var_nesting) gt nn then begin
					firstzero_prev = prev_var_nesting[nn].nest_cnt ne 0
					prev_nest_id = prev_var_nesting[nn].nest_id
				endif
			endif
			firstzero = firstzero_k and (firstzero_prev or $
				(k_nest_id ne prev_nest_id))
		endelse
		if (var_inq[k].nest_level ge nn) and (firstzero) then begin
			if (nn lt var_inq[k].nest_level) then begin
				ss = CREATE_STRUCT( nesting[nn+1].nest_name, *(nesting[nn+1].str_ptr) )
				ptr_free, nesting[nn+1].str_ptr
			endif else begin
				if ( var_inq[k].make_ptr eq 1 ) then begin
					; tag element needs to be a pointer
					; account for unlim dim but not common dim
					if (var_inq[k].make_single_ptr eq 0) and (ndims gt 0) and (num_elements eq 1) then $
						ss = CREATE_STRUCT( var_inq[k].var_name, ptrarr(dim_size[last_dim]) ) $
					else $
						ss = CREATE_STRUCT( var_inq[k].var_name, null_ptr )
				endif else begin
					if ( var_inq[k].make_single_ptr eq 1 ) then $
						ss = CREATE_STRUCT( var_inq[k].var_name, null_ptr ) $
					else begin
						;
						;	initialize the data from the type & dimensions
						;
						dims = var_inq[k].dim
						atype = var_inq[k].type
						ndim_array = var_inq[k].ndim_arr
						theData = 0
						if ndim_array eq 0 then begin
							if atype eq 7 then $
								theData = '' $
							else $
								theData = fix(0, type=atype, /print)
						endif else begin
							cur_sizes = dim_size[ dims[0:ndim_array] ]
							if (atype eq 6) or (atype eq 7) or (atype eq 9) then begin
								theData = make_array(type=atype, $
									dim=cur_sizes[1:ndim_array])
							endif else begin
								theData = make_array(type=atype, $
									dim=cur_sizes[0:ndim_array-1])
							endelse
						endelse
						ss = CREATE_STRUCT( var_inq[k].var_name, theData )
						theData = 0
					endelse
				endelse
			endelse
			k1 = k
			for kk=k+1, finq.nvars-1 do begin
				theData = 0
				if (var_inq[kk].type eq 0) then $
					continue
				;if (var_inq[kk].ptr eq null_ptr) then $
				;	continue
				last_dim = -1
				k1_nesting = *var_inq[k1].nesting_ptr
				kk_nesting = *var_inq[kk].nesting_ptr
				ndims = var_inq[kk].ndims
				if ndims gt 0 then $
					last_dim = var_inq[kk].dim[ ndims ]
				kk_nest_id = 0
				kk_nest_cnt = 0
				if n_elements(kk_nesting) gt nn then begin
					kk_nest_id = kk_nesting[nn].nest_id
					kk_nest_cnt = kk_nesting[nn].nest_cnt
				endif
				k_nest_id = 0
				if n_elements(nesting) gt nn then $
					k_nest_id = nesting[nn].nest_id
				k1_nest_cnt = 0
				if n_elements(k1_nesting) gt nn then $
					k1_nest_cnt = k1_nesting[nn].nest_cnt
				if ( var_inq[kk].nest_level ge nn ) and ( kk_nest_id eq k_nest_id ) and $
						( kk_nest_cnt eq (k1_nest_cnt + 1) ) then begin
					if (nn lt var_inq[kk].nest_level) then begin
						ss = CREATE_STRUCT( ss, kk_nesting[nn+1].nest_name, *(kk_nesting[nn+1].str_ptr) )
						ptr_free, kk_nesting[nn+1].str_ptr
					endif else begin
						if ( var_inq[kk].make_ptr eq 1 ) then begin
							; tag element needs to be a pointer
							; account for unlim dim but not common dim
							if (var_inq[kk].make_single_ptr eq 0) and (ndims gt 0) and (num_elements eq 1) then $
								ss = CREATE_STRUCT( ss, var_inq[kk].var_name, ptrarr(dim_size[last_dim]) ) $
							else $
								ss = CREATE_STRUCT( ss, var_inq[kk].var_name, null_ptr )
						endif else begin
							if ( var_inq[kk].make_single_ptr eq 1 ) then $
								ss = CREATE_STRUCT( ss, var_inq[kk].var_name, null_ptr ) $
							else begin
								;
								;	initialize the data from the type & dimensions
								;
								dims = var_inq[kk].dim
								atype = var_inq[kk].type
								ndim_array = var_inq[kk].ndim_arr
								theData = 0
								if ndim_array eq 0 then begin
									if atype eq 7 then $
										theData = '' $
									else $
										theData = fix(0, type=atype, /print)
								endif else begin
									cur_sizes = dim_size[ dims[0:ndim_array] ]
									if (atype eq 6) or (atype eq 7) or (atype eq 9) then begin
										theData = make_array(type=atype, $
											dim=cur_sizes[1:ndim_array])
									endif else begin
										theData = make_array(type=atype, $
											dim=cur_sizes[0:ndim_array-1])
									endelse
								endelse
								ss = CREATE_STRUCT( ss, var_inq[kk].var_name, theData )
								theData = 0
							endelse
						endelse
					endelse
					k1 = kk
					k1_nesting = 0
					kk_nesting = 0
				endif
			endfor
			theData = 0
			;
			;	store new structure as PTR
			;	if BASE structure, then replicate for all data reading later
			(*var_inq[k].nesting_ptr)[nn].str_ptr = PTR_NEW( ss )
			if (nn eq 0) then begin
				data = replicate( ss, num_elements )
			endif
			if (debug_mode gt 0) then begin
				if (nn gt 0) then print, k, nn, '  Structure defined for ', nesting[nn].nest_name $
				else print, k, nn, '  Base Structure defined as '
				help, ss, /struct
			endif
		endif
		prev_var_ind = k
		prev_var_nesting = 0
		nesting = 0
	endfor
endfor
ss = 0

if (debug_mode gt 0) then begin
	print, ' '
	print, '"data" array size is ', strtrim(num_elements,2)
	stop, 'Check out structure definitions in data...'
endif

;
;	Once structures are defined, then read the netCDF variables into "data"
;	8.	NCDF_VARGET: Read the data from the variables.
;
if used_var_cnt gt 0 then for i=0, used_var_cnt-1 do begin
	k = used_var_ind[i]

	; get nest counts
	nc = (*var_inq[k].nesting_ptr).nest_cnt
	value = 0
	nest_level = var_inq[k].nest_level
	var_type = var_inq[k].type
	err_result = 0
	catch, error_status
	if error_status ne 0 then begin
		status = BAD_FILE_DATA
		print, "WARNING: couldn't read variable ", var_inq[k].name, $
			" from file:"
		print, "    ", !error_state.msg
		catch, /cancel
		continue
	endif
	NCDF_VARGET, fid, k, value
	catch, /cancel
	val_ndims = size(value, /n_dim)
;	if !version.release ge 6.0 and (((val_ndims ge 1) and $
;		(num_elements eq 1)) or ((val_ndims gt 1) and $
;		(num_elements gt 1))) then begin
		; idl 6+ won't cause an error on NCDF_VARGET even when
		; the file is truncated.  Check and warn if the value
		; is all zeros
;		dum = where(value ne 0, cnt)
;		if cnt eq 0 then $
;			print, get_routine_name() + '  field ', var_inq[k].name, $
;				' contains all zeros - might be missing/invalid data'
;	endif
	if ( var_type eq 7 ) then begin
		value = string( temporary(value) )
	endif else if (var_type eq 6) or $
		(var_type eq 9) then begin
		val_dims = size(value, /dimensions)
		if val_dims[0] eq 0 then begin
			if var_type eq 6 then $
				value = complex(temporary(value)) $
			else $
				value = dcomplex(temporary(value))
		endif else begin
			; build a complex number array from a float/double array
			step = val_dims[0]
			n_comp = n_elements(value) / step
			comp_ind = indgen(n_comp) * step
			if (n_elements(val_dims) eq 1) then $
				val_dims = [1] $
			else $
				val_dims = val_dims[1:n_elements(val_dims)-1]
			if var_type eq 6 then begin
				tmpval = complexarr(val_dims)
				if step eq 1 then $
					tmpval[comp_ind] = complex(value[comp_ind]) $
				else $
					tmpval[indgen(n_comp)] = complex(value[comp_ind], value[comp_ind+1])
			endif else begin
				tmpval = dcomplexarr(val_dims)
				if step eq 1 then $
					tmpval[comp_ind] = dcomplex(value[comp_ind]) $
				else $
					tmpval[indgen(n_comp)] = dcomplex(value[comp_ind], value[comp_ind+1])
			endelse
			value = tmpval
			tmpval = 0
			comp_ind = 0
		endelse
		endif
	if ((var_inq[k].unlim_dim ge 0) or (common_dim ge 0)) $
		and (ind_ndim eq 1) then begin
		; extract subarray for the given indices
		if var_inq[k].dim[var_inq[k].ndims-1] ne dim_unlimited then begin
			print, 'WARNING: Not subindexing variable ', var_inq[k].name, $
				' on dimension ', dim_name[dim_unlimited], $
				': not the final dimension for the variable.'
		endif else case size(value, /n_dimensions) of
			1: value = temporary(value[indices])
			2: value = temporary(value[*,indices])
			3: value = temporary(value[*,*,indices])
			4: value = temporary(value[*,*,*,indices])
			5: value = temporary(value[*,*,*,*,indices])
			6: value = temporary(value[*,*,*,*,*,indices])
			7: value = temporary(value[*,*,*,*,*,*,indices])
		endcase
	endif
	if (var_inq[k].make_ptr eq 1) then begin
		if (var_inq[k].dont_contract ne 0) and (num_elements eq 1) then $
			value = ptr_new(temporary(value)) $
		else begin
			; the last dimension becomes the pointer array
			fill = 0
			if (var_inq[k].natts gt 0) then begin
				atts = *var_inq[k].atts
				ind = (where(atts.name eq 'fillValue'))[0]
				if (ind ge 0) then $
					fill = atts[ind].val
				atts = 0
			endif
			new_val = 0
			ptr_dims = var_inq[k].struct_dims
			if (ptr_dims[0] eq 0) and (num_elements eq 1) then $
				n_ptr_dims = 0 $
			else begin
				if have_struct eq 0 then begin
					if (ptr_dims[0] eq 1) and (num_elements eq 1) then $
						n_ptr_dims = 0 $
					else $
						n_ptr_dims = 1
				endif else begin
					n_ptr_dims = var_inq[k].n_struct_dims
					if (n_ptr_dims ne 1) or (ptr_dims[0] ne 1) then begin
						if (n_ptr_dims ge 0) then $
							ptr_dims = ptr_dims[0:n_ptr_dims-1]
						new_val = make_array(dim=ptr_dims, type=10)
					endif
				endelse
			endelse
			if (size(fill, /type) eq 7) then begin
				contract_arrays, value, new_val, fill, n_ptr_dims=n_ptr_dims
			endif else begin
				contract_arrays, value, new_val, n_ptr_dims=n_ptr_dims
			endelse
			value = new_val
			new_val = 0
		endelse
	endif
	if (num_elements eq 1) and (var_inq[k].make_single_ptr gt 0) then $
		value = ptr_new(temporary(value))

	error_status = 0
	catch, error_status
	if error_status ne 0 then begin
		type_conversion_error:
		print, "WARNING: Couldn't populate field ", var_inq[k].name, ":"
		print, "    ", !error_state.msg
		catch, /cancel
		on_ioerror, null
		if keyword_set(halt_on_error) then $
			stop, 'Halting - check VALUE, DATA.', var_inq[k].name, ' before continuing'
		continue
	endif
	on_ioerror, type_conversion_error
	if (ind_ndim gt 0 and ind_nelem lt n_elements(data)) then begin
		; indices parameter provided, but there are fewer indices than elements in data
		n_ind = n_ind-1
		case nest_level of
			0:  data[0:n_ind].(nc[0]) = value
			1:  data[0:n_ind].(nc[0]).(nc[1]) = value
			2:  data[0:n_ind].(nc[0]).(nc[1]).(nc[2]) = value
			3:  data[0:n_ind].(nc[0]).(nc[1]).(nc[2]).(nc[3]) = value
			4:  data[0:n_ind].(nc[0]).(nc[1]).(nc[2]).(nc[3]).(nc[4]) = value
			else: begin
				; more than 4 nest levels; build up the command and execute
				cmd = 'data[0:n_ind]' + strjoin(string(lindgen(nest_level+1), $
					format='(".(nc[", I0, "])")')) + ' = value'
				dum = execute(cmd)
			end
		endcase
	endif else begin
		case nest_level of
			0:  data.(nc[0]) = value
			1:  data.(nc[0]).(nc[1]) = value
			2:  data.(nc[0]).(nc[1]).(nc[2]) = value
			3:  data.(nc[0]).(nc[1]).(nc[2]).(nc[3]) = value
			4:  data.(nc[0]).(nc[1]).(nc[2]).(nc[3]).(nc[4]) = value
			else: begin
				; more than 4 nest levels; build a string command and execute
				cmd = 'data' + strjoin(string(lindgen(nest_level+1), $
					format='(".(nc[", I0, "])")')) + ' = value'
				dum = execute(cmd)
			end
		endcase
	endelse
	catch, /cancel
	on_ioerror, null
	value = 0
endfor

;
;	Close the netCDF file
;	9.	NCDF_CLOSE: Close the file.
;
NCDF_CLOSE, fid

cleanup:
;
;	Free up Pointers before exiting
;
for k=0, finq.nvars-1 do begin
	if PTR_VALID( var_inq[k].ptr ) then PTR_FREE, var_inq[k].ptr
	if PTR_VALID( var_inq[k].nesting_ptr ) then begin
		nesting_ptr = var_inq[k].nesting_ptr
		for jj=0,n_elements(*nesting_ptr)-1 do begin
			if PTR_VALID( (*nesting_ptr)[jj].str_ptr ) then $
				PTR_FREE, (*nesting_ptr)[jj].str_ptr
		endfor
		PTR_FREE, nesting_ptr
	endif
	if PTR_VALID( var_inq[k].atts ) then PTR_FREE, var_inq[k].atts
endfor

;	remove temporary file, if temp_dir used
if use_temp_file then begin
	if (debug_mode gt 0) then print, 'Removing temporary file: ', full_file 
	file_delete, full_file
endif

return
end