Index: /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/budget.sh
===================================================================
--- /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/budget.sh	(revision 4045)
+++ /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/budget.sh	(revision 4045)
@@ -0,0 +1,11 @@
+#--2D emissions - year 2000
+file=/prodigfs/project/input4MIPs/CMIP6/CMIP/PNNL-JGCRI/CEDS-2017-05-18/atmos/mon/CO2-em-anthro/gn/v20170519/CO2-em-anthro_input4MIPs_emissions_CMIP_CEDS-2017-05-18_gn_200001-201412.nc
+
+#--compute annual mean emissions in Tg/yr
+cdo outputf,%12.10g,10 -divc,1.e9 -mulc,86400. -mulc,365 -fldsum -mul -timmean -seltimestep,1/12 $file -gridarea $file
+
+#--3D emissions - year 2000
+file=/prodigfs/project/input4MIPs/CMIP6/CMIP/PNNL-JGCRI/CEDS-2017-08-30/atmos/mon/CO2-em-AIR-anthro/gn/v20170907/CO2-em-AIR-anthro_input4MIPs_emissions_CMIP_CEDS-2017-08-30_gn_200001-201412.nc
+
+#--compute annual mean emissions in Tg/yr
+cdo outputf,%12.10g,10 -divc,1.e9 -mulc,86400. -mulc,365 -fldsum -mul -timmean -vertsum -seltimestep,1/12 $file -gridarea $file
Index: /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_144x143
===================================================================
--- /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_144x143	(revision 4045)
+++ /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_144x143	(revision 4045)
@@ -0,0 +1,12 @@
+gridtype  = lonlat
+gridsize  = 20592
+xname     = lon
+xunits    = degrees
+yname     = lat
+yunits    = degrees
+xsize     = 144
+ysize     = 143
+xfirst    = -180
+xinc      = 2.5
+yfirst    = 90
+yinc      = -1.267605633802817
Index: /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_256x257
===================================================================
--- /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_256x257	(revision 4045)
+++ /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_256x257	(revision 4045)
@@ -0,0 +1,14 @@
+gridtype  = lonlat
+gridsize  = 65792
+xname     = lon
+xlongname = Longitude
+xunits    = degrees_east
+yname     = lat
+ylongname = Latitude
+yunits    = degrees_north
+xsize     = 256
+ysize     = 257
+xfirst    = -180
+xinc      = 1.40625
+yfirst    = 90
+yinc      = -0.703125
Index: /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_48x36
===================================================================
--- /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_48x36	(revision 4045)
+++ /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_48x36	(revision 4045)
@@ -0,0 +1,12 @@
+gridtype  = lonlat
+gridsize  = 1776
+xname     = lon
+xunits    = degrees
+yname     = lat
+yunits    = degrees
+xsize     = 48
+ysize     = 37
+xfirst    = -180
+xinc      = 7.5
+yfirst    = 90
+yinc      = -5
Index: /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_96x95
===================================================================
--- /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_96x95	(revision 4045)
+++ /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/grid-lmdz-lonlat_96x95	(revision 4045)
@@ -0,0 +1,12 @@
+gridtype  = lonlat
+gridsize  = 9216
+xname     = lon
+xunits    = degrees
+yname     = lat
+yunits    = degrees
+xsize     = 96
+ysize     = 96
+xfirst    = -180
+xinc      = 3.75
+yfirst    = 90
+yinc      = -1.89474
Index: /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/netcdf.pro
===================================================================
--- /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/netcdf.pro	(revision 4045)
+++ /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/netcdf.pro	(revision 4045)
@@ -0,0 +1,72 @@
+pro netcdfshow,filename
+;
+fileID = ncdf_open(filename)
+
+;
+;find the number of variables
+fileinq_struct = ncdf_inquire(fileID)
+nvars = fileinq_struct.nvars
+;
+print,filename
+;
+print,filename
+print, nvars, ' variables in file'
+;
+for varndx = 0,nvars-1 do begin
+
+   ; get the name, datatype, dims, number of attributes
+     varstruct = ncdf_varinq(fileID,varndx)
+   ; print the variable index, name, datatype, dims
+     print, '--------------------------------------------------------'
+     print,varndx,'    ',varstruct.name, '     ',varstruct.datatype
+     print,'dims = ',varstruct.dim
+   ; loop through attributes
+     for attndx = 0, varstruct.natts-1 do begin
+
+       ; get attribute name, then use it to get the value
+         attname = ncdf_attname(fileID,varndx,attndx)
+         ncdf_attget,fileID,varndx,attname,value
+
+       ; print name and value to file
+         print,attname, '     ',string(value)
+     endfor
+endfor ; attribute loop
+end
+;
+pro netcdfread_era, filename, variable_name, data_variable, scale_factor, add_offset, dims
+; This procedure will read netCDF data and place it in an IDL variable 
+; INPUT: filename - a string variable that includes the filepath
+;        variable_name - a string that must match exactly that produced by        
+;                        ncdfshow.pro
+; OUTPUT: data_variable - a user supplied variable for the data
+;         dims - a vector of the dimensions
+
+; get fileID, variable ID
+  fileID = ncdf_open(filename)
+  varID = ncdf_varid(fileID,variable_name)
+
+; get the data and dimensions
+  ncdf_varget,fileID, varID, data_variable
+  ncdf_attget,fileID, varID, 'scale_factor', scale_factor
+  ncdf_attget,fileID, varID, 'add_offset',   add_offset
+  dims = size(data_variable,/dimensions)
+end
+;---------------------
+;---LMDZ
+pro netcdfread, filename,variable_name, data_variable, dims
+; This procedure will read netCDF data and place it in an IDL variable 
+; INPUT: filename - a string variable that includes the filepath
+;        variable_name - a string that must match exactly that produced by        
+;                        ncdfshow.pro
+; OUTPUT: data_variable - a user supplied variable for the data
+;         dims - a vector of the dimensions
+
+; get fileID, variable ID
+  fileID = ncdf_open(filename)
+  varID = ncdf_varid(fileID,variable_name)
+
+; get the data and dimensions
+  ncdf_varget,fileID, varID, data_variable
+  dims = size(data_variable,/dimensions)
+end
+;---------------------
Index: /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/netcdfwrite.pro
===================================================================
--- /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/netcdfwrite.pro	(revision 4045)
+++ /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/netcdfwrite.pro	(revision 4045)
@@ -0,0 +1,258 @@
+; 1.11 01/10/07
+;
+; Takes an IDL structure and writes to NetCDF much as netcdf_read
+; but from the opposite direction. 
+;
+; Dimension variables are identified as the first vector in the structure 
+; of a certain length, and linked to all following vectors and multi-D 
+; arrays with dimensions the same size. 
+;
+; The approach is simple and powerful but with one big flaw - arrays
+; with non-unique dimension sizes (e.g. fltarr[36,36,5] rather than
+; fltarr[36,37,5]) cannot be written unless:
+;
+;   a) A quick and dirty solution that breaks the NetCDF dimension
+;      links: the DIM_ALL keyword is used.
+;
+;   b) A tedious but exact solution: the dimensions are explicilty
+;      specified with the DIMENSIONS keyword.
+;
+; It will be possible to actually define which are dimensions etc... 
+; ...in a later version... perhaps by defining "d_latitude" or "d_pressure"
+; as indicating a dimension? Or if this is not flexible enough
+; just write your own netcdf file following the Chapter 6 of the IDL 
+; Scientific Data Formats manual
+;
+; Parameters
+; ----------
+;
+; I:    file     string containing filename and path to be written
+; I:    field    IDL structure to be written to file, arranged as 
+;                in the example below
+;
+; Keywords
+; --------
+;
+; CLOBBER - overwrite any existing file (otherwise, we fail with error
+;       status if there's an exisitng file there)
+; 
+; DIM_ALL - treat all vectors as dimensions. This is needed in order to
+;       write arrays non-unique dimension sizes, e.g. fltarr[36,36,5]
+;
+; DIMENSIONS - specify the exact dimensions and links. See below for
+;      an example.
+;
+; ATTRIBUTES - pass in a structure containing the attributes for each 
+;      dimension and variable, in the order they are found in the data 
+;      structure we are writing to file. The structure should contain a
+;      strarr for each attribute. The attribute name when written to
+;      file will be the variable name in the struture, in lowercase.
+;
+; Example
+; -------
+;
+; To add "units" and "long_name" attributes to each variable 
+; for a partial compatability with the NetCDF standard, e.g. when 
+; writing the 3D array/structure ozone:
+;
+;      ozone = {lon:fltarr(96),lat:fltarr(73),pressure:fltarr(50),$
+;         ozone:fltarr(96,73,50)}
+;
+;      attributes = {units:strarr(4),long_name:strarr(4)}
+;      attributes.units = ['degrees_east','degrees_north','Pa','kg kg-1']
+;      attributes.long_name = ['Longitude','Latitude','Pressure',$
+;         'Ozone']
+;
+;      netcdfwrite, 'filename', ozone, attributes=attributes
+;
+; On the unix command line "ncdump -h filename" produces:
+;
+;      netcdf filename {
+;      dimensions:
+;         LON = 96 ;
+;         LAT = 73 ;
+;         PRESSURE = 50 ;
+;      variables:
+;         float LON(LON) ;
+;                 LON:units = "degrees_east" ;
+;                 LON:long_name = "Longitude" ;
+;         float LAT(LAT) ;
+;                 LAT:units = "degrees_north" ;
+;                 LAT:long_name = "Latitude" ;
+;         float PRESSURE(PRESSURE) ;
+;                 PRESSURE:units = "Pa" ;
+;                 PRESSURE:long_name = "Pressure" ;
+;         float OZONE(PRESSURE, LAT, LON) ;
+;                 OZONE:units = "kg kg-1" ;
+;                 OZONE:long_name = "Ozone" ;
+;      }
+;
+; Note that NCDUMP reverses the order of the dimensions in multi-
+; dimensional arrays compared to IDL.
+;
+; To explicitly specify the dimensions and linkages (not needed in most
+; cases including this example: it's only needed if the rule of matching 
+; dimensions by length breaks down because there are two dimensions of the
+; same length):
+;
+;      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]]
+;
+;      netcdfwrite, 'filename', ozone, attributes=attributes, $
+;        dimensions=dimensions
+;
+;
+;
+; AJG 22/5/2003
+;
+
+; ------------------------------------------------------------------
+;
+; Add a non-string variable to the NetCDF file
+;
+; INPUTS
+;
+; var_type            - IDL type code - see IDL help for SIZE function
+; ncid, name, dim_ids - parameters passed to ncdf_vardef (see IDL help here too)
+;
+; RETURNS
+;
+; NetCDF variable ID
+;
+function add_var, var_type, ncid, name, dim_ids
+
+case var_type of
+  ; byte
+  1: var_id = ncdf_vardef(ncid, name, dim_ids, /BYTE)
+
+  ; int
+  2: var_id = ncdf_vardef(ncid, name, dim_ids, /SHORT)
+
+  ; long
+  3: var_id = ncdf_vardef(ncid, name, dim_ids, /LONG)
+
+  ; float
+  4: var_id = ncdf_vardef(ncid, name, dim_ids, /FLOAT)
+
+  ; double
+  5: var_id = ncdf_vardef(ncid, name, dim_ids, /DOUBLE)
+
+  ; catch all (will probably catch us out some-day...)
+  else: var_id = ncdf_vardef(ncid, name, dim_ids, /FLOAT) 
+endcase
+
+return, var_id
+
+end
+
+; ------------------------------------------------------------------
+;
+; Main function
+;
+
+pro  netcdfwrite, file, field, clobber=clobber, dim_all=dim_all, $
+        attributes=attributes, dimensions=dimensions
+
+   on_error, 2 ; on error, return to caller (could be improved
+               ; at a later date to catch IO errors and close the
+	       ; file correctly).
+
+   if ~keyword_set(clobber) then clobber = 0
+
+   tag_names = tag_names(field)
+   nTags = N_TAGS(field)
+
+   if keyword_set(attributes) then begin
+     ; Attributes are being supplied. Check they match the
+     ; data structure also supplied.
+     att_names = tag_names(attributes)
+     nAtts = n_tags(attributes)
+     nInAtt = lonarr(nAtts)
+     for i=0,nAtts-1 do nInAtt[i] = n_elements(attributes.(i))
+     if ~array_equal(nInAtt,replicate(nTags,nAtts)) then message, $
+       'Each attribute must have a member for each data structure element'
+   endif
+  
+   ncid = ncdf_create(file, clobber = clobber)   
+
+   ; Run through all data "tags" to identify "dimensions" and
+   ; "variables". For this to work, the data structure must always 
+   ; have dimensions before normal variables  
+   struct_dims, field, dimension_sizes, is_string_max, var_type, $
+     nDimensions, dim_num, iDimLink, dim_all=dim_all, dimensions=dimensions
+
+   dim_ids = intarr(nTags)
+   var_ids = intarr(nTags)
+   string_id = 0 ; a unique ID for any hidden string dimensions
+
+   ; initialise file, defining the dimensions and variables
+   ncdf_control, ncid, /fill
+   for iTag=0, nTags-1 do begin
+   
+      if dim_num[iTag] eq 0 then begin
+      
+         ; Each dimension will have an associated variable
+	 ; containing the values at each point
+         dim_ids[iTag] = ncdf_dimdef(ncid, tag_names[iTag], $
+	    dimension_sizes[iTag])
+
+         var_ids[iTag] = add_var(var_type[iTag], ncid, tag_names[iTag], $
+           dim_ids[iTag])
+	    
+      endif else begin      
+         if is_string_max[iTag] eq 0 then begin
+
+           var_ids[iTag] = add_var(var_type[iTag], ncid, tag_names[iTag], $
+             dim_ids[reform(iDimlink[iTag,0:dim_num[iTag]-1])])
+	    
+	 endif else begin
+	 
+            ; Deal with the string array special case
+
+            ; * Note that a string variable will gain an extra, hidden
+            ; dimension when we write it to file because of NetCDF 
+            ; limitations.
+	    
+	    ; Create a new dimension for the array of characters which forms
+	    ; a string, with its length being that of the longest string
+	    ; in the original array of strings. Don't bother defining
+	    ; a variable for the dimension.
+	    string_id = string_id+1
+	    string_dim = ncdf_dimdef(ncid, $
+	       'STRING'+string(format='(i0)', string_id), $
+	       is_string_max[iTag])
+
+            ; Combine the "IDL" dimensions with the extra hidden
+	    ; string dimension (which must be the first, "fastest varying")
+	    string_dim_ids = $
+	       [string_dim, dim_ids[reform(iDimlink[iTag,0:dim_num[iTag]-1])]]            
+            var_ids[iTag] = ncdf_vardef(ncid, tag_names[iTag], $
+               string_dim_ids, /CHAR)	 
+	             
+	 endelse
+      endelse
+      
+      ; Write any associated atributes
+      if keyword_set(attributes) then begin
+        for iAtt = 0, nAtts-1 do begin
+          ncdf_attput, ncid, var_ids[iTag], strlowcase(att_names[iAtt]), $
+            attributes.(iAtt)[iTag] 
+        endfor
+      endif
+
+   endfor
+   ncdf_control, ncid, /endef
+
+   ; write the data 
+   for iTag=0, nTags-1 do begin
+   
+      ncdf_varput, ncid, var_ids[iTag], field.(iTag)
+      
+   endfor
+
+   ; close the file
+   ncdf_close, ncid
+   
+end
Index: /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/regrid.pro
===================================================================
--- /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/regrid.pro	(revision 4045)
+++ /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/regrid.pro	(revision 4045)
@@ -0,0 +1,141 @@
+pro regrid
+filename='/data/obolmd/CMIP6/CO2_EMISSIONS_48x36/flux_2D_CO2_1867.nc'
+print, filename
+NETCDFREAD,filename,'flux',flux,dimflux
+NETCDFREAD,filename,'lat',lat,dimlat0
+NETCDFREAD,filename,'lon',lon,dimlon0
+NETCDFREAD,filename,'TIME',time,dimtime0
+dimlat=dimlat0(0)
+dimlon=dimlon0(0)
+dimtime=dimtime0(0)
+print, 'dim flux=', dimflux
+A = float([	[3./4., 1./8., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1./8.],$
+		[1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8.],$
+		[1./8., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1./8., 3./4.]   ])
+A = float([	[1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]   ])
+fluxinit=flux
+fluxinit=flux
+flux_check=flux
+for lo=0,dimlon-1 do begin
+for la=0,dimlat-1 do begin
+flux_check(lo,la,*) = invert(A) ## transpose(fluxinit(lo,la,*))
+endfor
+endfor			
+m_bloq = make_array(dimlon,dimlat,12,value=0)		; Matrice boolÃ©enne "mois Ã  bloquer ou non"
+if total(where(flux_check lt 0)) ne -1 then m_bloq(where(flux_check lt 0)) = 1
+; Correction/adaptation de la matrice S&Z en fonction du masque boolÃ©en m_bloq
+for lo=0,dimlon-1 do begin
+for la=0,dimlat-1 do begin
+whereneg = where(flux_check(lo,la,*) lt 0)		; (12 pts max) Identification de potentiels points Ã  problÃšmes, corrigÃ©s nÃ©gativement
+nbannul=n_elements(whereneg)*(total(whereneg) ne -1)
+flux_corr=flux_check(lo,la,*)				; CrÃ©ation d'un vecteur pour recevoir les valeurs corrigÃ©es, initialisÃ© Ã  flux_check au cas oÃ¹ on n'ait rien Ã  faire d'autre qu'une seule itÃ©ration
+A2 = A							; Je repars de la matrice A initiale, ce pour chaque point de grille
+;							; Potentiellement plusieurs passages pour Ã©liminer toutes les valeurs nÃ©gatives
+while nbannul ne 0 do begin				; Si l'on a effectivement des Ã©missions corrigÃ©es nÃ©gativement...
+m_bloq(lo,la,whereneg) = 1				; Update de la matrice m_bloq
+for m=0,11 do begin
+if m eq 11 then begin					; Pour plus de facilitÃ©, mois prÃ©cÃ©dents et suivants codÃ©s ici
+p=10
+s=0
+endif else if m eq 0 then begin
+p=11
+s=1
+endif else begin
+p = m-1
+s = m+1
+endelse
+if m_bloq(lo,la,m) then begin				; Je traite les mois bloquÃ©s en eux-mÃªmes
+A2(p,m) = 0.
+A2(m,m) = 1.
+A2(s,m) = 0.
+endif							; Fin du cas si l'on est sur un mois bloquÃ©
+if ~m_bloq(lo,la,m) then begin				; Je traite les mois non bloquÃ©s, pour ceux adjacents Ã  un mois bloquÃ©
+if m_bloq(lo,la,p) and m_bloq(lo,la,s) then begin   		; Mois encadrÃ© de deux mois bloquÃ©s
+A2(p,m) = 1./4.
+A2(m,m) = 1./2.
+A2(s,m) = 1./4.
+endif else if m_bloq(lo,la,p) and ~m_bloq(lo,la,s) then begin	; Mois prÃ©cÃ©dent bloquÃ© (uniquement)
+A2(p,m) = 2./8.
+A2(m,m) = 5./8.
+A2(s,m) = 1./8.
+endif else if ~m_bloq(lo,la,p) and m_bloq(lo,la,s) then begin	; Mois suivant bloquÃ© (uniquement)
+A2(p,m) = 1./8.
+A2(m,m) = 5./8.
+A2(s,m) = 2./8.
+endif
+endif								; Fin du cas mois non bloquÃ©
+endfor								; Fin de la boucle sur les mois, balayage de la matrice
+flux_corr = invert(A2) ## transpose(fluxinit(lo,la,*))		; RÃ©-itÃ©ration de la multiplication matricielle avec la matrice A modifiÃ©e (A2)
+whereneg = where(flux_corr lt 0)				; (12 pts max) RÃ©-identification de potentiels nouveaux points Ã  problÃšmes, corrigÃ©s nÃ©gativement
+nbannul=n_elements(whereneg)*(total(whereneg) ne -1)
+endwhile							; Fin du cas oÃ¹ l'on avait des problÃšmes d'Ã©missions corrigÃ©es nÃ©gativement
+; *** IMPORTANT ! *** Pour signaler les mois bloquÃ©s, on prend la convention suivante :
+;	valeur nÃ©gative ou nulle <=> mois bloquÃ©
+;	valeur positive <=> mois Ã  interpolation classique
+flux(lo,la,*) = flux_corr					; En sortie de la boucle while, normalement flux_corr est complÃštement positif
+endfor								; Fin boucle lat
+endfor								; Fin boucle lon
+nbnegtotal = n_elements(where(m_bloq eq 1)) * (total(where(m_bloq eq 1)) ne -1)
+if nbnegtotal ne 0 then flux(where(m_bloq eq 1)) = -flux(where(m_bloq eq 1))		; Je force Ã  des valeurs nÃ©gatives
+;
+month_in_year=12
+nbpoint=1682
+flux2=fltarr(nbpoint,month_in_year)
+flux2(*,*)=0.0
+;
+for l=0,month_in_year-1 do begin
+flux2(0,l)=TOTAL(flux(*,0,l))/float(dimlon)
+flux2(nbpoint-1,l)=TOTAL(flux(*,dimlat-1,l))/float(dimlon)
+endfor
+;
+pt=1
+for j=1,dimlat-2  do begin
+for i=0,dimlon-1  do begin
+for l=0,month_in_year-1 do begin
+  flux2(pt,l)=flux(i,j,l)
+endfor
+pt=pt+1
+endfor
+endfor
+;
+;saving netcdf file
+;
+fluxstruct={vector:fltarr(nbpoint),time:fltarr(month_in_year), $
+            flx_CO2:fltarr(nbpoint,month_in_year) }
+;
+fluxstruct.vector=float(indgen(nbpoint)+1)
+;;fluxstruct.time=float(indgen(month_in_year)+1)
+fluxstruct.time=[15, 45, 75, 105, 135, 165, 195, 225, 255, 285, 315, 345]
+fluxstruct.flx_CO2=flux2
+;
+attributes = {units:strarr(3),long_name:strarr(3)}
+attributes.units = ['vector','days since 1960-01-01','flux']
+attributes.long_name = ['vector', 'time', 'flux']
+;
+dimensions = {isdim:intarr(3), links:intarr(2,3)}
+       dimensions.isdim =  [1,1,0]  ; (1=dimension, 0=variable)
+       dimensions.links = [ [-1,-1],[-1,-1],[0,1]    ]
+;
+netcdfwrite,'/data/obolmd/CMIP6/CO2_EMISSIONS_48x36/flux_1D_CO2_1867.nc',fluxstruct,clobber=1, attributes=attributes, dimensions=dimensions
+;
+end
Index: /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/regrid_v1.job
===================================================================
--- /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/regrid_v1.job	(revision 4045)
+++ /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/regrid_v1.job	(revision 4045)
@@ -0,0 +1,281 @@
+#PBS -N process_co2
+#PBS -S /bin/bash
+#PBS -q week # there exist: short, day, days3, week...
+#PBS -k eo  # to write the output of stdin
+### Max memory
+#PBS -l vmem=10gb   # virtual memory
+#PBS -l  mem=10gb
+
+#--updated on 4/5/2017 with improved Sheng & Zwiers algorithm, ThL
+#--corrected some interpolation preprocessing (compared to v4)
+#--updated on 9/5/2017 with output separation BB/anthro for SO2, NOx and NH3
+#--corrected on 22/6/2017 for BB: undef values zeroed before remapping
+#
+#--input directory for anthropogenic (non-BB) emissions
+dirinPNNL_2D="/prodigfs/project/input4MIPs/CMIP6/CMIP/PNNL-JGCRI/CEDS-2017-05-18/"
+dirinPNNL_3D="/prodigfs/project/input4MIPs/CMIP6/CMIP/PNNL-JGCRI/CEDS-2017-08-30/"
+
+#--LMDz grid information
+#grid="144x143"
+#nbpoint=$((144*141+2))
+grid="48x36"
+nbpoint=$((48*35+2))
+
+gridfile="/home/oboucher/CMIP6/GRID/grid-lmdz-lonlat_"${grid}
+
+#--output directory
+dirout="/data/obolmd/CMIP6/CO2_EMISSIONS_"${grid}"/"
+if [ ! -d ${dirout} ] ; then mkdir -p ${dirout} ; fi
+
+#--year
+for year in {1850..2014}
+#for year in {2000..2000}
+do
+
+#--species 
+for species in "CO2"
+do
+
+#--finding correct file for PNNL data
+if [ $year -ge 1750 -a $year -lt 1800 ]; then
+year1=1750
+year2=1799
+elif [ $year -ge 1800 -a $year -lt 1850 ]; then 
+year1=1800
+year2=1849
+elif [ $year == 1850 ]; then
+year1=1850
+year2=1850
+elif [ $year -ge 1851 -a $year -lt 1900 ]; then 
+year1=1851
+year2=1899
+elif [ $year -ge 1900 -a $year -lt 1950 ]; then 
+year1=1900
+year2=1949
+elif [ $year -ge 1950 -a $year -lt 2000 ]; then 
+year1=1950
+year2=1999
+elif [ $year -ge 2000 -a $year -lt 2015 ]; then 
+year1=2000
+year2=2014
+else
+echo 'Houston we have a problem for the PNNL data'
+exit 1
+fi
+
+#--input file PNNL surface emissions
+filename2D=${dirinPNNL_2D}/atmos/mon/${species}-em-anthro/gn/v20170519/${species}-em-anthro_input4MIPs_emissions_CMIP_CEDS-2017-05-18_gn_${year1}01-${year2}12.nc
+#--input file PNNL air emissions
+filename3D=${dirinPNNL_3D}/atmos/mon/${species}-em-AIR-anthro/gn/v20170907/${species}-em-AIR-anthro_input4MIPs_emissions_CMIP_CEDS-2017-08-30_gn_${year1}01-${year2}12.nc
+
+#--output files
+filenameout1=${dirout}flux_${species}_${year}.nc
+filenameout2=${dirout}flux_2D_${species}_${year}.nc
+filenameout3=${dirout}flux_1D_${species}_${year}.nc
+
+echo ${filename2D} ${filename3D} ${filenameout1} ${filenameout2} ${filenameout3}
+rm -f ${filenameout1} ${filenameout2} ${filenameout3}
+
+#--unfortunately idl not happy with PNNL netcdf files so need to ferretize files
+#--I also sum over sectors for the 2D file and over altitude for the 3D file and I extract the correct year as well
+rm -f rewrite.jnl
+cat << eod > rewrite.jnl
+use "${filename2D}"
+use "${filename3D}"
+set region/t=16-jan-${year}:16-dec-${year}
+let flux=${species}_EM_ANTHRO[d=1,k=@sum]+${species}_EM_AIR_ANTHRO[d=2,k=@sum]
+save/clobber/file="${filenameout1}" flux
+eod
+
+#--run ferret script
+ferret << eod 
+go rewrite.jnl
+exit
+eod
+rm -f rewrite.jnl ferret.jnl
+
+#--remap to LMDz grid
+echo cdo remapcon,${gridfile} -chname,FLUX,flux  ${filenameout1} ${filenameout2}
+cdo remapcon,${gridfile} -chname,FLUX,flux ${filenameout1} ${filenameout2}
+
+#--Improved Sheng & Zwiers algorithm + transform into vector
+rm -f regrid.pro
+cat << eod >> regrid.pro
+pro regrid
+filename='${filenameout2}'
+print, filename
+NETCDFREAD,filename,'flux',flux,dimflux
+NETCDFREAD,filename,'lat',lat,dimlat0
+NETCDFREAD,filename,'lon',lon,dimlon0
+NETCDFREAD,filename,'TIME',time,dimtime0
+dimlat=dimlat0(0)
+dimlon=dimlon0(0)
+dimtime=dimtime0(0)
+print, 'dim flux=', dimflux
+A = float([	[3./4., 1./8., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1./8.],$
+		[1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 0., 1./8., 3./4., 1./8.],$
+		[1./8., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1./8., 3./4.]   ])
+A = float([	[1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],$
+		[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]   ])
+fluxinit=flux
+fluxinit=flux
+flux_check=flux
+for lo=0,dimlon-1 do begin
+for la=0,dimlat-1 do begin
+flux_check(lo,la,*) = invert(A) ## transpose(fluxinit(lo,la,*))
+endfor
+endfor			
+m_bloq = make_array(dimlon,dimlat,12,value=0)		; Matrice boolÃ©enne "mois Ã  bloquer ou non"
+if total(where(flux_check lt 0)) ne -1 then m_bloq(where(flux_check lt 0)) = 1
+; Correction/adaptation de la matrice S&Z en fonction du masque boolÃ©en m_bloq
+for lo=0,dimlon-1 do begin
+for la=0,dimlat-1 do begin
+whereneg = where(flux_check(lo,la,*) lt 0)		; (12 pts max) Identification de potentiels points Ã  problÃšmes, corrigÃ©s nÃ©gativement
+nbannul=n_elements(whereneg)*(total(whereneg) ne -1)
+flux_corr=flux_check(lo,la,*)				; CrÃ©ation d'un vecteur pour recevoir les valeurs corrigÃ©es, initialisÃ© Ã  flux_check au cas oÃ¹ on n'ait rien Ã  faire d'autre qu'une seule itÃ©ration
+A2 = A							; Je repars de la matrice A initiale, ce pour chaque point de grille
+;							; Potentiellement plusieurs passages pour Ã©liminer toutes les valeurs nÃ©gatives
+while nbannul ne 0 do begin				; Si l'on a effectivement des Ã©missions corrigÃ©es nÃ©gativement...
+m_bloq(lo,la,whereneg) = 1				; Update de la matrice m_bloq
+for m=0,11 do begin
+if m eq 11 then begin					; Pour plus de facilitÃ©, mois prÃ©cÃ©dents et suivants codÃ©s ici
+p=10
+s=0
+endif else if m eq 0 then begin
+p=11
+s=1
+endif else begin
+p = m-1
+s = m+1
+endelse
+if m_bloq(lo,la,m) then begin				; Je traite les mois bloquÃ©s en eux-mÃªmes
+A2(p,m) = 0.
+A2(m,m) = 1.
+A2(s,m) = 0.
+endif							; Fin du cas si l'on est sur un mois bloquÃ©
+if ~m_bloq(lo,la,m) then begin				; Je traite les mois non bloquÃ©s, pour ceux adjacents Ã  un mois bloquÃ©
+if m_bloq(lo,la,p) and m_bloq(lo,la,s) then begin   		; Mois encadrÃ© de deux mois bloquÃ©s
+A2(p,m) = 1./4.
+A2(m,m) = 1./2.
+A2(s,m) = 1./4.
+endif else if m_bloq(lo,la,p) and ~m_bloq(lo,la,s) then begin	; Mois prÃ©cÃ©dent bloquÃ© (uniquement)
+A2(p,m) = 2./8.
+A2(m,m) = 5./8.
+A2(s,m) = 1./8.
+endif else if ~m_bloq(lo,la,p) and m_bloq(lo,la,s) then begin	; Mois suivant bloquÃ© (uniquement)
+A2(p,m) = 1./8.
+A2(m,m) = 5./8.
+A2(s,m) = 2./8.
+endif
+endif								; Fin du cas mois non bloquÃ©
+endfor								; Fin de la boucle sur les mois, balayage de la matrice
+flux_corr = invert(A2) ## transpose(fluxinit(lo,la,*))		; RÃ©-itÃ©ration de la multiplication matricielle avec la matrice A modifiÃ©e (A2)
+whereneg = where(flux_corr lt 0)				; (12 pts max) RÃ©-identification de potentiels nouveaux points Ã  problÃšmes, corrigÃ©s nÃ©gativement
+nbannul=n_elements(whereneg)*(total(whereneg) ne -1)
+endwhile							; Fin du cas oÃ¹ l'on avait des problÃšmes d'Ã©missions corrigÃ©es nÃ©gativement
+; *** IMPORTANT ! *** Pour signaler les mois bloquÃ©s, on prend la convention suivante :
+;	valeur nÃ©gative ou nulle <=> mois bloquÃ©
+;	valeur positive <=> mois Ã  interpolation classique
+flux(lo,la,*) = flux_corr					; En sortie de la boucle while, normalement flux_corr est complÃštement positif
+endfor								; Fin boucle lat
+endfor								; Fin boucle lon
+nbnegtotal = n_elements(where(m_bloq eq 1)) * (total(where(m_bloq eq 1)) ne -1)
+if nbnegtotal ne 0 then flux(where(m_bloq eq 1)) = -flux(where(m_bloq eq 1))		; Je force Ã  des valeurs nÃ©gatives
+;
+month_in_year=12
+nbpoint=${nbpoint}
+flux2=fltarr(nbpoint,month_in_year)
+flux2(*,*)=0.0
+;
+for l=0,month_in_year-1 do begin
+flux2(0,l)=TOTAL(flux(*,0,l))/float(dimlon)
+flux2(nbpoint-1,l)=TOTAL(flux(*,dimlat-1,l))/float(dimlon)
+endfor
+;
+pt=1
+for j=1,dimlat-2  do begin
+for i=0,dimlon-1  do begin
+for l=0,month_in_year-1 do begin
+  flux2(pt,l)=flux(i,j,l)
+endfor
+pt=pt+1
+endfor
+endfor
+;
+;saving netcdf file
+;
+fluxstruct={vector:fltarr(nbpoint),time:fltarr(month_in_year), $
+            flx_${species}:fltarr(nbpoint,month_in_year) }
+;
+fluxstruct.vector=float(indgen(nbpoint)+1)
+;;fluxstruct.time=float(indgen(month_in_year)+1)
+fluxstruct.time=[15, 45, 75, 105, 135, 165, 195, 225, 255, 285, 315, 345]
+fluxstruct.flx_${species}=flux2
+;
+attributes = {units:strarr(3),long_name:strarr(3)}
+attributes.units = ['vector','days since 1960-01-01','flux']
+attributes.long_name = ['vector', 'time', 'flux']
+;
+dimensions = {isdim:intarr(3), links:intarr(2,3)}
+       dimensions.isdim =  [1,1,0]  ; (1=dimension, 0=variable)
+       dimensions.links = [ [-1,-1],[-1,-1],[0,1]    ]
+;
+netcdfwrite,'${filenameout3}',fluxstruct,clobber=1, attributes=attributes, dimensions=dimensions
+;
+end
+eod
+
+#
+#--calling IDL
+#
+/opt/idl-6.4/idl/bin/idl << eod
+.r netcdf
+.r netcdfwrite
+.r struct_dims
+.r regrid
+regrid
+exit
+eod
+#
+
+#--end loop on species
+done
+
+#--deleting output file if already there
+#rm -f ${dirout}co2ff_lmdz_cmip6_${year}.nc
+
+#--rename a few things
+#cdo expr,'flx_co2=FLX_CO2' ${dirout}flux_vector_CO2_${year}.nc ${dirout}co2ff_lmdz_cmip6_${year}.nc
+#ncrename -d VECTOR,vector -v VECTOR,vector ${dirout}co2ff_lmdz_cmip6_${year}.nc
+#ncrename -d TIME,time     -v TIME,time     ${dirout}co2ff_lmdz_cmip6_${year}.nc
+
+#--cleaning up 
+#rm -f ${dirout}flux*_${year}.nc
+rm -f ${filenameout1}
+
+#--end loop on years
+done
+
+#--cleaning the mess
+rm -f ferret*
+rm -f regrid.pro
+rm -f rewrite.jnl
Index: /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/struct_autoconvert.pro
===================================================================
--- /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/struct_autoconvert.pro	(revision 4045)
+++ /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/struct_autoconvert.pro	(revision 4045)
@@ -0,0 +1,59 @@
+;-------------------------------------------------------------
+;+
+; NAME:
+;       STRUCT_AUTOCONVERT
+; PURPOSE:
+;       Convert a structure wih all string items to numeric.
+; CATEGORY:
+; CALLING SEQUENCE:
+;       out = struct_autoconvert(in)
+; INPUTS:
+;       in = Structure with all string items.   in
+; KEYWORD PARAMETERS:
+; OUTPUTS:
+;       out = Converted structure with numeric  out
+;         strings converted to numeric values.
+; COMMON BLOCKS:
+; NOTES:
+;       Notes: Non-numeric items are kept as strings.
+;       Integer items are returned as long ints.
+;       Floating items are returned as double.
+; MODIFICATION HISTORY:
+;       R. Sterner, 2006 Sep 07
+;
+; Copyright (C) 2006, Johns Hopkins University/Applied Physics Laboratory
+; This software may be used, copied, or redistributed as long as it is not
+; sold and this copyright notice is reproduced on each copy made.  This
+; routine is provided as is without any express or implied warranties
+; whatsoever.  Other limitations apply as described in the file disclaimer.txt.
+;-
+;-------------------------------------------------------------
+	function struct_autoconvert, in, help=hlp
+ 
+	if (n_params(0) lt 1) or keyword_set(hlp) then begin
+	  print,' Convert a structure wih all string items to numeric.'
+	  print,' out = struct_autoconvert(in)'
+	  print,'   in = Structure with all string items.   in'
+	  print,'   out = Converted structure with numeric  out'
+	  print,'     strings converted to numeric values.'
+	  print,' Notes: Non-numeric items are kept as strings.'
+	  print,' Integer items are returned as long ints.'
+	  print,' Floating items are returned as double.'
+	  return,''
+	endif
+ 
+	tags = tag_names(in)		; Tags in incoming strutcure.
+	n = n_elements(tags)		; Number of tags.
+	typ = bytarr(n)			; To hold data types.
+	tbl = [7,3,5]			; Convert isnumber output to type.
+ 
+	for i=0,n-1 do typ(i)=tbl(isnumber(in.(i)))	; Find data types.
+ 
+	out = create_struct(tags(0), $			; Create output struct.
+	  to_datatype(in.(0),typ(0)))
+	for i=1,n-1 do out = create_struct(out, $
+	  tags(i),to_datatype(in.(i),typ(i)))
+ 
+	return, out
+ 
+	end
Index: /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/struct_dims.pro
===================================================================
--- /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/struct_dims.pro	(revision 4045)
+++ /TOOLS/CMIP6_FORCING/CO2_EMISSIONS/struct_dims.pro	(revision 4045)
@@ -0,0 +1,206 @@
+; 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
+
