Index: /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_144x143
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_144x143	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_144x143	(revision 4007)
@@ -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/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_256x257
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_256x257	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_256x257	(revision 4007)
@@ -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/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_48x36
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_48x36	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_48x36	(revision 4007)
@@ -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/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_96x95
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_96x95	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid-lmdz-lonlat_96x95	(revision 4007)
@@ -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/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid_vert_lmdz.F90
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid_vert_lmdz.F90	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/grid_vert_lmdz.F90	(revision 4007)
@@ -0,0 +1,69 @@
+ PROGRAM toot
+ IMPLICIT NONE 
+
+!--ICAO profile
+!--flight level given ICAO table
+!--41 levels from ICAO table and 1 dummy level on bottom, 2 on top
+  INTEGER, PARAMETER :: flevel=44, klev=79
+!--feetICAO(flevel)=surface, feetICAO(1)=haute altitude
+  REAL :: feetICAO(flevel)
+!--pressure ICAO in hPa - from top to surface
+  REAL, DIMENSION(flevel), PARAMETER ::               &
+                          pICAO=(/0.,125.,188.,197.,  &
+  206.,217.,227.,238.,250.,262.,274.,287.,301.,315.,  &
+  329.,344.,360.,376.,393.,410.,428.,446.,466.,485.,  &
+  506.,527.,549.,572.,595.,619.,644.,670.,697.,724.,  &
+  753.,782.,812.,843.,875.,908.,942.,977.,1013.,1050./)
+
+  REAL :: feetLMDz, altLMDz(klev)
+  REAL, PARAMETER :: feet2meters=0.3048
+  INTEGER :: k, kk
+
+  REAL,DIMENSION(klev),PARAMETER :: pplay=(/
+&
+    101202.1, 100942.4, 100653.4, 100331.9, 99974.2, 99576.55,
+&
+    99134.67, 98643.84, 98098.98, 97494.48, 96824.3, 96081.91, 95260.24,
+&
+    94351.73, 93348.34, 92241.55, 91022.46, 89681.82, 88210.2, 86598.11,
+&
+    84836.23, 82915.69, 80828.4, 78567.48, 76127.79, 73506.45, 70703.54,
+&
+    67722.81, 64572.39, 61265.38, 57820.47, 54262.17, 50620.74,
+46931.69,  &
+    43234.69, 39571.99, 35986.37, 32518.85, 29206.36, 26079.83,
+23162.71,  &
+    20470.42, 18010.36, 15782.71, 13781.58, 11996.42, 10413.44,
+9016.865,  &
+    7790.02, 6716.147, 5778.999, 4963.233, 4254.644, 3640.276, 3108.45,
+&
+    2648.733, 2251.867, 1909.687, 1615.018, 1361.581, 1143.891,
+957.1722,  &
+    797.2696, 660.5754, 543.9607, 444.7146, 360.4912, 289.2627,
+229.2787,  &
+    179.0303, 137.2198, 102.7344, 74.62312, 52.07758, 34.41441,
+21.05668,  &
+    11.50741, 5.296797, 1.483435 /)
+
+!--initialise ICAO flight levels
+    DO k=3,flevel
+      feetICAO(k)=FLOAT(flevel-k-1)*1000.
+    ENDDO
+    feetICAO(2)= 47000.
+    feetICAO(1)=300000.
+
+    DO k=1,klev
+        DO kk=1,flevel-1
+            IF (pplay(k).LT.pICAO(kk+1)*100. .AND.
+pplay(k).GE.pICAO(kk)*100.) THEN
+               feetLMDz=feetICAO(kk)+(feetICAO(kk+1)-feetICAO(kk)) &
+                            *(pplay(k)-pICAO(kk)*100.)/(pICAO(kk+1)*100.-pICAO(kk)*100.)
+                 altLMDz(k)=feet2meters*feetLMDz*1.e-3 !-en km
+              ENDIF
+       ENDDO
+    ENDDO
+
+    print *,'altLMDz=', altLMDz
+
+    END
+
Index: /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/vct-lmdz
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/vct-lmdz	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/GRID/vct-lmdz	(revision 4007)
@@ -0,0 +1,20 @@
+ 0     0. 0.000
+ 1    64. 0.000
+ 2   196. 0.000
+ 3   457. 0.000
+ 4   939. 0.000
+ 5  1764. 0.000
+ 6  3074. 0.000
+ 7  5018. 0.000
+ 8  7716. 0.000
+ 9 11226. 0.000
+10 15509. 0.000
+11 20088. 0.007 
+12 22615. 0.062
+13 20901. 0.203
+14 16236. 0.396
+15 11063. 0.586
+16  6780. 0.742
+17  3702. 0.856
+18  1590. 0.937
+19     0. 1.000
Index: /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/README
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/README	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/README	(revision 4007)
@@ -0,0 +1,3 @@
+regrid.sh  = current version, needs to be run in a job
+so on ciclad : qsub regrid.job
+
Index: /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/README
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/README	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/README	(revision 4007)
@@ -0,0 +1,7 @@
+ThL 7/2018
+
+Input emission files for scenarios are given by species, with the 2015-2100 period being sampled only at years 2015, 2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090 and 2100.
+Therefore, a time interpolation is necessary.
+
+All time interpolations are performed in one go prior to any file creation (IDL script). The script then operates upon a 2015-2100 continuous input file.
+
Index: /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/interpol_time.pro
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/interpol_time.pro	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/interpol_time.pro	(revision 4007)
@@ -0,0 +1,167 @@
+pro interpol_time
+
+; Les dimensions des fichiers
+nlon=720
+nlat=360
+nt=1032
+
+une_annee=[31,30,29,31,30,31,30,31,31,30,31,30]
+
+jours_in = [15, 45, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349, 1840, 1870, $
+    1899, 1930, 1960, 1991, 2021, 2052, 2083, 2113, 2144, 2174, 5490, 5520, $
+    5549, 5580, 5610, 5641, 5671, 5702, 5733, 5763, 5794, 5824, 9140, 9170, $
+    9199, 9230, 9260, 9291, 9321, 9352, 9383, 9413, 9444, 9474, 12790, 12820, $
+    12849, 12880, 12910, 12941, 12971, 13002, 13033, 13063, 13094, 13124, $
+    16440, 16470, 16499, 16530, 16560, 16591, 16621, 16652, 16683, 16713, $
+    16744, 16774, 20090, 20120, 20149, 20180, 20210, 20241, 20271, 20302, $
+    20333, 20363, 20394, 20424, 23740, 23770, 23799, 23830, 23860, 23891, $
+    23921, 23952, 23983, 24013, 24044, 24074, 27390, 27420, 27449, 27480, $ 
+    27510, 27541, 27571, 27602, 27633, 27663, 27694, 27724, 31040, 31070, $ 
+    31099, 31130, 31160, 31191, 31221, 31252, 31283, 31313, 31344, 31374]
+
+jours_out=make_array(1032)
+jours_out(0) = 15
+for cpt=1,1031 do jours_out(cpt) = jours_out(cpt-1) + une_annee(cpt mod 12)
+
+idin = ncdf_open('/home/tlurton/Scen_EMISS_AER/REGRID/BC-em-bb_nc3.nc')
+varid = ncdf_varid(idin,'BC_em_openburning')
+lonid = ncdf_varid(idin,'lon')
+latid = ncdf_varid(idin,'lat')
+lonbndid = ncdf_varid(idin,'lon_bnds')
+latbndid = ncdf_varid(idin,'lat_bnds')
+timbndid = ncdf_varid(idin,'time_bnds')
+ncdf_varget, idin, varid, emiss_init
+ncdf_attget, idin, varid, "units", emiss_unit
+ncdf_attget, idin, varid, "long_name", emiss_longname
+ncdf_varget, idin, lonid, lon
+ncdf_varget, idin, latid, lat
+ncdf_varget, idin, lonbndid, lonbnds
+ncdf_varget, idin, latbndid, latbnds
+ncdf_varget, idin, timbndid, timbnds_in
+ncdf_close, idin
+
+timbnds_out = make_array(2,1032)
+timbnds_out(*,0) = [0,31]
+for cpt=1,1031 do timbnds_out(0,cpt) = timbnds_out(0,cpt-1) + une_annee(cpt mod 12)
+for cpt=1,1030 do timbnds_out(1,cpt) = timbnds_out(0,cpt+1)
+timbnds_out(1,1031) = 31390
+
+emiss_interpol=make_array(nlon,nlat,1032,value=0.)
+
+indices_out_mois = make_array(86,12,value=0)
+for m=0,11 do begin
+   for n=0,85 do indices_out_mois(n,m) = 12*n + m
+endfor
+
+jours_out_mois = make_array(86,12,value=0)
+for n=0,85 do begin
+   for m=0,11 do jours_out_mois(n,m) = jours_out(indices_out_mois(n,m))
+endfor
+
+jours_in_mois = make_array(n_elements(jours_in)/12,12,value=0)
+for m=0,11 do begin
+   for n=0,n_elements(jours_in)/12-1 do jours_in_mois(n,m) = jours_in(n*12 + m)
+endfor
+
+indices_in_mois = make_array(n_elements(jours_in)/12,12,value=0)
+for m=0,11 do begin
+   for n=0,n_elements(jours_in)/12-1 do indices_in_mois(n,m) = where(jours_in eq jours_in_mois(n,m))
+endfor
+
+; Interpolation mois par mois...
+for i=0,nlon-1 do begin
+   for j=0,nlat-1 do begin
+      for m=0,11 do begin
+         emiss_interpol(i,j,indices_out_mois(*,m)) = interpol(emiss_init(i,j,indices_in_mois(*,m)),jours_in_mois(*,m),jours_out_mois(*,m))
+      endfor
+   endfor
+endfor
+
+plot, jours_out, emiss_interpol(404,192,*)
+
+idout = ncdf_create('BC-em-bb_interp.nc', /clobber)
+dimlonid = ncdf_dimdef(idout,'lon',nlon)
+dimlatid = ncdf_dimdef(idout,'lat',nlat)
+dimtimid = ncdf_dimdef(idout,'time',1032)
+dimbndid = ncdf_dimdef(idout,'bound',2)
+varlonid = ncdf_vardef(idout, 'lon', [dimlonid], /double)
+varlatid = ncdf_vardef(idout, 'lat', [dimlatid], /double)
+vartimid = ncdf_vardef(idout, 'time', [dimtimid], /double)
+varlonbndsid = ncdf_vardef(idout, 'lon_bnds', [dimbndid,dimlonid], /double)
+varlatbndsid = ncdf_vardef(idout, 'lat_bnds', [dimbndid,dimlatid], /double)
+vartimbndsid = ncdf_vardef(idout, 'time_bnds', [dimbndid,dimtimid], /double)
+varemissid = ncdf_vardef(idout, 'BC_em_openburning', [dimlonid,dimlatid,dimtimid], /float)
+ncdf_control, idout, /endef
+ncdf_varput, idout, varlonid, lon
+ncdf_varput, idout, varlatid, lat
+ncdf_varput, idout, varlonbndsid, lonbnds
+ncdf_varput, idout, varlatbndsid, latbnds
+ncdf_varput, idout, vartimbndsid, timbnds_out
+ncdf_varput, idout, varemissid, emiss_interpol
+ncdf_varput, idout, vartimid, jours_out
+ncdf_control, idout, /redef
+ncdf_attput, idout, varlonid, "units", "degrees_east"
+ncdf_attput, idout, varlonid, "long_name", "longitude"
+ncdf_attput, idout, varlonid, "axis", "X"
+ncdf_attput, idout, varlonid, "bounds", "lon_bnds"
+ncdf_attput, idout, varlonid, "modulo", 360.
+ncdf_attput, idout, varlonid, "realtopology", "circular"
+ncdf_attput, idout, varlonid, "standard_name", "longitude"
+ncdf_attput, idout, varlonid, "topology", "circular"
+ncdf_attput, idout, varlatid, "units", "degrees_north"
+ncdf_attput, idout, varlatid, "long_name", "latitude"
+ncdf_attput, idout, varlatid, "axis", "Y"
+ncdf_attput, idout, varlatid, "bounds", "lat_bnds"
+ncdf_attput, idout, varlatid, "realtopology", "linear"
+ncdf_attput, idout, varlatid, "standard_name", "latitude"
+ncdf_attput, idout, varemissid, "units", string(emiss_unit)
+ncdf_attput, idout, varemissid, "_FillValue", 1.e+20
+ncdf_attput, idout, varemissid, "cell_methods", "time: mean"
+ncdf_attput, idout, varemissid, "long_name", string(emiss_longname)
+ncdf_attput, idout, varemissid, "missing_value", 1.e+20
+ncdf_attput, idout, vartimid, "units", "days since 2015-01-01 0:0:0"
+ncdf_attput, idout, vartimid, "long_name", "time"
+ncdf_attput, idout, vartimid, "calendar", '365_day'
+ncdf_attput, idout, vartimid, "axis", "T"
+ncdf_attput, idout, vartimid, "bounds", "time_bnds"
+ncdf_attput, idout, vartimid, "realtopology", "linear"
+ncdf_attput, idout, vartimid, "standard_name", "time"
+idin = ncdf_open('/home/tlurton/Scen_EMISS_AER/REGRID/BC-em-bb_nc3.nc') ; Je rouvre le fichier
+res = ncdf_attcopy(idin, "_NCProperties", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "Conventions", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "activity_id", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "comment", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "contact", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "creation_date", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "data_structure", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "dataset_category", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "external_variables", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "frequency", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "further_info_url", idout, /in_global, /out_global)
+ncdf_attput, idout, /global, "grid", '2.5 x 1.25 degree longitude x latitude'
+res = ncdf_attcopy(idin, "grid_label", idout, /in_global, /out_global)
+ncdf_attput, idout, /global, "nominal_resolution", '250 x 125 km longitude x latitude'
+res = ncdf_attcopy(idin, "institution", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "institution_id", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "mip_era", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "product", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "realm", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "references", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "source", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "source_id", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "table_id", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "target_mip", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "title", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "variable_id", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "license", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "data_usage_tips", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "reporting_unit", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "tracking_id", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "history", idout, /in_global, /out_global)
+res = ncdf_attcopy(idin, "NCO", idout, /in_global, /out_global)
+ncdf_close, idin
+ncdf_close, idout
+
+stop
+
+end
Index: /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/netcdf.pro
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/netcdf.pro	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/netcdf.pro	(revision 4007)
@@ -0,0 +1,73 @@
+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/SCENARIOS/AER_TROP_EMISSIONS/REGRID/netcdfwrite.pro
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/netcdfwrite.pro	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/netcdfwrite.pro	(revision 4007)
@@ -0,0 +1,259 @@
+; 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/SCENARIOS/AER_TROP_EMISSIONS/REGRID/regrid_scen.job
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/regrid_scen.job	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/regrid_scen.job	(revision 4007)
@@ -0,0 +1,11 @@
+#PBS -N make_AER_Sc_v1
+#PBS -S /bin/bash
+#PBS -q week # there exist: short, day, days3, week...
+#PBS -k eo  # to write the output of stdin
+#PBS -l vmem=10gb   # virtual memory
+#PBS -l mem=10gb
+
+#replace the following line with your own running directory
+cd /home/tlurton/Scen_EMISS_AER/REGRID
+
+./regrid_scen.sh
Index: /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/regrid_scen.sh
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/regrid_scen.sh	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/regrid_scen.sh	(revision 4007)
@@ -0,0 +1,869 @@
+#--reworked 7/2018 for scenario emissions, ThL
+#--updated on 5/2/2018 with new paths, ThL
+#--updated on 4/5/2017 with improved Sheng & Zwiers algorithm, ThL
+#--corrected some interpolation preprocessing (compared to v4)
+#--updated on 5/5/2017 with correction factor (30/46) on NOx vs. NO
+#--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
+#--corrected on 26/09/2017 for NOx units: PNNL dataset is kg NO2, VUA is kg NO. 
+#  INCA expects NO in its AER version.
+#
+#INCA  conc_dms flx_nox flx_bc flx_pom flx_bbbc flx_bbpom flx_so2 flx_so4 flx_h2s flx_nh3
+#CMIP6 species  NOx     BC     OC                         SO2                     NH3
+# + NMVOC CO
+
+#--INCA example file where dms_conc can be reused
+fileINCAex="/home/oboucher/CMIP6/AER_EMISSIONS/INCAfile/sflx_lmdz_phy_1997.nc"
+
+#--input directory for emissions
+dirin="/prodigfs/project/input4MIPs/CMIP6/ScenarioMIP/IAMC/"
+
+#--LMDz grid information
+grid="144x143"
+gridfile="../GRID/grid-lmdz-lonlat_"${grid}
+nbpoint=$((144*141+2))
+
+for scen in "ssp119" "ssp126" "ssp245" "ssp370" "ssp434" "ssp460" "ssp534-over" "ssp585"
+do
+
+#--output directory
+dirout="/data/"${USER}"/CMIP6/AEROSOL/ScenarioMIP/${scen}/v1/"
+if [ ! -d ${dirout} ] ; then mkdir -p ${dirout} ; fi
+
+if [ ${scen} == "ssp119" ] ; then prefix="IMAGE" ; fi
+if [ ${scen} == "ssp126" ] ; then prefix="IMAGE" ; fi
+if [ ${scen} == "ssp245" ] ; then prefix="MESSAGE-GLOBIOM" ; fi
+if [ ${scen} == "ssp370" ] ; then prefix="AIM" ; fi
+if [ ${scen} == "ssp434" ] ; then prefix="GCAM4" ; fi
+if [ ${scen} == "ssp460" ] ; then prefix="GCAM4" ; fi
+if [ ${scen} == "ssp534-over" ] ; then prefix="REMIND-MAGPIE" ; fi
+if [ ${scen} == "ssp585" ] ; then prefix="REMIND-MAGPIE" ; fi
+
+#--Unique file for the whole period, so fixed variables:
+year1=2015
+year2=2100
+
+#--Loop on species (anthro)
+for species in "BC" "NOx" "OC" "SO2" "NH3"
+do
+
+echo "... "${year}" : Dealing with "${species}" anthro input file..."
+
+if [ ${species} = "BC"  ] ; then speciesinca="bc"  ; fi
+if [ ${species} = "NOx" ] ; then speciesinca="nox" ; fi
+if [ ${species} = "OC"  ] ; then speciesinca="pom" ; fi
+if [ ${species} = "SO2" ] ; then speciesinca="so2" ; fi
+if [ ${species} = "NH3" ] ; then speciesinca="nh3" ; fi
+
+if [ ${species} = "BC"  ] ; then speciesUp="BC"  ; fi
+if [ ${species} = "NOx" ] ; then speciesUp="NOX" ; fi
+if [ ${species} = "OC"  ] ; then speciesUp="OC" ; fi
+if [ ${species} = "SO2" ] ; then speciesUp="SO2" ; fi
+if [ ${species} = "NH3" ] ; then speciesUp="NH3" ; fi
+
+#--input file anthro
+filename=${dirin}IAMC-${prefix}-${scen}-1-1/atmos/mon/${species}_em_anthro/gn/v20180628/${species}-em-anthro_input4MIPs_emissions_ScenarioMIP_IAMC-${prefix}-${scen}-1-1_gn_${year1}01-${year2}12.nc
+
+#--input file summed over sectors
+filesum=${dirout}${species}-em-anthro_input4MIPs_emissions_ScenarioMIP_IAMC-${prefix}-${scen}-1-1_gn_${year1}01-${year2}12_sum.nc
+#--reworked input file
+filetmp=${dirout}${species}-em-anthro_input4MIPs_emissions_ScenarioMIP_IAMC-${prefix}-${scen}-1-1_gn_${year1}01-${year2}12_tmp.nc
+#--time-interpolated input file
+fileintanthro=${dirout}${species}-em-anthro_input4MIPs_emissions_ScenarioMIP_IAMC-${prefix}-${scen}-1-1_gn_${year1}01-${year2}12_interp.nc
+
+if [ -f ${fileintanthro} ]
+then
+   echo "...Nothing to do!"
+else
+
+rm -f ${filesum} ${filetmp} ${fileintanthro}
+
+#--re-formatting the input file
+#--summing over sectors
+rm -f sumsectors.jnl
+cat << eod > sumsectors.jnl
+set memory/size=100
+use "${filename}"
+save/clobber/file="${filesum}" ${species}_EM_ANTHRO[k=@sum]
+eod
+#--run ferret script
+ferret << eod 
+go sumsectors.jnl
+exit
+eod
+
+#--reformatting input file for IDL to be able to read
+ncks -3 ${filesum} ${filetmp}
+#--cleaning up
+rm -f ${filesum}
+
+#--IDL code for time interpolation, to fill in 9-year gaps in input file
+rm -f time_interp_over_gaps.pro
+cat << eod >> time_interp_over_gaps.pro
+pro time_interp_over_gaps
+
+; Les dimensions des fichiers
+nlon=720
+nlat=360
+nt=1032
+
+une_annee=[31,30,29,31,30,31,30,31,31,30,31,30]
+les_mois = [31,28,31,30,31,30,31,31,30,31,30,31]
+
+jours_in = [15, 45, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349, 1840, 1870, $
+    1899, 1930, 1960, 1991, 2021, 2052, 2083, 2113, 2144, 2174, 5490, 5520, $
+    5549, 5580, 5610, 5641, 5671, 5702, 5733, 5763, 5794, 5824, 9140, 9170, $
+    9199, 9230, 9260, 9291, 9321, 9352, 9383, 9413, 9444, 9474, 12790, 12820, $
+    12849, 12880, 12910, 12941, 12971, 13002, 13033, 13063, 13094, 13124, $
+    16440, 16470, 16499, 16530, 16560, 16591, 16621, 16652, 16683, 16713, $
+    16744, 16774, 20090, 20120, 20149, 20180, 20210, 20241, 20271, 20302, $
+    20333, 20363, 20394, 20424, 23740, 23770, 23799, 23830, 23860, 23891, $
+    23921, 23952, 23983, 24013, 24044, 24074, 27390, 27420, 27449, 27480, $ 
+    27510, 27541, 27571, 27602, 27633, 27663, 27694, 27724, 31040, 31070, $ 
+    31099, 31130, 31160, 31191, 31221, 31252, 31283, 31313, 31344, 31374]
+
+jours_out=make_array(1032)
+jours_out(0) = 15
+for cpt=1,1031 do jours_out(cpt) = jours_out(cpt-1) + une_annee(cpt mod 12)
+
+idin = ncdf_open('${filetmp}')
+varid = ncdf_varid(idin,'${speciesUp}_EM_ANTHRO')
+lonid = ncdf_varid(idin,'LON')
+latid = ncdf_varid(idin,'LAT')
+;lonbndid = ncdf_varid(idin,'LON_bnds')
+;latbndid = ncdf_varid(idin,'LAT_bnds')
+timbndid = ncdf_varid(idin,'TIME_bnds')
+ncdf_varget, idin, varid, emiss_init
+ncdf_attget, idin, varid, "units", emiss_unit
+ncdf_attget, idin, varid, "long_name", emiss_longname
+ncdf_varget, idin, lonid, lon
+ncdf_varget, idin, latid, lat
+;ncdf_varget, idin, lonbndid, lonbnds
+;ncdf_varget, idin, latbndid, latbnds
+ncdf_varget, idin, timbndid, timbnds_in
+ncdf_close, idin
+
+timbnds_out = make_array(2,1032)
+timbnds_out(*,0) = [0,31]
+timbnds_out(*,1) = [31,59]
+for cpt=2,1031 do timbnds_out(0,cpt) = timbnds_out(0,cpt-1) + les_mois((cpt-1) mod 12)
+for cpt=1,1030 do timbnds_out(1,cpt) = timbnds_out(0,cpt+1)
+timbnds_out(1,1031) = 31390
+
+emiss_interpol=make_array(nlon,nlat,1032,value=0.)
+
+indices_out_mois = make_array(86,12,value=0)
+for m=0,11 do begin
+   for n=0,85 do indices_out_mois(n,m) = 12*n + m
+endfor
+
+jours_out_mois = make_array(86,12,value=0)
+for n=0,85 do begin
+   for m=0,11 do jours_out_mois(n,m) = jours_out(indices_out_mois(n,m))
+endfor
+
+jours_in_mois = make_array(n_elements(jours_in)/12,12,value=0)
+for m=0,11 do begin
+   for n=0,n_elements(jours_in)/12-1 do jours_in_mois(n,m) = jours_in(n*12 + m)
+endfor
+
+indices_in_mois = make_array(n_elements(jours_in)/12,12,value=0)
+for m=0,11 do begin
+   for n=0,n_elements(jours_in)/12-1 do indices_in_mois(n,m) = where(jours_in eq jours_in_mois(n,m))
+endfor
+
+; Interpolation mois par mois...
+for i=0,nlon-1 do begin
+   for j=0,nlat-1 do begin
+      for m=0,11 do begin
+         emiss_interpol(i,j,indices_out_mois(*,m)) = interpol(emiss_init(i,j,indices_in_mois(*,m)),jours_in_mois(*,m),jours_out_mois(*,m))
+      endfor
+   endfor
+endfor
+
+idout = ncdf_create('${fileintanthro}', /clobber)
+dimlonid = ncdf_dimdef(idout,'lon',nlon)
+dimlatid = ncdf_dimdef(idout,'lat',nlat)
+dimtimid = ncdf_dimdef(idout,'time',1032)
+dimbndid = ncdf_dimdef(idout,'bound',2)
+varlonid = ncdf_vardef(idout, 'lon', [dimlonid], /double)
+varlatid = ncdf_vardef(idout, 'lat', [dimlatid], /double)
+vartimid = ncdf_vardef(idout, 'time', [dimtimid], /double)
+;varlonbndsid = ncdf_vardef(idout, 'lon_bnds', [dimbndid,dimlonid], /double)
+;varlatbndsid = ncdf_vardef(idout, 'lat_bnds', [dimbndid,dimlatid], /double)
+vartimbndsid = ncdf_vardef(idout, 'time_bnds', [dimbndid,dimtimid], /double)
+varemissid = ncdf_vardef(idout, '${species}_em_anthro', [dimlonid,dimlatid,dimtimid], /float)
+ncdf_control, idout, /endef
+ncdf_varput, idout, varlonid, lon
+ncdf_varput, idout, varlatid, lat
+;ncdf_varput, idout, varlonbndsid, lonbnds
+;ncdf_varput, idout, varlatbndsid, latbnds
+ncdf_varput, idout, vartimbndsid, timbnds_out
+ncdf_varput, idout, varemissid, emiss_interpol
+ncdf_varput, idout, vartimid, jours_out
+ncdf_control, idout, /redef
+ncdf_attput, idout, varlonid, "units", "degrees_east"
+ncdf_attput, idout, varlonid, "long_name", "longitude"
+ncdf_attput, idout, varlonid, "axis", "X"
+ncdf_attput, idout, varlonid, "bounds", "lon_bnds"
+ncdf_attput, idout, varlonid, "modulo", 360.
+ncdf_attput, idout, varlonid, "realtopology", "circular"
+ncdf_attput, idout, varlonid, "standard_name", "longitude"
+ncdf_attput, idout, varlonid, "topology", "circular"
+ncdf_attput, idout, varlatid, "units", "degrees_north"
+ncdf_attput, idout, varlatid, "long_name", "latitude"
+ncdf_attput, idout, varlatid, "axis", "Y"
+ncdf_attput, idout, varlatid, "bounds", "lat_bnds"
+ncdf_attput, idout, varlatid, "realtopology", "linear"
+ncdf_attput, idout, varlatid, "standard_name", "latitude"
+ncdf_attput, idout, varemissid, "units", string(emiss_unit)
+ncdf_attput, idout, varemissid, "_FillValue", 1.e+20
+ncdf_attput, idout, varemissid, "cell_methods", "time: mean"
+ncdf_attput, idout, varemissid, "long_name", string(emiss_longname)
+ncdf_attput, idout, varemissid, "missing_value", 1.e+20
+ncdf_attput, idout, vartimid, "units", "days since 2015-01-01 0:0:0"
+ncdf_attput, idout, vartimid, "long_name", "time"
+ncdf_attput, idout, vartimid, "calendar", '365_day'
+ncdf_attput, idout, vartimid, "axis", "T"
+ncdf_attput, idout, vartimid, "bounds", "time_bnds"
+ncdf_attput, idout, vartimid, "realtopology", "linear"
+ncdf_attput, idout, vartimid, "standard_name", "time"
+ncdf_close, idout
+end
+eod
+#
+#--calling IDL
+#
+/opt/idl-6.4/idl/bin/idl << eod
+.r time_interp_over_gaps.pro
+time_interp_over_gaps
+exit
+eod
+#
+
+#--adding lat_bnds and lon_bnds variables
+ncks -A -v lat_bnds,lon_bnds ${filename} ${fileintanthro}
+
+#--cleaning up
+rm -f ${filetmp}
+
+#--end test whether file tmp already exists
+fi
+
+
+
+#--now dealing with the BB input file
+
+echo "... "${year}" : Dealing with "${species}" openburning input file..."
+
+#--input file BB
+filename=${dirin}IAMC-${prefix}-${scen}-1-1/atmos/mon/${species}_em_openburning/gn/v20180628/${species}-em-openburning_input4MIPs_emissions_ScenarioMIP_IAMC-${prefix}-${scen}-1-1_gn_${year1}01-${year2}12.nc
+#--reworked input file
+filetmp=${dirout}${species}-em-openburning_input4MIPs_emissions_ScenarioMIP_IAMC-${prefix}-${scen}-1-1_gn_${year1}01-${year2}12_tmp.nc
+#--time-interpolated input file
+fileintbb=${dirout}${species}-em-openburning_input4MIPs_emissions_ScenarioMIP_IAMC-${prefix}-${scen}-1-1_gn_${year1}01-${year2}12_interp.nc
+
+if [ -f ${fileintbb} ]
+then
+   echo "...Nothing to do!"
+else
+
+rm -f ${filetmp} ${fileintbb}
+
+#--no need to sum over sectors, as none in BB files
+
+#--reformatting input file for IDL to be able to read
+ncks -3 ${filename} ${filetmp}
+
+#--IDL code for time interpolation, to fill in 9-year gaps in input file
+rm -f time_interp_over_gaps.pro
+cat << eod >> time_interp_over_gaps.pro
+pro time_interp_over_gaps
+
+; Les dimensions des fichiers
+nlon=720
+nlat=360
+nt=1032
+
+une_annee=[31,30,29,31,30,31,30,31,31,30,31,30]
+les_mois = [31,28,31,30,31,30,31,31,30,31,30,31]
+
+jours_in = [15, 45, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349, 1840, 1870, $
+    1899, 1930, 1960, 1991, 2021, 2052, 2083, 2113, 2144, 2174, 5490, 5520, $
+    5549, 5580, 5610, 5641, 5671, 5702, 5733, 5763, 5794, 5824, 9140, 9170, $
+    9199, 9230, 9260, 9291, 9321, 9352, 9383, 9413, 9444, 9474, 12790, 12820, $
+    12849, 12880, 12910, 12941, 12971, 13002, 13033, 13063, 13094, 13124, $
+    16440, 16470, 16499, 16530, 16560, 16591, 16621, 16652, 16683, 16713, $
+    16744, 16774, 20090, 20120, 20149, 20180, 20210, 20241, 20271, 20302, $
+    20333, 20363, 20394, 20424, 23740, 23770, 23799, 23830, 23860, 23891, $
+    23921, 23952, 23983, 24013, 24044, 24074, 27390, 27420, 27449, 27480, $ 
+    27510, 27541, 27571, 27602, 27633, 27663, 27694, 27724, 31040, 31070, $ 
+    31099, 31130, 31160, 31191, 31221, 31252, 31283, 31313, 31344, 31374]
+
+jours_out=make_array(1032)
+jours_out(0) = 15
+for cpt=1,1031 do jours_out(cpt) = jours_out(cpt-1) + une_annee(cpt mod 12)
+
+idin = ncdf_open('${filetmp}')
+varid = ncdf_varid(idin,'${species}_em_openburning')      ; Here, small caps
+lonid = ncdf_varid(idin,'lon')
+latid = ncdf_varid(idin,'lat')
+;lonbndid = ncdf_varid(idin,'lon_bnds')
+;latbndid = ncdf_varid(idin,'lat_bnds')
+timbndid = ncdf_varid(idin,'time_bnds')
+ncdf_varget, idin, varid, emiss_init
+ncdf_attget, idin, varid, "units", emiss_unit
+ncdf_attget, idin, varid, "long_name", emiss_longname
+ncdf_varget, idin, lonid, lon
+ncdf_varget, idin, latid, lat
+;ncdf_varget, idin, lonbndid, lonbnds
+;ncdf_varget, idin, latbndid, latbnds
+ncdf_varget, idin, timbndid, timbnds_in
+ncdf_close, idin
+
+timbnds_out = make_array(2,1032)
+timbnds_out(*,0) = [0,31]
+timbnds_out(*,1) = [31,59]
+for cpt=2,1031 do timbnds_out(0,cpt) = timbnds_out(0,cpt-1) + une_annee((cpt-1) mod 12)
+for cpt=1,1030 do timbnds_out(1,cpt) = timbnds_out(0,cpt+1)
+timbnds_out(1,1031) = 31390
+
+emiss_interpol=make_array(nlon,nlat,1032,value=0.)
+
+indices_out_mois = make_array(86,12,value=0)
+for m=0,11 do begin
+   for n=0,85 do indices_out_mois(n,m) = 12*n + m
+endfor
+
+jours_out_mois = make_array(86,12,value=0)
+for n=0,85 do begin
+   for m=0,11 do jours_out_mois(n,m) = jours_out(indices_out_mois(n,m))
+endfor
+
+jours_in_mois = make_array(n_elements(jours_in)/12,12,value=0)
+for m=0,11 do begin
+   for n=0,n_elements(jours_in)/12-1 do jours_in_mois(n,m) = jours_in(n*12 + m)
+endfor
+
+indices_in_mois = make_array(n_elements(jours_in)/12,12,value=0)
+for m=0,11 do begin
+   for n=0,n_elements(jours_in)/12-1 do indices_in_mois(n,m) = where(jours_in eq jours_in_mois(n,m))
+endfor
+
+; Interpolation mois par mois...
+for i=0,nlon-1 do begin
+   for j=0,nlat-1 do begin
+      for m=0,11 do begin
+         emiss_interpol(i,j,indices_out_mois(*,m)) = interpol(emiss_init(i,j,indices_in_mois(*,m)),jours_in_mois(*,m),jours_out_mois(*,m))
+      endfor
+   endfor
+endfor
+
+idout = ncdf_create('${fileintbb}', /clobber)
+dimlonid = ncdf_dimdef(idout,'lon',nlon)
+dimlatid = ncdf_dimdef(idout,'lat',nlat)
+dimtimid = ncdf_dimdef(idout,'time',1032)
+dimbndid = ncdf_dimdef(idout,'bound',2)
+varlonid = ncdf_vardef(idout, 'lon', [dimlonid], /double)
+varlatid = ncdf_vardef(idout, 'lat', [dimlatid], /double)
+vartimid = ncdf_vardef(idout, 'time', [dimtimid], /double)
+;varlonbndsid = ncdf_vardef(idout, 'lon_bnds', [dimbndid,dimlonid], /double)
+;varlatbndsid = ncdf_vardef(idout, 'lat_bnds', [dimbndid,dimlatid], /double)
+vartimbndsid = ncdf_vardef(idout, 'time_bnds', [dimbndid,dimtimid], /double)
+varemissid = ncdf_vardef(idout, '${species}_em_openburning', [dimlonid,dimlatid,dimtimid], /float)
+ncdf_control, idout, /endef
+ncdf_varput, idout, varlonid, lon
+ncdf_varput, idout, varlatid, lat
+;ncdf_varput, idout, varlonbndsid, lonbnds
+;ncdf_varput, idout, varlatbndsid, latbnds
+ncdf_varput, idout, vartimbndsid, timbnds_out
+ncdf_varput, idout, varemissid, emiss_interpol
+ncdf_varput, idout, vartimid, jours_out
+ncdf_control, idout, /redef
+ncdf_attput, idout, varlonid, "units", "degrees_east"
+ncdf_attput, idout, varlonid, "long_name", "longitude"
+ncdf_attput, idout, varlonid, "axis", "X"
+ncdf_attput, idout, varlonid, "bounds", "lon_bnds"
+ncdf_attput, idout, varlonid, "modulo", 360.
+ncdf_attput, idout, varlonid, "realtopology", "circular"
+ncdf_attput, idout, varlonid, "standard_name", "longitude"
+ncdf_attput, idout, varlonid, "topology", "circular"
+ncdf_attput, idout, varlatid, "units", "degrees_north"
+ncdf_attput, idout, varlatid, "long_name", "latitude"
+ncdf_attput, idout, varlatid, "axis", "Y"
+ncdf_attput, idout, varlatid, "bounds", "lat_bnds"
+ncdf_attput, idout, varlatid, "realtopology", "linear"
+ncdf_attput, idout, varlatid, "standard_name", "latitude"
+ncdf_attput, idout, varemissid, "units", string(emiss_unit)
+ncdf_attput, idout, varemissid, "_FillValue", 1.e+20
+ncdf_attput, idout, varemissid, "cell_methods", "time: mean"
+ncdf_attput, idout, varemissid, "long_name", string(emiss_longname)
+ncdf_attput, idout, varemissid, "missing_value", 1.e+20
+ncdf_attput, idout, vartimid, "units", "days since 2015-01-01 0:0:0"
+ncdf_attput, idout, vartimid, "long_name", "time"
+ncdf_attput, idout, vartimid, "calendar", '365_day'
+ncdf_attput, idout, vartimid, "axis", "T"
+ncdf_attput, idout, vartimid, "bounds", "time_bnds"
+ncdf_attput, idout, vartimid, "realtopology", "linear"
+ncdf_attput, idout, vartimid, "standard_name", "time"
+ncdf_close, idout
+end
+eod
+#
+#--calling IDL
+#
+/opt/idl-6.4/idl/bin/idl << eod
+.r time_interp_over_gaps.pro
+time_interp_over_gaps
+exit
+eod
+#
+
+#--adding lat_bnds and lon_bnds variables
+ncks -A -v lat_bnds,lon_bnds ${filename} ${fileintbb}
+
+#--cleaning up
+rm -f ${filetmp}
+
+#--end test if interpolated file already exists
+fi
+
+#--end loop on species
+done
+
+
+
+#--anthro and BB input files fileintanthro and fileintbb are now ready, for the whole period and for the current scenario
+#--new double loop on species and years
+
+#--loop on years
+for year in {2015..2100}
+do
+
+#--loop on species
+for species in "BC" "NOx" "OC" "SO2" "NH3"
+do
+
+if [ ${species} = "BC"  ] ; then speciesinca="bc"  ; fi
+if [ ${species} = "NOx" ] ; then speciesinca="nox" ; fi
+if [ ${species} = "OC"  ] ; then speciesinca="pom" ; fi
+if [ ${species} = "SO2" ] ; then speciesinca="so2" ; fi
+if [ ${species} = "NH3" ] ; then speciesinca="nh3" ; fi
+
+#--dealing with anthro
+
+#--input file
+fileintanthro=${dirout}${species}-em-anthro_input4MIPs_emissions_ScenarioMIP_IAMC-${prefix}-${scen}-1-1_gn_${year1}01-${year2}12_interp.nc
+#--output files
+filenameout1=${dirout}flux_${speciesinca}_${year}.nc
+filenameout1b=${dirout}flux_0_${speciesinca}_${year}.nc
+filenameout2=${dirout}flux_lmdz_${speciesinca}_${year}.nc
+filenameout3=${dirout}flux_vector_${speciesinca}_${year}.nc
+
+echo ${fileintanthro} ${filenameout1} ${filenameout2} ${filenameout3}
+rm -f ${filenameout1} ${filenameout2} ${filenameout3}
+
+#--extracting the correct year
+rm -f select_yr.jnl
+cat << eod > select_yr.jnl
+set memory/size=100
+use "${fileintanthro}"
+set region/t=16-jan-${year}:16-dec-${year}
+save/clobber/file="${filenameout1}" ${species}_EM_ANTHRO
+eod
+#--run ferret script
+ferret << eod 
+go select_yr.jnl
+exit
+eod
+
+#--replace undef with 0
+cdo setmisstoc,0.0 ${filenameout1} ${filenameout1b}
+
+#--remap to LMDz grid
+#--OC to POM conversion factor
+#--otherwise change to capital letters if not (eg NOx)
+if [ ${species} == "OC"  ] ; then 
+echo cdo remapcon,${gridfile} -chname,`echo ${species}_EM_ANTHRO | awk '{print toupper($0)}'`,flux -mulc,1.4 ${filenameout1b} ${filenameout2}
+cdo remapcon,${gridfile} -chname,`echo ${species}_EM_ANTHRO | awk '{print toupper($0)}'`,flux -mulc,1.4 ${filenameout1b} ${filenameout2}
+else
+echo cdo remapcon,${gridfile} -chname,`echo ${species}_EM_ANTHRO | awk '{print toupper($0)}'`,flux  ${filenameout1b} ${filenameout2}
+cdo remapcon,${gridfile} -chname,`echo ${species}_EM_ANTHRO | awk '{print toupper($0)}'`,flux ${filenameout1b} ${filenameout2}
+fi
+
+#--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.]   ])
+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_${speciesinca}: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_${speciesinca}=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
+#
+
+#--cleaning the mess
+rm -f ferret*
+rm -f regrid.pro
+
+
+#--now dealing with BB
+
+#--input file
+fileintbb=${dirout}${species}-em-openburning_input4MIPs_emissions_ScenarioMIP_IAMC-${prefix}-${scen}-1-1_gn_${year1}01-${year2}12_interp.nc
+#--output files
+filenameout1=${dirout}flux_${speciesinca}bb_${year}.nc
+filenameout1b=${dirout}flux_0_${speciesinca}bb_${year}.nc
+filenameout2=${dirout}flux_lmdz0_${speciesinca}bb_${year}.nc
+filenameout3=${dirout}flux_vector_${speciesinca}bb_${year}.nc
+
+echo ${fileintbb} ${filenameout1} ${filenameout1b} ${filenameout2} ${filenameout3}
+rm -f ${filenameout1} ${filenameout1b} ${filenameout2} ${filenameout3}
+
+#--unfortunately idl not happy with input netcdf files so need to ferretize files
+#--I extract the correct year as well; no need to sum up sectors, as none in BB files
+rm -f rewrite.jnl
+cat << eod > rewrite.jnl
+set memory/size=100
+use "${fileintbb}"
+set region/t=16-jan-${year}:16-dec-${year}
+save/clobber/file="${filenameout1}" ${species}_EM_OPENBURNING
+eod
+#--run ferret script
+ferret << eod
+go rewrite.jnl
+exit
+eod
+
+#--replace undef with 0 
+cdo setmisstoc,0.0 ${filenameout1} ${filenameout1b}
+
+#--remap to LMDz grid
+#--OC to POM conversion factor
+#--as ferret returns NOX, treat NOx NOX inconsistency in names by converting to upper case
+if [ ${species} != "OC" ] ; then
+echo cdo remapcon,${gridfile} -chname,`echo ${species}_EM_OPENBURNING | awk '{print toupper($0)}'`,flux ${filenameout1b} ${filenameout2}
+cdo remapcon,${gridfile} -chname,`echo ${species}_EM_OPENBURNING | awk '{print toupper($0)}'`,flux ${filenameout1b} ${filenameout2}
+else
+echo cdo remapcon,${gridfile} -chname,${species}_EM_OPENBURNING,flux -mulc,1.4 ${filenameout1b} ${filenameout2}
+cdo remapcon,${gridfile} -chname,${species}_EM_OPENBURNING,flux -mulc,1.4 ${filenameout1b} ${filenameout2}
+fi
+
+#--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.]   ])
+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 (Ã  part 1 seule correction matricielle)
+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_bb${speciesinca}: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_bb${speciesinca}=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
+
+#--unfortunately idl use capital letters for variable names so need to change to small letters for now 
+rm -f ${dirout}flux_vector_h2s_${year}.nc ${dirout}flux_vector_so4_${year}.nc
+cdo expr,'flx_h2s=0.*FLX_SO2' ${dirout}flux_vector_so2_${year}.nc ${dirout}flux_vector_h2s_${year}.nc
+cdo expr,'flx_so4=0.*FLX_SO2' ${dirout}flux_vector_so2_${year}.nc ${dirout}flux_vector_so4_${year}.nc
+
+rm -f ${dirout}flux_vector_${year}.nc
+
+#--combining everything into a single file with some final preprocessing
+rm -f ${dirout}flux_vector_noxtot_${year}.nc ${dirout}flux_vector_so2tot_${year}.nc ${dirout}flux_vector_nh3tot_${year}.nc
+
+#--deleting output file if already there
+rm -f ${dirout}sflx_lmdz_cmip6_${year}.nc
+#--merging all files into a single one
+#--anthro NOx is NO2 so 30/46 scaling factor to change to NO 
+#--openburning NOx is NO so no change in unit
+cdo merge -expr,'flx_bc=FLX_BC' ${dirout}flux_vector_bc_${year}.nc -expr,'flx_bbbc=FLX_BBBC' ${dirout}flux_vector_bcbb_${year}.nc -expr,'flx_pom=FLX_POM' ${dirout}flux_vector_pom_${year}.nc -expr,'flx_bbpom=FLX_BBPOM' ${dirout}flux_vector_pombb_${year}.nc -expr,'flx_nox=30.*FLX_NOX/46.' ${dirout}flux_vector_nox_${year}.nc -expr,'flx_bbnox=FLX_BBNOX;' ${dirout}flux_vector_noxbb_${year}.nc -expr,'flx_so2=FLX_SO2' ${dirout}flux_vector_so2_${year}.nc -expr,'flx_bbso2=FLX_BBSO2' ${dirout}flux_vector_so2bb_${year}.nc -expr,'flx_nh3=FLX_NH3' ${dirout}flux_vector_nh3_${year}.nc -expr,'flx_bbnh3=FLX_BBNH3' ${dirout}flux_vector_nh3bb_${year}.nc  ${dirout}flux_vector_h2s_${year}.nc ${dirout}flux_vector_so4_${year}.nc -selname,conc_dms ${fileINCAex} ${dirout}sflx_lmdz_cmip6_${year}.nc
+
+ncrename -d VECTOR,vector -v VECTOR,vector ${dirout}sflx_lmdz_cmip6_${year}.nc
+ncrename -d TIME,time     -v TIME,time     ${dirout}sflx_lmdz_cmip6_${year}.nc
+
+#--cleaning up 
+#rm -f ${dirout}flux*_${year}.nc
+#--to be uncommented in final version
+
+#--end loop on years
+done
+
+rm -f ${fileintanthro} ${fileintbb}
+
+#--end loop on scenarios
+done
+
+#--cleaning the mess
+rm -f ferret*
+rm -f regrid.pro
+rm -f time_interp_over_gaps.pro sumsectors.jnl select_yr.jnl rewrite.jnl
Index: /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/struct_autoconvert.pro
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/struct_autoconvert.pro	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/struct_autoconvert.pro	(revision 4007)
@@ -0,0 +1,60 @@
+;-------------------------------------------------------------
+;+
+; 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/SCENARIOS/AER_TROP_EMISSIONS/REGRID/struct_dims.pro
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/struct_dims.pro	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/struct_dims.pro	(revision 4007)
@@ -0,0 +1,207 @@
+; 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
+
+
Index: /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/time_interp_over_gaps.pro
===================================================================
--- /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/time_interp_over_gaps.pro	(revision 4007)
+++ /TOOLS/CMIP6_FORCING/SCENARIOS/AER_TROP_EMISSIONS/REGRID/time_interp_over_gaps.pro	(revision 4007)
@@ -0,0 +1,130 @@
+pro time_interp_over_gaps
+
+; Les dimensions des fichiers
+nlon=720
+nlat=360
+nt=1032
+
+une_annee=[31,30,29,31,30,31,30,31,31,30,31,30]
+les_mois = [31,28,31,30,31,30,31,31,30,31,30,31]
+
+jours_in = [15, 45, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349, 1840, 1870, $
+    1899, 1930, 1960, 1991, 2021, 2052, 2083, 2113, 2144, 2174, 5490, 5520, $
+    5549, 5580, 5610, 5641, 5671, 5702, 5733, 5763, 5794, 5824, 9140, 9170, $
+    9199, 9230, 9260, 9291, 9321, 9352, 9383, 9413, 9444, 9474, 12790, 12820, $
+    12849, 12880, 12910, 12941, 12971, 13002, 13033, 13063, 13094, 13124, $
+    16440, 16470, 16499, 16530, 16560, 16591, 16621, 16652, 16683, 16713, $
+    16744, 16774, 20090, 20120, 20149, 20180, 20210, 20241, 20271, 20302, $
+    20333, 20363, 20394, 20424, 23740, 23770, 23799, 23830, 23860, 23891, $
+    23921, 23952, 23983, 24013, 24044, 24074, 27390, 27420, 27449, 27480, $ 
+    27510, 27541, 27571, 27602, 27633, 27663, 27694, 27724, 31040, 31070, $ 
+    31099, 31130, 31160, 31191, 31221, 31252, 31283, 31313, 31344, 31374]
+
+jours_out=make_array(1032)
+jours_out(0) = 15
+for cpt=1,1031 do jours_out(cpt) = jours_out(cpt-1) + une_annee(cpt mod 12)
+
+idin = ncdf_open('/data/tlurton/CMIP6/AEROSOL/ScenarioMIP/ssp119/v1/NH3-em-openburning_input4MIPs_emissions_ScenarioMIP_IAMC-IMAGE-ssp119-1-1_gn_201501-210012_tmp.nc')
+varid = ncdf_varid(idin,'NH3_em_openburning')      ; Here, small caps
+lonid = ncdf_varid(idin,'lon')
+latid = ncdf_varid(idin,'lat')
+;lonbndid = ncdf_varid(idin,'lon_bnds')
+;latbndid = ncdf_varid(idin,'lat_bnds')
+timbndid = ncdf_varid(idin,'time_bnds')
+ncdf_varget, idin, varid, emiss_init
+ncdf_attget, idin, varid, "units", emiss_unit
+ncdf_attget, idin, varid, "long_name", emiss_longname
+ncdf_varget, idin, lonid, lon
+ncdf_varget, idin, latid, lat
+;ncdf_varget, idin, lonbndid, lonbnds
+;ncdf_varget, idin, latbndid, latbnds
+ncdf_varget, idin, timbndid, timbnds_in
+ncdf_close, idin
+
+timbnds_out = make_array(2,1032)
+timbnds_out(*,0) = [0,31]
+timbnds_out(*,1) = [31,59]
+for cpt=2,1031 do timbnds_out(0,cpt) = timbnds_out(0,cpt-1) + une_annee((cpt-1) mod 12)
+for cpt=1,1030 do timbnds_out(1,cpt) = timbnds_out(0,cpt+1)
+timbnds_out(1,1031) = 31390
+
+emiss_interpol=make_array(nlon,nlat,1032,value=0.)
+
+indices_out_mois = make_array(86,12,value=0)
+for m=0,11 do begin
+   for n=0,85 do indices_out_mois(n,m) = 12*n + m
+endfor
+
+jours_out_mois = make_array(86,12,value=0)
+for n=0,85 do begin
+   for m=0,11 do jours_out_mois(n,m) = jours_out(indices_out_mois(n,m))
+endfor
+
+jours_in_mois = make_array(n_elements(jours_in)/12,12,value=0)
+for m=0,11 do begin
+   for n=0,n_elements(jours_in)/12-1 do jours_in_mois(n,m) = jours_in(n*12 + m)
+endfor
+
+indices_in_mois = make_array(n_elements(jours_in)/12,12,value=0)
+for m=0,11 do begin
+   for n=0,n_elements(jours_in)/12-1 do indices_in_mois(n,m) = where(jours_in eq jours_in_mois(n,m))
+endfor
+
+; Interpolation mois par mois...
+for i=0,nlon-1 do begin
+   for j=0,nlat-1 do begin
+      for m=0,11 do begin
+         emiss_interpol(i,j,indices_out_mois(*,m)) = interpol(emiss_init(i,j,indices_in_mois(*,m)),jours_in_mois(*,m),jours_out_mois(*,m))
+      endfor
+   endfor
+endfor
+
+idout = ncdf_create('/data/tlurton/CMIP6/AEROSOL/ScenarioMIP/ssp119/v1/NH3-em-openburning_input4MIPs_emissions_ScenarioMIP_IAMC-IMAGE-ssp119-1-1_gn_201501-210012_interp.nc', /clobber)
+dimlonid = ncdf_dimdef(idout,'lon',nlon)
+dimlatid = ncdf_dimdef(idout,'lat',nlat)
+dimtimid = ncdf_dimdef(idout,'time',1032)
+dimbndid = ncdf_dimdef(idout,'bound',2)
+varlonid = ncdf_vardef(idout, 'lon', [dimlonid], /double)
+varlatid = ncdf_vardef(idout, 'lat', [dimlatid], /double)
+vartimid = ncdf_vardef(idout, 'time', [dimtimid], /double)
+;varlonbndsid = ncdf_vardef(idout, 'lon_bnds', [dimbndid,dimlonid], /double)
+;varlatbndsid = ncdf_vardef(idout, 'lat_bnds', [dimbndid,dimlatid], /double)
+vartimbndsid = ncdf_vardef(idout, 'time_bnds', [dimbndid,dimtimid], /double)
+varemissid = ncdf_vardef(idout, 'NH3_em_openburning', [dimlonid,dimlatid,dimtimid], /float)
+ncdf_control, idout, /endef
+ncdf_varput, idout, varlonid, lon
+ncdf_varput, idout, varlatid, lat
+;ncdf_varput, idout, varlonbndsid, lonbnds
+;ncdf_varput, idout, varlatbndsid, latbnds
+ncdf_varput, idout, vartimbndsid, timbnds_out
+ncdf_varput, idout, varemissid, emiss_interpol
+ncdf_varput, idout, vartimid, jours_out
+ncdf_control, idout, /redef
+ncdf_attput, idout, varlonid, "units", "degrees_east"
+ncdf_attput, idout, varlonid, "long_name", "longitude"
+ncdf_attput, idout, varlonid, "axis", "X"
+ncdf_attput, idout, varlonid, "bounds", "lon_bnds"
+ncdf_attput, idout, varlonid, "modulo", 360.
+ncdf_attput, idout, varlonid, "realtopology", "circular"
+ncdf_attput, idout, varlonid, "standard_name", "longitude"
+ncdf_attput, idout, varlonid, "topology", "circular"
+ncdf_attput, idout, varlatid, "units", "degrees_north"
+ncdf_attput, idout, varlatid, "long_name", "latitude"
+ncdf_attput, idout, varlatid, "axis", "Y"
+ncdf_attput, idout, varlatid, "bounds", "lat_bnds"
+ncdf_attput, idout, varlatid, "realtopology", "linear"
+ncdf_attput, idout, varlatid, "standard_name", "latitude"
+ncdf_attput, idout, varemissid, "units", string(emiss_unit)
+ncdf_attput, idout, varemissid, "_FillValue", 1.e+20
+ncdf_attput, idout, varemissid, "cell_methods", "time: mean"
+ncdf_attput, idout, varemissid, "long_name", string(emiss_longname)
+ncdf_attput, idout, varemissid, "missing_value", 1.e+20
+ncdf_attput, idout, vartimid, "units", "days since 2015-01-01 0:0:0"
+ncdf_attput, idout, vartimid, "long_name", "time"
+ncdf_attput, idout, vartimid, "calendar", '365_day'
+ncdf_attput, idout, vartimid, "axis", "T"
+ncdf_attput, idout, vartimid, "bounds", "time_bnds"
+ncdf_attput, idout, vartimid, "realtopology", "linear"
+ncdf_attput, idout, vartimid, "standard_name", "time"
+ncdf_close, idout
+end
