source: TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/struct_dims.pro @ 4007

; 1.1 01/10/07

; Work out the "dimensions" to a structure for writing
; to netcdf file, based on some simplistic rules which 
; fail if more than one dimension has the same length
;
; Where by these rules it is not obvious which should be
; dimensions and which variables relate to which dimension
; there are two options - either a quick and dirty fix, or 
; tediously indicating and linking dimensions by hand:
;
;   DIM_ALL - just flags all vectors as dimensions in the
;       netcdf file. The NetCDF linkage between dimensions
;       and variables will be broken and should be ignored. 
;       This is fine if you're planning on no more than 
;       reading the data into an IDL structure again, but 
;       will produce confusing results with ncdump
;
;   DIMENSIONS - A structure indicating dimensions and
;       the dimensions applying to variables, up to
;       a maximum of five. E.g. when there are four
;       elements in the IDL structure (the netcdf_write.pro
;       example):
; 
;       dimensions = {isdim:intarr(4), links:intarr(5,4)}
;       dimensions.isdim =  [1,1,1,0]  ; (1=dimension, 0=variable)
;       dimensions.links = [[-1,-1,-1,-1,-1],[-1,-1,-1,-1,-1],$
;          [-1,-1,-1,-1,-1],[0,1,2,-1,-1]]
;
; Run through all data "tags" to identify "dimensions" and
; "variables". For this to work, the data structure must always 
; have dimensions before normal variables  
;
; AJG 31/7/06 extracted from netcdf_write.pro
;
pro struct_dims, field, dimension_sizes, is_string_max, var_type, $
  nDimensions, dim_num, iDimLink, dim_all=dim_all, dimensions=dimensions

; Output codes from the IDL SIZE routine - numerical types
; that we may convert to float in our simplistic universe
; though there may be some truncation of information
; (e.g. 64bit int reduced to float). Simplistic universe now
; includes treatment of doubles, though (essential for Julian time)
valid_float = [1,2,3,4,5,12,13,14,15]

; We can make a special (and difficult) exception for string data, 
; too, but only string vectors at the moment. See NetCDF user's 
; guide, "reading and writing character string values". Also
; IDL online help "String data in NetCDF files"
valid_string = [7]

; Find tag names and number
tag_names = tag_names(field)
nTags = N_TAGS(field)

; Output data: information on what we think are dimensions and
; variables
dimension_sizes = lonarr(nTags) ; Made long not int for big files!
is_string_max = intarr(nTags) ; 0=not a string; >0=max string length
var_type = intarr(nTags) ; Variable type based on the output of SIZE
nDimensions = 0
dim_num = intarr(nTags)
iDimlink = intarr(nTags,5) ; up to 5D currently 

; First pass - identify dimensions. All dimensions must
; be vectors. Also identifies variables, and flags string data
for iTag=0, nTags-1 do begin 
   
  tag_size = size(field.(iTag))
  var_type[iTag] = size(field.(iTag),/type) 
           
  if tag_size[0] le 1 then begin

    ; We have found a vector (tag_size[0]=1) or a scalar (=0)
    ; which we will treat as a 1D 1 element vector. The 
    ; following rules apply: If it's a numerical vector, it 
    ; is a dimension only if it's the first of that size. 
    ; Following vectors of the same size will be treated as
    ; variables. (unless keyword DIM_ALL is set)

    ; If scalar pretend it's a 1D 1 element vector
    if tag_size[0] eq 0 then begin
      tag_size[0] = 1
      tag_size[1] = 1
    endif

    iTemp = where( var_type[iTag] eq valid_float, count)
    if count gt 0 then begin
      
      ; We have a numerical value that will be treated as FLOAT or DOUBLE
      
      ; Is this a dimension?
      bDim = 0

      ; (a) If it's the first vector of this length found
      iTemp = where( tag_size[1] eq dimension_sizes, count)
      if count eq 0 then bDim = 1

      ; (b) If the DIM_ALL keyword is set (all vectors treated as dimensions)
      if keyword_set(DIM_ALL) then bDim = 1

      ; (c) If indicated by the DIMENSIONS structure
      if keyword_set(DIMENSIONS) then begin
        if dimensions.isdim[iTag] then bDim = 1 
      endif

      if bDim then begin

        ; We have found a new dimension
        dimension_sizes[iTag] = tag_size[1]
        dim_num[iTag] = 0 ; indicates this IS a dimension
        nDimensions = nDimensions + 1

      endif else begin
            
        ; This is a variable
        dim_num[iTag] = tag_size[0] ; number of dimensions linked           

      endelse

    endif else begin
         
      iTemp = where( tag_size[2] eq valid_string, count)
      if count eq 1 then begin
            
        ; This is a variable, but it's a string array
        ; which needs special treatment
        dim_num[iTag] = tag_size[0] ; number of dimensions linked*
        is_string_max[iTag] = max(strlen(field.(iTag)))  

        ; * Note that a string variable will gain an extra, hidden
        ; dimension when we write it to file because of NetCDF 
        ; limitations.

        ; NB also that a string array cannot be a dimension.
                 
      endif else begin 
        message, "Data type not supported for "+tag_names[iTag]
      endelse
    endelse

  endif else begin

    ; This is a 2D or more variable
    dim_num[iTag] = tag_size[0] ; number of dimensions applicable
         
  endelse 
endfor    
   
; Second pass - link variables to dimensions
for iTag=0, nTags-1 do begin 
   
  tag_size = size(field.(iTag))
      
  ; Is it a variable?      
  if dim_num[iTag] ne 0 then begin

    ; If scalar pretend it's a 1D 1 element vector
    if tag_size[0] eq 0 then begin
      tag_size[0] = 1
      tag_size[1] = 1
    endif
          
    ; find which dimensions apply to this variable by finding
    ; vectors of the same length as the array dimensions
    for iDim = 0, tag_size[0]-1 do begin

      if keyword_set(DIMENSIONS) then begin

        ; Caller has explicitly specified the dimensions and links
        iDimlink[iTag,iDim] = dimensions.links[iDim,iTag]

      endif else begin

        ; Assign dimensions on a first come first served basis. 
        ; by matching sizes. This works fine until you have two
        ; dimensions the same size. Use DIM_ALL or DIMENSIONS in that case
        iDimSizeEqual = where(dimension_sizes eq tag_size[iDim+1], count)
        if count ne 0 then begin
          if count eq 1 then begin
            iDimlink[iTag,iDim] = iDimSizeEqual[0]
          endif else begin
            ; Work out which dimension applies if more than one is of
            ; identical length: do this by order.
            iTagSizeEqual = where(tag_size eq tag_size[iDim+1], tag_ct)
            if tag_ct eq 1 then begin
              iDimlink[iTag,iDim] = iDimSizeEqual[0]
            endif else begin
              tag_order = where(iTagSizeEqual eq iDim+1)
              iDimlink[iTag,iDim] = iDimSizeEqual[tag_order]
            endelse
          endelse
        endif else begin
          message, "Matching dimension not found for "+tag_names[iTag]
        endelse   
     
      endelse

    endfor

  endif        
    
endfor

end