;+
; NAME:
;	write_netCDF.pro
;
; PURPOSE:
;	Write netCDF file given a structure variable
;
; CATEGORY:
;	All levels of processing
;
; CALLING SEQUENCE:  
;	write_netCDF, data, filename, status, path=dir_path, att_file=att_filename, $
;		name_unlimited=name_unlimited, size_unlimited=size_unlimited, $
;		/unlimited, /no_unlimited, /clobber
;
; INPUTS:
;	data = structure variable of input data
;	filename = filename for new netCDF file
;	path = optional directory path for the attributes definition file
;	att_file = optional filename for the attributes definition file
;	name_unlimited = optional name for the unlimited dimension
;	size_unlimited = optional integer, used only if the unlimited
;		dimension is otherwise unspecified.  The unlimited dimension
;		will be set to this size, and any variables with a last
;		dimension of that size will use the unlimited dimension.
;		(If this parameter is unset, the final dimension that shows up 
;		most often in at least 1/3 of the number of actual variables
;		will be the unlimited dimension)
;	unlimited = optional flag - the structure's dimension will be unlimited
;	no_unlimited = optional flag - no unlimited dimension will be defined
;		(overrides /unlimited flag)
;	clobber = optional flag for creating netCDF file;
;		clobber means an existing file will be overwritten
;
;	An external *.att file is used to define attributes (where * = "data" structure name)
;
; OUTPUTS:  
;	status = result status: 0 = OK_STATUS, -1 = BAD_PARAMS, -2 = BAD_FILE,
;			-3 = BAD_FILE_DATA, -4 = FILE_ALREADY_OPENED
;
;	A netCDF file is created and written.
;
; COMMON BLOCKS:
;	None
;
; PROCEDURE:
;	Check for valid input parameters.
;	Use the structure name and optional 'path' variable for the Attributes filename
;		OR use the optional 'att_file' parameter for this filename.
;	If this Attributes definition file exists, then read in the attributes and
;		any other useful information (like dimension labels).
;	Open the netCDF file.
;	Define netCDF dimensions based on the structure and any given dimension labels.
;	Use the structure's tag names for defining the variable names in the netCDF.
;	Transfer any existing attributes into the netCDF file.
;	Once netCDF variables and attributes are defined, then write the structure's data to netCDF file.
;	Close the netCDF file.
;
;	NetCDF IDL Procedures / Process:
;	1. NCDF_CREATE: Call this procedure to begin creating a new file. The new file is put into define mode.
;	2. NCDF_DIMDEF: Create dimensions for the file.
; 	3. NCDF_VARDEF: Define the variables to be used in the file.
;	4. NCDF_ATTPUT: Optionally, use attributes to describe the data.  Global attributes also allowed.
;	4. NCDF_CONTROL, /ENDEF: Leave define mode and enter data mode.
;	5. NCDF_VARPUT: Write the appropriate data to the netCDF file.
;	6. NCDF_CLOSE: Close the file.
;
; MODIFICATION HISTORY:
;	9/20/99		Tom Woods		Original release code, Version 1.00
;
;-

pro	write_netCDF, data, filename, status, path=path, att_file=att_file, $
		name_unlimited=name_unlimited, size_unlimited=size_unlimited, $
		unlimited=unlimited, no_unlimited=no_unlimited, clobber=clobber

;
;	Supported IDL types (and their codes) for netCDF 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 write_netCDF will convert a LONG64 (14) type to netCDF DOUBLE,
;	with some resulting loss of precision for very large values 
;	(more than 15 digits.)
;	Unsigned integers are stored as twos-complement signed integers,
;	and write_netCDF will convert UINT (12) to INT or ULONG (13) to LONG
;	with no loss of data upon re-read; ULONG64 (15) will be converted
;	to netCDF DOUBLE, with similar loss of precision as for LONG64.
;

;
;	Generic "status" values
;
OK_STATUS = 0
BAD_PARAMS = -1
BAD_FILE = -2
BAD_FILE_DATA = -3
FILE_ALREADY_OPENED = -4
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: write_netCDF, data, filename, status, path=dir_path, ' + $
		'att_file=att_filename, name_unlimited=name_unlimited, ' + $
		'size_unlimited=size_unlimited, ' + $
		'/unlimited, /no_unlimited, /clobber'
	return
endif
dsize = size(data)
if (dsize[0] ne 1) or (dsize[2] ne 8) then begin
	print, 'ERROR: write_netCDF requires the data to be a 1-dimensional structure array'
	return
endif

num_elements = dsize[1]

if keyword_set(unlimited) then $
	unlimited = 1 $
else $
	unlimited = 0

if keyword_set(no_unlimited) then begin
	no_unlimited = 1
	unlimited = 0
endif else $
	no_unlimited = 0

struct_unlimited = 0
if (num_elements gt 1) or unlimited then struct_unlimited = 1
if no_unlimited then struct_unlimited = 0

if (n_params(0) lt 2) then begin
	filename = ''
	read, 'Enter filename for the new netCDF file : ', filename
	if (strlen(filename) lt 1) then return
endif

dir_path = !CIPS_NETCDF_ATT_PATH
att_filename = ''
if keyword_set(path) then dir_path = path
IF keyword_set(att_file) THEN BEGIN
  att_filename = att_file
ENDIF ELSE BEGIN
  att_filename = dir_path + "/" + att_filename
ENDELSE
null_ptr = ptr_new()

name_val = { name: '', val: '' }
var_att = { var: '', att: null_ptr }
att_def = 0
var_adim = { var: '', dim: null_ptr }
dim_def = 0

;	variables for tracking unlimited dimension info
uses_unlimited = 0
unlim_dim_name = ''
if struct_unlimited then $
	unlim_dim_size = num_elements $
else $
	unlim_dim_size = 0
all_have_unlimited = 1
none_have_unlimited = 1

;
;	Use the structure name and optional 'path' variable for the Attributes filename
;		OR use the optional 'att_file' parameter for this filename
;	If this Attributes definition file exists, then transfer those attributes into the netCDF file
;		OR else don't write any attributes to the netCDF file (except for fillValues)
;	Variable names are case insensitive, but variable names in the netCDF file
;	will follow the letter casing used in the attribute file.
;
on_ioerror, bad_att_file
openr,alun, att_filename, /get_lun
cur_varid = -1		; GLOBAL default start
astr = ''

vname = ''
cur_var_att = 0
while not( eof(alun) ) do begin
	readf, alun, astr
	astr = strtrim(strcompress(astr),2)
	; variable and attribute names must start with a letter or underscore
	pos_let = (stregex(astr, '[A-Za-z_]'))[0]
	alen = strlen(astr)
	if ( alen gt 0 ) and (pos_let eq 0) then begin
		; the first non-whitespace is a valid name start character
		if (strpos( astr, '=' ) lt 0) then begin
			pos_colon = strpos( astr, ':' )
			if ( pos_colon gt 0 ) then begin
				if (vname ne '') then begin
					cur_avar = var_att
					if (size(cur_var_att, /n_dimensions) gt 0) then begin
						cur_avar.var = vname
						cur_avar.att = ptr_new(cur_var_att)
						if (size(att_def, /n_dimensions) eq 0) then $
							att_def = [ cur_avar ] $
						else $
							att_def = [ temporary(att_def), cur_avar ]
					endif
				endif
				cur_var_att = 0
				;
				;	got a variable name (could be GLOBAL)
				;
				vname = strtrim(strmid(astr,0,pos_colon),2)
				; look for dimension labels after the colon
				dims = strtrim(strmid(astr,pos_colon+1))
				if (dims ne '') then begin
					; save the dimension labels, seen as '(dim1,dim2)'
					dims = strmid(dims,1,strlen(dims)-2)
					; the first dim label for GLOBAL is for
					; the unlimited dim
					if (strupcase(vname) eq 'GLOBAL') then begin
						unlim_dim_name = (strsplit(dims, ',', /extract))[0]
					endif else begin
						cur_dvar = var_adim
						cur_dvar.var = vname
						dim_labels = strsplit(dims, ',', /extract)
						if unlim_dim_name ne '' then begin
							if dim_labels[n_elements(dim_labels)-1] $
								eq unlim_dim_name then $
								none_have_unlimited = 0 $
							else $
								all_have_unlimited = 0
						endif
						cur_dvar.dim = ptr_new(dim_labels)
						if (size(dim_def, /n_dimensions) eq 0) then $
							dim_def = [ cur_dvar ] $
						else $
							dim_def = [ temporary(dim_def), cur_dvar ]
					endelse
				endif else if unlim_dim_name ne '' then $
					all_have_unlimited = 0
			endif
		endif else begin
			;
			;	attribute definition (name = text)
			;
			adef = str_sep( astr, "=" )
			n_adef = n_elements(adef)
			aname = strtrim( adef[0], 2 )
			atext = adef[1]
			jj = 2   ;  merge any other text together
			while (jj lt n_adef) do begin
				atext = atext + '=' + adef[jj]
				jj = jj + 1
			endwhile
			atext = strtrim( atext, 2 )
			;
			;	check for special substitutions:  $FILE, $DATE, $TIME
			;	also check for "$" at end of atext which means it needs to read another line
			;
			pos_dollar = strpos( atext, "$", /reverse_search )
			if ( pos_dollar ge 0 ) then begin
				alen = strlen(atext)
				while ( pos_dollar eq (alen-1) ) and not( eof(alun) ) do begin
					readf, alun, astr
					atext = strtrim(strmid(atext, 0, alen-1),2)
					atext = atext + ' ' + strtrim( strcompress( astr ), 2 )
					alen = strlen(atext)
					pos_dollar = strpos( atext, "$", /reverse_search )
				endwhile
				up_atext = strupcase( atext )
				repeat begin
					pos1 = strpos( up_atext, "$FILE" )
					if (pos1 ge 0) then begin
						newtext = strmid( atext, 0, pos1 ) + filename
						if ( pos1 lt (alen-5) ) then newtext = newtext + strmid( atext, pos1+5, alen-pos1-5 )
						atext = newtext
						alen = strlen(atext)
						up_atext = strupcase( atext )
					endif
				endrep until pos1 lt 0
				repeat begin
					pos1 = strpos( up_atext, "$DATE" )
					if (pos1 ge 0) then begin
						newtext = strmid( atext, 0, pos1 ) + systime()
						if ( pos1 lt (alen-5) ) then newtext = newtext + strmid( atext, pos1+5, alen-pos1-5 )
						atext = newtext
						alen = strlen(atext)
						up_atext = strupcase( atext )
					endif
				endrep until pos1 lt 0
				repeat begin
					pos1 = strpos( up_atext, "$TIME" )
					if (pos1 ge 0) then begin
						newtext = strmid( atext, 0, pos1 ) + systime()
						if ( pos1 lt (alen-5) ) then newtext = newtext + strmid( atext, pos1+5, alen-pos1-5 )
						atext = newtext
						alen = strlen(atext)
						up_atext = strupcase( atext )
					endif
				endrep until pos1 lt 0
			endif
			up_atext = ''

			tmp_att = name_val
			tmp_att.name = aname
			tmp_att.val = atext
			if (size(cur_var_att, /n_dimensions) eq 0) then $
				cur_var_att = [ tmp_att ] $
			else $
				cur_var_att = [ temporary(cur_var_att), tmp_att ]
		endelse
	endif
endwhile
; process the last variable's attribute info
if (vname ne '') then begin
	cur_avar = var_att
	if (size(cur_var_att, /n_dimensions) gt 0) then begin
		cur_avar.var = vname
		cur_avar.att = ptr_new(cur_var_att)
		if (size(att_def, /n_dimensions) eq 0) then $
			att_def = [ cur_avar ] $
		else $
			att_def = [ temporary(att_def), cur_avar ]
	endif
endif

close, alun
free_lun, alun
goto, end_att_file

bad_att_file:
	print, 'WARNING: write_netCDF could not find attributes file = ', att_filename

end_att_file:
	on_ioerror, NULL

;	If the attribute file provided a label for the unlimited dim,
;	assume it gets used - the whole structure is unlimited if
;	exactly all or none of the variables use the unlimited dim label
if (struct_unlimited eq 0) and (not no_unlimited) and $
	(unlim_dim_name ne '') then begin
	uses_unlimited = 1
	;if (none_have_unlimited gt 0) or (all_have_unlimited gt 0) then $
	;	struct_unlimited = 1
endif

if (size(name_unlimited, /type) eq 7) then $
	unlim_dim_name = name_unlimited

;
;	Do survey of variables and nested structures
;	to verify limitation on dimensions of arrays and nested structures
;	
;	LIMITATIONS:  4 dimensions on arrays and 6 nested structures
;
;	Use internal name structure for tracking any nested structures
;
temp_def = { name : ' ', isVar : 0B, tag_index : 0L, var_size : lonarr(10), $
	nest_level : 0, struct_index : null_ptr, dim_index : lonarr(10), $
	var_ptr: null_ptr, string_maxlen : 0L, dim_unlimited : 0, is_null : 0, $
	last_dim_size : 0, is_single : 0, total_dims : 0, tag_ind_str: '' }
var_def = temp_def

;
;	define first structure entry into "var_def" for the "data" structure
;
var_def[0].name = tag_names( data, /structure_name )
var_def[0].isVar = 0
var_def[0].tag_index = 0
var_def[0].var_size = size( data )
var_def[0].nest_level = 0
var_def[0].struct_index = ptr_new( [ 0L ] )
temp_dim = lonarr(10) - 1
var_def[0].dim_index = temp_dim
var_def[0].var_ptr = ptr_new( data[0] )
if num_elements gt 1 then $
	var_def[0].last_dim_size = num_elements
var_def[0].is_single = (num_elements eq 1)
var_def[0].total_dims = 1

next_var = 1
level_index = ptr_new( [ 1L ] )
extra_var = n_tags( data )
nest_level = 0

has_string = 0
has_complex = 0

while (extra_var gt 0) do begin
	;
	; each level of nested structures are appended to var_def
	;
	var_def = [ var_def, replicate( temp_def, extra_var ) ]
	if (nest_level gt 0) then j_start = (*level_index)[nest_level-1] else j_start = 0
	j_end = (*level_index)[nest_level] - 1
	extra_var = 0	
	for j=j_start, j_end do begin
		;
		; only process structure definitions
		;
		if ( var_def[j].isVar eq 0 ) then begin
			temp_index = *var_def[j].struct_index
			ti = var_def[temp_index].tag_index
			theData = *var_def[j].var_ptr
			theWholeData = 0
			wholedatasize = 0
			nn = var_def[j].nest_level
			if ( nn gt 0 ) then begin
				parentName = var_def[ temp_index[nn] ].name + '.'
			endif else begin
				parentName = ''
			endelse
			tnames = tag_names( theData )
			k_total = n_tags( theData ) - 1
			for k=0, k_total do begin
				theName = parentName + tnames[k]
				theVar = theData[0].(k)
				nelem = n_elements(theVar)
				single_data = var_def[j].is_single
				tempsize = size( theVar )
				; process fillValue attribute, if any
				att_var_ind = -1
				att_ind = -1
				atts = 0
				fill = 0
				string_max = 0
				var_type = size(theVar, /type)
				if (size(att_def, /n_dimensions) gt 0) then begin
					; get the fillValue from attributes, if any
					att_var_ind = (where(strupcase(att_def.var) $
						eq strupcase(theName)))[0]
					if (att_var_ind ge 0) then begin
						theName = att_def[att_var_ind].var
						atts = *att_def[att_var_ind].att
						att_ind = (where(atts.name eq 'fillValue'))[0]
						if (att_ind ge 0) then $
							fill = atts[att_ind].val
					endif
				endif
				is_null = 0
				var_single_data = single_data and (nelem eq 1)
				while (var_type eq 10) and (tempsize[0] eq 0) and $
					(not is_null) do begin
					; drill down through nested single pointers
					if theVar[0] ne null_ptr then begin
						theVar = *(temporary(theVar[0]))
						tempsize = size( theVar )
						var_type = size(theVar, /type)
						var_single_data = n_elements(theVar) eq 1
					endif else is_null = 1
				endwhile
				if (not single_data) and ((var_type eq 7) or (var_type eq 10)) $
					and ((size(theWholeData))[0] eq 0) then begin
					; for pointer and string types, need the whole set of data
					; for the field to find the best fit dimensions
					case var_def[j].nest_level of
						0 : theWholeData = data
						1 : theWholeData = data.(ti[1])
						2 : theWholeData = data.(ti[1]).(ti[2])
						3 : theWholeData = data.(ti[1]).(ti[2]).(ti[3])
						4 : theWholeData = data.(ti[1]).(ti[2]).(ti[3]).(ti[4])
						else : begin
							; need to build a command to execute
							cmd = 'theWholeData = data' + strjoin(string( $
								lindgen(var_def[j].nest_level) + 1L, $
								format='(".(ti[", I0, "])")'))
							dum = execute(cmd)
						end
					endcase
					wholedatasize = size( theWholeData )
				endif
				if (not single_data) and ((var_type eq 7) or (var_type eq 10)) then $
					theVar = theWholeData.(k)
				if (var_type eq 10) then begin
					; pointer type
					vind = where(theVar ne null_ptr)  ; indices of "valid" pointers
					if (vind[0] ge 0) then begin
						var_type = size(*theVar[vind[0]], /type)
						; should only be pointer to array of non-structs
						if (var_type eq 8) then begin
							print, "ERROR: can't have pointer to a structure"
							goto, cleanup_no_close
						endif
					endif else is_null = 1
					array_bounds, theVar, tmpDim, fill, string_max
					if (tmpDim[0] eq 0) then begin
						print, "ERROR: array types don't match for ", $
							"pointer array field ", theName
						goto, cleanup_no_close
					endif
					if (not single_data) and (not is_null) then begin
						; remove dimensions that belong to theWholeData
						n_data_dim = wholedatasize[0]
						n_tmp_dim = n_elements(tmpDim) - n_data_dim
						tmpDim = (temporary(tmpDim))[0:n_tmp_dim-1]
					endif
					if (size(fill, /type) eq 7) then begin
						; set the fillValue attribute
						if (att_ind ge 0) then begin
							atts[att_ind].val = fill
							*att_def[att_var_ind].att = atts
						endif else begin
							; make a new attribute for fillValue
							new_attr = name_val
							new_attr.name = 'fillValue'
							new_attr.val = fill
							if (size(atts, /n_dimensions) gt 0) then begin
								atts = [ atts, new_attr ]
								*att_def[att_var_ind[0]].att = atts
							endif else begin
								new_var_att = var_att
								new_var_att.var = theName
								new_var_att.att = ptr_new([ new_attr ])
								if (size(att_def, /n_dimensions) gt 0) then $
									att_def = [ att_def, new_var_att ] $
								else $
									att_def = [ new_var_att ]
							endelse
						endelse
					endif
					; don't really need theVar anymore,
					; so save some time by not processing it
					if is_null then begin
						tmpDim = [1]
						var_type = 1
						n_tmpdims = 1
					endif else begin
						n_tmpdims = n_elements(tmpDim)
						if (n_tmpdims gt 1) and $
							(tmpDim[n_tmpdims-1] eq 1) then begin
							last_ind = max(where(tmpDim ne 1))
							if last_ind lt 0 then $
								tmpDim = [1] $
							else $
								tmpDim = tmpDim[0:last_ind]
							n_tmpdims = n_elements(tmpDim)
						endif
					endelse
					if tmpDim[0] eq 0 then $
						tempsize = [ 0, var_type ] $
					else $
						tempsize = [ n_tmpdims, tmpDim, var_type ]
					var_single_data = single_data and (max(tmpDim) le 1)
				endif else begin
					fill = 0
					if (size(theVar, /type) eq 7) then $
						string_max = max(strlen(theVar))
				endelse
				; remove fillValue if not needed
				if (size(fill, /type) ne 7) then begin
					if (att_ind[0] ge 0) then begin
						keep_ind = where(atts.name ne 'fillValue')
						if (keep_ind[0] ge 0) then $
							*att_def[att_var_ind[0]].att = atts[keep_ind] $
						else begin
							ptr_free, att_def[att_var_ind[0]].att
							att_def[att_var_ind[0]].att = null_ptr
							att_def = att_def[where(strupcase(att_def.var) ne theName)]
						endelse
					endif
				endif
				is_complex = 0
				if var_type eq 7 then $
					has_string = 1 $
				else if (var_type eq 6) or (var_type eq 9) then $
					is_complex = 1
				has_complex = has_complex or is_complex

				var_def[next_var].name = theName
				var_def[next_var].isVar = 1
				var_def[next_var].tag_index = k
				var_def[next_var].nest_level = nest_level
				var_def[next_var].tag_ind_str = $
					strjoin(string([ ti, k ], format='(I5)'), ',')
				var_def[next_var].dim_index = temp_dim
				if (not single_data) then $
					var_def[next_var].last_dim_size = var_def[j].last_dim_size $
				else if (not var_single_data) then $
					var_def[next_var].last_dim_size = tempsize(tempsize[0])
				var_def[next_var].is_single = var_single_data
				if var_def[j].is_single then $
					var_def[next_var].total_dims = 1 > tempsize[0] $
				else if (var_type eq 8) and (n_elements(theVar) eq 1) then $
					var_def[next_var].total_dims = $
						var_def[j].total_dims $
				else $
					var_def[next_var].total_dims = $
						var_def[j].total_dims + tempsize[0]
				if (var_def[next_var].total_dims gt 7) then begin
					; complex #s need an extra dim for real/imag parts
					print, 'ERROR:  write_netCDF has a limitation of 7 total dimensions for its variables'
					goto, cleanup_no_close
				endif
				var_def[next_var].var_size = tempsize
				var_def[next_var].string_maxlen = string_max
				var_def[next_var].is_null = is_null
				;
				;	if structure, then need to set it up special
				;
				if (var_type eq 8) then begin
					var_def[next_var].isVar = 0
					var_def[next_var].var_ptr = ptr_new( theVar[0] )
					var_def[next_var].nest_level = nest_level + 1
					var_def[next_var].struct_index = ptr_new([ temp_index, next_var ])
					extra_var = extra_var + n_tags( theVar[0] )
				endif else $
					var_def[next_var].struct_index = ptr_new( temp_index )
				theVar = 0
				next_var = next_var + 1
			endfor
			theData = 0
			theWholeData = 0
			wholedatasize = 0
		endif
	endfor
	;
	;	get ready for next level of nested structures
	;
	nest_level = nest_level + 1
	*level_index = [ *level_index, next_var ]
endwhile

num_var = next_var		; the maximum number of variables for netCDF file (size of var_def)
if (num_var ne n_elements(var_def)) then begin
	print, 'WARNING: write_netCDF has error in pre-parsing for variable definitions'
endif

if (extra_var gt 0) then begin
	print, 'ERROR:  write_netCDF was not able to process all nested structure variables'
	print, 'ABORTING....'
	; NCDF_CONTROL, fid, /ABORT
	goto, cleanup_no_close
endif

;
;	if unlimited dim not explicitly used, then find a good last dim size
;	(if any) to use for the unlimited dim
;
if (not no_unlimited) and (not struct_unlimited) and $
	(not uses_unlimited) then begin
	if size(size_unlimited, /type) ne 0 then begin
		; if the given size is 1, make the struct unlimited
		if size_unlimited eq 1 then $
			unlim_dim_size = 0 $
		else begin
			; use the given size as the unlimited dim,
			; if any vars use it for their last dim
			dum = where(var_def[var_inds].last_dim_size eq $
				size_unlimited, cnt)
			if cnt ne 0 then $
				unlim_dim_size = long(size_unlimited) $
			else begin
				unlim_dim_size = -1
				print, 'Given size for unlimited dim (', $
					strtrim(size_unlimited,2), $
					') is not used by any variables.', $
					'  No unlimited dim will be used.'
			endelse
		endelse
	endif else begin
		; See if a certain percentage of vars use the same last 
		; dimension size.  If so, make that the unlimited dim
		var_inds = where(var_def.isVar, n_var_inds)
		last_dims = uniq(var_def[var_inds].last_dim_size, $
			sort(var_def[var_inds].last_dim_size))
		last_dims = var_def[var_inds[temporary(last_dims)]].last_dim_size
		n_last_dims = n_elements(last_dims)
		if n_last_dims eq 1 then begin
			; all actual variables have the same last_dim
			all_have_unlimited = 1
			unlim_dim_size = last_dims[0]
		endif else begin
			; look for a dim used in the majority of vars
			unlim_dim_size = 0
			best_cnt = 0
			for i=0,n_last_dims-1 do begin
				dum = where(var_def[var_inds].last_dim_size eq $
					last_dims[i], dim_cnt)
				; from any non-zero last dim sizes in more than 1/3
				; of the variables, choose the largest majority holder
				if (last_dims[i] gt 1) and (dim_cnt ge best_cnt) and $
					(((dim_cnt * 3) gt n_var_inds)) then begin
					unlim_dim_size = last_dims[i]
					best_cnt = dim_cnt
					uses_unlimited = 1
				endif
			endfor
		endelse
	endelse
	if unlim_dim_size eq 0 then begin
		uses_unlimited = 1
		struct_unlimited = 1
	endif else if unlim_dim_size gt 0 then begin
		uind = where(var_def[var_inds].last_dim_size eq unlim_dim_size)
		var_def[var_inds[uind]].dim_unlimited = 1
		uind = 0
	endif
endif

;
;	Open the netCDF file - option to CLOBBER any existing file
;
status = BAD_FILE
if keyword_set(clobber) then fid = NCDF_CREATE( filename, /CLOBBER ) $
else fid = NCDF_CREATE( filename, /NOCLOBBER )
status = OK_STATUS

if (debug_mode gt 0) then print, 'Opened NETCDF file...'

;
;	Define the netCDF dimensions
;	Use the size() function to make dimensions
;	Define the dimension of the structure itself as UNLIMITED (in case want to append to this file)
;

if (debug_mode gt 0) then begin
	print, ' '
	print, 'Number of structures / variables = ', num_var
	print, ' '
	print, 'Defining dimensions and variables...'
	print, '    Index   Dimensions   Data-Type   Name'
	print, '    -----   ----------   ---------   ----'
endif

;
;	record needed dimensions for later definition
;	assume no more than 8 dimensions per variable
;	set val to 0 to indicate the dimension needs to be defined
;
dim_id = replicate( { val: -1, size: 0, label: '', label_num: 0, $
	default_label: 0, string_dim: 0 }, num_var * 8 + 4 )

;	indices into dim_id for common-use dimensions
;	these will be defined if needed
unlimited_dim_ind = 0
if no_unlimited then begin
	unlim_dim_name = ''
	dim_id[unlimited_dim_ind].label = ''
endif else if unlim_dim_name eq '' then $
	dim_id[unlimited_dim_ind].label = 'structure_elements' $
else $
	dim_id[unlimited_dim_ind].label = unlim_dim_name

single_dim_ind = 1
dim_id[single_dim_ind].size = 1
dim_id[single_dim_ind].label = 'single_DIM'

string_dim_ind = 2
dim_id[string_dim_ind].label = 'string'

complex_dim_ind = 3
dim_id[complex_dim_ind].size = 2
dim_id[complex_dim_ind].label = 'complex_number'

forced_struct_unlimited = struct_unlimited and $
	(not (none_have_unlimited or all_have_unlimited))

cur_dim_ind = 4
for k=0,num_var-1 do begin
	var_size = var_def[k].var_size
	if ((k eq 0) or (var_def[k].is_null)) then begin
		if struct_unlimited then begin
			dim_id[unlimited_dim_ind].val = 0
			dim_id[unlimited_dim_ind].size = num_elements
			var_def[k].dim_index[0] = unlimited_dim_ind
		endif
		if not ((k eq 0) and no_unlimited) then $
			continue
	endif
	jnext = 0
	var_type = var_size[var_size[0] + 1]
	; front-load dimension indices for string or complex dimensions
	if (var_type eq 6) or (var_type eq 9) then begin
		var_def[k].dim_index[jnext] = complex_dim_ind
		jnext = jnext + 1
	endif else if (var_type eq 7) then begin
		var_def[k].dim_index[jnext] = string_dim_ind
		jnext = jnext + 1
	endif
	var_dims = 0
	if (var_size[0] gt 0) then begin
		var_dims = var_size[1:var_size[0]]
		for j=1,var_size[0] do begin
			id_ind = -1
			cur_size = var_size[j]
			if (cur_size ge 1) and (var_def[k].isVar or $
				(var_size[0] gt 1) or (var_size[1] gt 1)) then begin
				if cur_size eq 1 then $
					id_ind = single_dim_ind $
				else begin
					; allocate a default dimension definition
					id_ind = cur_dim_ind
					if (k eq 0) then $
						dim_label = dim_id[unlimited_dim_ind].label $
					else $
						dim_label = 'dim' + strtrim(j,2) + $
							'_' + var_def[k].name
					dim_id[cur_dim_ind].label = dim_label
					dim_id[cur_dim_ind].default_label = 1
					dim_id[cur_dim_ind].size = cur_size
					cur_dim_ind = cur_dim_ind + 1
				endelse
				var_def[k].dim_index[jnext] = id_ind
			endif
			jnext = jnext + 1
		endfor
	endif
	;
	;	append dimensions for any structure dimensions
	;
	if (var_def[k].nest_level gt 0) and (var_def[k].isVar) then begin
		for j=var_def[k].nest_level,1,-1 do begin
			ii = 0
			jj = (*var_def[k].struct_index)[j]
			while (var_def[jj].dim_index[ii] ge 0) do begin
				var_def[k].dim_index[jnext] = var_def[jj].dim_index[ii]
				jnext = jnext + 1
				ii = ii + 1
			endwhile
			if ii gt 0 then begin
				if var_dims[0] eq 0 then begin
					var_dims = [ var_def[jj].var_size[1:ii] ]
					;jnext = jnext - 1
				endif else $
					var_dims = [ var_dims, var_def[jj].var_size[1:ii] ]
			endif
		endfor
	endif
	;
	;	append the dimension for "data" structure array for each variable
	;	IF has common dimension
	;
	if var_def[k].isVar and (struct_unlimited or $
		(no_unlimited and (num_elements gt 1))) then begin
		var_def[k].dim_index[jnext] = var_def[0].dim_index[0]
		jnext = jnext + 1
		;if (num_elements gt 1) then $
		var_dims = [ var_dims, num_elements ]
	endif
	n_var_dims = n_elements(var_dims)

	;	if dim_unlimited set, make the last dim index unlimited
	if var_def[k].dim_unlimited then $
		var_def[k].dim_index[jnext-1] = unlimited_dim_ind

	;	match up any dimension labels from the attribute file
	dim_labels = 0
	dim_labels_ind = -1
	if (size(dim_def, /n_dimensions) gt 0) then $
		dim_labels_ind = (where(strupcase(dim_def.var) $
		eq strupcase(var_def[k].name)))[0]
	n_dim_labels = 0
	if (dim_labels_ind ge 0) then begin
		dim_labels = *dim_def[dim_labels_ind].dim
		n_dim_labels = n_elements(dim_labels)
	endif else $
		dim_labels = [ '' ]
	if (n_dim_labels gt 0) then begin
		; shift for complex/string types
		if (var_type eq 6) or (var_type eq 7) or $
			(var_type eq 9) then $
			i = 1 $
		else $
			i = 0
		if none_have_unlimited and uses_unlimited then begin
			dim_labels = [temporary(dim_labels), unlim_dim_name]
			n_dim_labels = n_dim_labels + 1
		endif
		if (var_type eq 7) and ((n_dim_labels gt n_var_dims) or $
			((n_dim_labels eq n_var_dims) and $
			forced_struct_unlimited)) then begin
			i = 0
		endif
		for j=0,n_dim_labels-1 do begin
			dim_label = dim_labels[j]
			if dim_label eq '' then continue
			past_dim_ind = var_def[k].dim_index[j+i+1]
			; shouldn't have any dims after an unlimited dim
			bad_unlim_dim = (dim_label eq unlim_dim_name) and $
				(past_dim_ind ge 0)
			id_ind = where((dim_id.label eq dim_label) and $
				(not dim_id.default_label), label_cnt)
			ind0 = id_ind[0]
			if ind0 lt 0 then begin
				; haven't defined a dim with this label yet
				id_ind = var_def[k].dim_index[j+i]
				if (id_ind lt 0) or (id_ind eq single_dim_ind) then begin
					; if more labels than dims,
					; pad out with dims of size 1
					id_ind = cur_dim_ind
					dim_id[id_ind].label = dim_label
					dim_id[id_ind].size = 1
					var_def[k].dim_index[j+i] = id_ind
					cur_dim_ind = cur_dim_ind + 1
				endif else if (var_type eq 7) and $
					(id_ind eq string_dim_ind) then begin
					id_ind = cur_dim_ind
					dim_id[id_ind].label = dim_label
					dim_id[id_ind].size = var_def[k].string_maxlen
					dim_id[id_ind].string_dim = 1
					var_def[k].string_maxlen = 0
					var_def[k].dim_index[j+i] = id_ind
					cur_dim_ind = cur_dim_ind + 1
				endif else if (dim_id[id_ind].default_label) $
					then begin
					dim_id[id_ind].label = dim_label
					dim_id[id_ind].default_label = 0
				endif
			endif else if (var_type eq 7) and (i eq 0) and (j eq 0) then begin
				; a suggested label for this var's string dim
				str_ind = where(dim_id[id_ind].string_dim)
				str_ind0 = str_ind[0]
				if str_ind0 ge 0 then begin
					; have a string dim by that name,
					; use it (resize if needed)
					id_ind = ind_ind[str_ind0]
					if dim_ind[id_ind].size lt $
						var_def[k].string_maxlen then $
						dim_ind[id_ind].size = var_def[k].string_maxlen
				endif else begin
					; no string dims by that name, create
					id_ind = cur_dim_ind
					dim_id[id_ind].label = dim_label
					dim_id[id_ind].size = var_def[k].string_maxlen
					dim_id[id_ind].string_dim = 1
					dim_id[id_ind].label_num = label_cnt + 1
				endelse
				var_def[k].string_maxlen = 0
				var_def[k].dim_index[j+i] = id_ind
			endif else begin
				if (var_type eq 7) and (i eq 0) then $
					j2 = j - 1 $
				else $
					j2 = j
				; use existing label dim_id
				size_ind = where((dim_id[id_ind].size eq var_dims[j2]) $
					and (not dim_id[id_ind].string_dim))
				size_ind0 = size_ind[0]
				unlim_size_match = 0
				if size_ind0 ge 0 then $
					unlim_size_match = (id_ind[size_ind0] $
						eq unlimited_dim_ind)
				if unlim_size_match and bad_unlim_dim then begin
					if n_elements(size_ind) gt 1 then $
						size_ind0 = size_ind[1] $
					else $
						size_ind0 = -1
				endif
				if (dim_id[ind0].size eq 0) and $
					(not bad_unlim_dim) then begin
					; this sets the unlimited dim size
					dim_id[ind0].size = var_dims[j2]
					var_def[k].dim_index[j+i] = ind0
				endif else if (size_ind0 lt 0) then begin
					; nothing of the right size,
					; or shouldn't be unlimited
					; relabel and add a label_num
					id_ind = var_def[k].dim_index[j+i]
					dim_id[id_ind].label = dim_label
					dim_id[id_ind].label_num = label_cnt + 1
					dim_id[id_ind].default_label = 0
				endif else begin
					var_def[k].dim_index[j+i] = id_ind[size_ind0]
				endelse
			endelse
		endfor
	endif
	
	;	now flag the dimensions to be defined
	j = 0
	while (var_def[k].dim_index[j] ge 0) do begin
		dim_id[var_def[k].dim_index[j]].val = 0
		j = j + 1
	endwhile
	
	if size(var_dims, /n_dimensions) gt 0 then $
		var_def[k].var_size = [n_elements(var_dims), var_dims, var_type] $
	else $
		var_def[k].var_size = [0, var_type]
endfor

; find some reasonable string length dimensions
string_max = 0
string_ind = where(var_def.string_maxlen gt 0, cnt)
if cnt ne 0 then begin
	n_str_ind = n_elements(string_ind)
	string_max = max(var_def[string_ind].string_maxlen)
	dim_id[string_dim_ind].size = string_max	; fix string length to max length used
endif

; now defined the needed dimensions
used_dim_ind = where(dim_id.val ge 0, used_dim_cnt)
if used_dim_cnt gt 0 then for i=0, used_dim_cnt-1 do begin
	cur_ind = used_dim_ind[i]
	cur_dim_id = dim_id[cur_ind]
	dim_label = cur_dim_id.label
	if (cur_dim_id.label_num gt 1) then $
		dim_label = dim_label + '_' + strtrim(cur_dim_id.label_num, 2)
	if (cur_ind eq unlimited_dim_ind) and (not no_unlimited) then $
		dim_id[cur_ind].val = NCDF_DIMDEF( fid, dim_label, /UNLIMITED ) $
	else $
		dim_id[cur_ind].val = NCDF_DIMDEF( fid, dim_label, cur_dim_id.size )
	cur_dim_id = 0
endfor

;
;	find the sorting order for var_def based on tag_index order
;
sorted_var_inds = sort(var_def.tag_ind_str)

;
;	Now define the netCDF variables
;	Use the structure's tag names for defining the variable names in the netCDF
;
first_var=0
for vdi=0,num_var-1 do begin
	;
	;  only process real variables (not structure definitions)
	;
	k = sorted_var_inds[vdi]
	if (var_def[k].isVar) then begin
		var_size = var_def[k].var_size
		data_type = var_size[ var_size[0] + 1 ]
		if (debug_mode gt 0) then print, k, var_size[0], data_type, '   ', var_def[k].name
		;
		;	now make dimension array
		;
		ii = 0
		the_dim = 0
		while (var_def[k].dim_index[ii] ge 0) do begin
			if (ii eq 0) then the_dim = [ dim_id[ var_def[k].dim_index[ii] ].val ] $
			else the_dim = [ the_dim, dim_id[ var_def[k].dim_index[ii] ].val ]
			ii = ii + 1
		endwhile
		;
		;	now make variable in a big case statement for different data type
		;
		if size(the_dim, /n_dimensions) eq 0 then $
			case data_type of
			1:	var_defid = NCDF_VARDEF( fid, var_def[k].name, /BYTE )
			2:	var_defid = NCDF_VARDEF( fid, var_def[k].name, /SHORT )
			3:	var_defid = NCDF_VARDEF( fid, var_def[k].name, /LONG )
			4:	var_defid = NCDF_VARDEF( fid, var_def[k].name, /FLOAT )
			5:	var_defid = NCDF_VARDEF( fid, var_def[k].name, /DOUBLE )
			6: 	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /FLOAT )
			7: 	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /CHAR )
			9: 	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /DOUBLE )
			;
			; NOTE: change this block of code when netCDF finally supports
			; 	unsigned integer types and 64-bit integers
			; meanwhile, store unsigned int (12) as SHORT,
			; 	unsigned long (13) as LONG,
			;	signed/unsigned long64 (14,15) as DOUBLE
			;
			12:	var_defid = NCDF_VARDEF( fid, var_def[k].name, /SHORT )
			13:	var_defid = NCDF_VARDEF( fid, var_def[k].name, /LONG )
			14: begin
				print, "WARNING: write_netCDF can't handle 64-bit integers, casting to DOUBLE"
				var_defid = NCDF_VARDEF( fid, var_def[k].name, /DOUBLE )
				end
			15: begin
				print, "WARNING: write_netCDF can't handle 64-bit integers, casting to DOUBLE"
				var_defid = NCDF_VARDEF( fid, var_def[k].name, /DOUBLE )
				end
			else: begin
				print, 'WARNING: write_netCDF error in variable type, assuming float'
				var_defid = NCDF_VARDEF( fid, var_def[k].name)	; assume it is /FLOAT ???
				end
		endcase else case data_type of 
			1:	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /BYTE )
			2:	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /SHORT )
			3:	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /LONG )
			4:	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /FLOAT )
			5:	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /DOUBLE )
			6: 	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /FLOAT )
			7: 	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /CHAR )
			9: 	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /DOUBLE )
			;
			; NOTE: change this block of code when netCDF finally supports
			; 	unsigned integer types and 64-bit integers
			; meanwhile, store unsigned int (12) as SHORT,
			; 	unsigned long (13) as LONG,
			;	signed/unsigned long64 (14,15) as DOUBLE
			;
			12:	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /SHORT )
			13:	var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /LONG )
			14: begin
				print, "WARNING: write_netCDF can't handle 64-bit integers, casting to double"
				var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /DOUBLE )
				end
			15: begin
				print, "WARNING: write_netCDF can't handle 64-bit integers, casting to double"
				var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim, /DOUBLE )
				end
			else: begin
				print, 'WARNING: write_netCDF error in variable type, assuming float'
				var_defid = NCDF_VARDEF( fid, var_def[k].name, the_dim )	; assume it is /FLOAT ???
				end
		endcase
		if (first_var eq 0) then begin
			var_id = replicate( var_defid, num_var )
			first_var = 1
		endif
		var_id[k] = var_defid
		if (var_defid lt 0) then begin
			status = BAD_FILE_DATA
			print, 'ERROR: only the rightmost dimension can be UNLIMITED', $
				' (field ', var_def[k].name, ')'
			print, 'Aborting...'
			goto, cleanup
		endif
	endif
endfor

; if (debug_mode gt 0) then stop, 'Check out the "var_id"...'

;
;	Now write the attributes that match with actual variables
;	to the file, if any attributes were read
;
if (size(att_def, /n_dimensions) gt 0) then begin
	for i=0,n_elements(att_def)-1 do begin
		if att_def[i].att eq null_ptr then continue
		vname = att_def[i].var
		atts = *att_def[i].att
		;
		;	got a variable name (could be GLOBAL)
		;		find out which variable index if it exists
		;
		if (strupcase(vname) eq 'GLOBAL') then begin
			cur_varid = -1
		endif else begin
			cur_varid = -2
			jj = 0
			test_vname = strupcase(vname)
			while ( jj lt num_var ) and (cur_varid lt 0) do begin
				if (test_vname eq strupcase(var_def[jj].name) ) then cur_varid = jj
				jj = jj + 1
			endwhile
			if (cur_varid lt -1) then begin
				print, 'WARNING: write_netCDF variable NOT found for attribute ', test_vname
			endif
		endelse
		if (cur_varid ge -1) then begin
			;
			;	now define attributes (either GLOBAL or as part of variable)
			;
			for j=0,n_elements(atts)-1 do begin
				aname = atts[j].name
				atext = atts[j].val
				if (atext eq '') then atext = "''"
				if ( cur_varid lt 0 ) then begin
					NCDF_ATTPUT, fid, /GLOBAL, aname, atext
				endif else begin
					NCDF_ATTPUT, fid, var_id[cur_varid], aname, atext
				endelse
			endfor
		endif
	endfor
endif

;
;	Once netCDF variables and attributes are defined, then write the structure's data to netCDF file
;
NCDF_CONTROL, fid, /ENDEF
for vdi=0,num_var-1 do begin
	;
	;  only process real variables (not structure definitions)
	;
	k = sorted_var_inds[vdi]
	if (var_def[k].is_null) then begin
		if struct_unlimited then $
			theData = bytarr(dim_id[unlimited_dim_ind].size) + 1 $
		else $
			theData = [ 1B ]
		NCDF_VARPUT, fid, var_id[k], theData
	endif else if var_def[k].isVar then begin
		ti = var_def[*var_def[k].struct_index].tag_index
		k_ti = var_def[k].tag_index
		case var_def[k].nest_level of
			0 : theData = data.(k_ti)
			1 : theData = data.(ti[1]).(k_ti)
			2 : theData = data.(ti[1]).(ti[2]).(k_ti)
			3 : theData = data.(ti[1]).(ti[2]).(ti[3]).(k_ti)
			4 : theData = data.(ti[1]).(ti[2]).(ti[3]).(ti[4]).(k_ti)
			else : begin
				; need to build a command to execute
				cmd = 'theData = data' + strjoin(string( $
					lindgen(var_def[k].nest_level) + 1L, $
					format='(".(ti[", I0, "])")')) + '.(k_ti)'
				dum = execute(cmd)
			end
		endcase
		dsize = size(theData, /n_dim)
		data_type = size(theData, /type)
		; reduce single-element pointer
		while (dsize eq 0) and (data_type eq 10) do begin
			theData = *theData
			data_type = size(theData, /type)
			dsize = size(theData, /n_dim)
		endwhile
		if (data_type eq 10) then begin
			; theData has pointers
			vind = where(theData ne null_ptr)
			if (vind[0] lt 0) then $
				tmpData = bytarr(1) $
			else $
				tmpData = *theData[vind[0]]
			data_type = size(tmpData, /type)
			var_size = var_def[k].var_size
			tmpSize = var_size[0]
			if (tmpSize gt 0) then begin
				tmpSize = var_size[indgen(tmpSize)+1]
				; make elem the fillValue
				elem = 0
				if (size(att_def, /n_dimensions) gt 0) then begin
					att_ind = where(strupcase(att_def.var) eq $
						strupcase(var_def[k].name), count)
					if (count gt 0) then begin
						if att_def[att_ind[0]].att ne null_ptr then begin
							atts = *att_def[att_ind[0]].att
							att_ind = where(atts.name eq 'fillValue', count)
							if (count gt 0) then $
								elem = atts[att_ind[0]].val
						endif
					endif
				endif
				expand_arrays, theData, tmpData, tmpSize, elem
			endif else begin
				if (dsize gt 1) then $
					for j=1,dsize-1 do $
						tmpData = [tmpData, [*theData[j]]]
			endelse
			theData = tmpData
			tmpData = 0
		endif
		; convert complex type to 2-element float type
		if (data_type eq 6) or (data_type eq 9) then begin
			tmpSize = size(theData, /dimensions)
			if tmpSize[0] eq 0 then $
				tmpSize = [2] $
			else $
				tmpSize = [2, temporary(tmpSize)]
			evens = indgen(n_elements(theData)) * 2
			case data_type of
				6: begin
					tmpData = fltarr(tmpSize)
					tmpData[evens] = float(theData)
				end
				9: begin
					tmpData = dblarr(tmpSize)
					tmpData[evens] = double(theData)
				end
			endcase
			tmpData[evens + 1] = imaginary(theData)
			theData = tmpData
			tmpData = 0
		endif
		NCDF_VARPUT, fid, var_id[k], theData
	endif
endfor

cleanup:
;
;	Close the netCDF file
;
NCDF_CLOSE, fid

cleanup_no_close:
;
;	clean up pointer heap before leaving
;

ind = where(ptr_valid(var_def.struct_index), cnt)
if cnt gt 0 then ptr_free, var_def[ind].struct_index
ind = where(ptr_valid(var_def.var_ptr), cnt)
if cnt gt 0 then ptr_free, var_def[ind].var_ptr

if (size(att_def, /n_dimensions) gt 0) then begin
	ind = where(ptr_valid(att_def.att), cnt)
	if cnt gt 0 then ptr_free, att_def[ind].att
endif

if (size(dim_def, /n_dimensions) gt 0) then begin
	ind = where(ptr_valid(dim_def.dim), cnt)
	if cnt gt 0 then ptr_free, dim_def[ind].dim
endif

if ptr_valid( level_index ) then ptr_free, level_index

return
end