Index: TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/README
===================================================================
--- TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/README	(revision 6267)
+++ TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/README	(revision 6267)
@@ -0,0 +1,69 @@
+Stratospheric aerosol files for LMDz can be created by running on ciclad the script 
+./volc_ECRAD.sh
+which itself prepares two more files and submit a job on ciclad which may take a couple of days to run.
+
+The stratospheric AOD data originates from data prepared by Beiping Luo ETZH for CMIP6. 
+See http://www.wcrp-climate.org/wgcm-cmip/wgcm-cmip6
+and ftp://iacftp.ethz.ch/pub_read/luo/CMIP6/
+The data were provided in June 2016 and were not versioned (but were later referred as v2). 
+A new version (v3) was made available on 15 September 2017 due to a glitch in the 1970s data. 
+The new version also modified the trend in the background aerosol.
+
+Beiping calculated for the solar radiation 8 bands (bands 1-6 for the SW 6 bands 
+model and 7-8 for the SW 2 bands model) and likewise for the 16 bands of the RRTM LW model. 
+
+The data are processed in IDL on ciclad for different resolutions 
+of the IPSL model and aerosol optical properties are generated for AOD, 
+omega and g per model layer (SW), and absorption AOD per model layer (LW).
+
+The stratospheric aerosol data extends into the troposphere but a mask is applied 
+in the LMDz model to remove tropospheric data. The data are interpolated on the vertical 
+from the CMIP6 altitude coordinate (z, resolution 0.5 km) to the model coordinate (sigma). 
+This is done using a monthly (10-year mean) climatology of the geopotential height. 
+zalt_LR_l79.nc = 10 yr climatology from a coupled IPSL-CM6.0.12 simulation
+zalt_zonmean_LR_l79.nc = zonal mean of the above file  (cdo zonmean zalt.nc zalt_zonmean.nc)
+zalt_zonmean_LR_l79_rev.nc = same as above but ferret'ized so that can be read in idl
+
+In September 2017 the processing idl routine was changed to rely on climatological LMDz altitudes
+that vary monthly and with latitude (based on GEOP variable Z=GEOP/RG/1000 [in km]). 
+Start from a file geop.nc with monthly GEOP variable (ideally from a 10-yr run)
+cdo ymonmean geop.nc geop_LR.nc
+cdo divc,9810. geop_LR.nc zalt_LR_l79.nc
+cdo zonmean zalt_LR_l79.nc zalt_zonmean_LR_l79.nc
+ferret'ize file to make it readable by ferret : zalt_zonmean_LR_l79_rev.nc 
+for this an example is provided in example_convert.jnl
+We have files for VLR, LR, MR, HR model versions.
+The variable is called GEOP but is now the altitude in km as required by the interpolation scheme.
+
+The data are fed on ciclad (currently in /data/$USER/CMIP6/VOLC/)
+
+Annual files: taulwstrat.2D.$year.nc and taulwstrat.2D.$year.nc
+to be used for historical simulations
+
+Average 1850-2014 conditions (no seasonal cycle): taulwstrat.2D.ave.nc and tauswstrat.2D.ave.nc
+to be used for piControl, pdControl and for scenario runs (past a 10 year transition from 2014 conditions)
+subject to changes in the protocol
+
+Following the CMIP6 protocol, a linear transition between 2014 and average conditions is 
+performed and should be used for the period 2015-2023. Average conditions should then be used for
+2024 onwards.
+
+CMIP6 stratospheric aerosols are activated in LMDz with flag_aerosol_strat=2
+
+Files to be renamed as taulwstrat.2D.nc and tauswstrat.2D.nc for use by LMDz 
+No additional interpolation between months is done. 
+
+v3 : original CMIP6 files 
+v4 : corrected CMIP6 files (Pinatubo) and extended in time beyond 2014
+v4_ECRAD : adapted to deal with ECRAD files.
+
+Olivier Boucher 
+27/07/2016
+Revised 10/09/2017
+Revised 04/10/2017
+Revised 15/04/2018
+Revised 18/09/2018
+Revised 11/12/2018
+Revised 06/2019 by Th. Lurton to add support for VLR-L39.
+Revised 03/07/2019
+Modified 09/2022 by Th. Lurton to support ECRAD files.
Index: TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/example_convert.jnl
===================================================================
--- TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/example_convert.jnl	(revision 6267)
+++ TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/example_convert.jnl	(revision 6267)
@@ -0,0 +1,3 @@
+use zalt_zonmean_VLR_l39.nc
+save/file=zalt_zonmean_VLR_l39_rev.nc GEOP
+quit
Index: TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/grid_lr
===================================================================
--- TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/grid_lr	(revision 6267)
+++ TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/grid_lr	(revision 6267)
@@ -0,0 +1,13 @@
+gridtype  = lonlat
+gridsize  = 190
+xname     = lon
+xlongname = longitude
+xunits    = degrees_east
+yname     = LAT
+ylongname = latitude
+yunits    = degrees_north
+xsize     = 2
+xvals     = 0 180
+ysize     = 95
+yfirst    = 90
+yinc      = -1.91489
Index: TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/grid_vlr
===================================================================
--- TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/grid_vlr	(revision 6267)
+++ TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/grid_vlr	(revision 6267)
@@ -0,0 +1,14 @@
+gridtype  = lonlat
+gridsize  = 9216
+xname     = lon
+xlongname = Longitude
+xunits    = degrees_east
+yname     = lat
+ylongname = Latitude
+yunits    = degrees_north
+xsize     = 96
+ysize     = 96
+xfirst    = -180
+xinc      = 3.75
+yfirst    = 90
+yinc      = -1.89474
Index: TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/netcdf.pro
===================================================================
--- TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/netcdf.pro	(revision 6267)
+++ TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/netcdf.pro	(revision 6267)
@@ -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/AER_STRAT/ECRAD/netcdfwrite.pro
===================================================================
--- TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/netcdfwrite.pro	(revision 6267)
+++ TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/netcdfwrite.pro	(revision 6267)
@@ -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/AER_STRAT/ECRAD/struct_autoconvert.pro
===================================================================
--- TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/struct_autoconvert.pro	(revision 6267)
+++ TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/struct_autoconvert.pro	(revision 6267)
@@ -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/AER_STRAT/ECRAD/struct_dims.pro
===================================================================
--- TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/struct_dims.pro	(revision 6267)
+++ TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/struct_dims.pro	(revision 6267)
@@ -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
+
Index: TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/volc_ECRAD.sh
===================================================================
--- TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/volc_ECRAD.sh	(revision 6267)
+++ TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/volc_ECRAD.sh	(revision 6267)
@@ -0,0 +1,601 @@
+dirpwd=`pwd`
+
+#original CMIP6 files
+#vv='v3'
+#corrected CMIP6 files + extended in time
+vv='v4_ECRAD'
+
+#--choose the resolution
+#--only LR and zoom works fully for now !
+#lmdz='VVLR_l19'
+#lmdz='VLR_l39'
+#lmdz='VLR_l79'
+lmdz='LR_l79'
+#lmdz='MR_l79'
+#--added for Frederique Cheruy
+#lmdz='zoom_128x89_l79'
+
+dirout='/data/'${USER}'/CMIP6/VOLC/'${lmdz}'_'${vv}'/'
+echo $dirout
+if [ ! -d ${dirout} ] ; then mkdir -p ${dirout} ; fi
+
+cat > process_volc_${lmdz}.pro << EOF
+pro regrid
+;--version sans longitude
+;--latitude - level - wavelength - time
+;--sep 2017: add lat and month z variations of model layers
+;--Olivier Boucher
+;
+lmdz='${lmdz}'
+;
+if (lmdz eq 'VVLR_l19') then begin
+dimlonlmdz=48
+dimlatlmdz=37
+dimz=19
+output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
+latfirst=90.
+latinc=-180./float(dimlatlmdz-1)
+latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
+endif
+
+if (lmdz eq 'VLR_l39') then begin
+dimlonlmdz=96
+dimlatlmdz=96
+dimz=39
+output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
+latfirst=90.
+latinc=-180./float(dimlatlmdz-1)
+latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
+endif
+
+if (lmdz eq 'VLR_l79') then begin
+dimlonlmdz=96
+dimlatlmdz=96
+dimz=79
+output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
+latfirst=90.
+latinc=-180./float(dimlatlmdz-1)
+latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
+endif
+
+if (lmdz eq 'LR_l79') then begin
+dimlonlmdz=144
+dimlatlmdz=143
+dimz=79
+output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
+latfirst=90.
+latinc=-180./float(dimlatlmdz-1)
+latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
+endif
+
+if (lmdz eq 'MR_l79') then begin
+dimlonlmdz=256
+dimlatlmdz=257
+dimz=79
+output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
+latfirst=90.
+latinc=-180./float(dimlatlmdz-1)
+latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
+endif
+
+if (lmdz eq 'zoom_128x89_l79') then begin
+dimlonlmdz=128
+dimlatlmdz=89
+dimz=79
+output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
+latitudelmdz=[90., 86.85495, 83.8430634, 80.9792252, 78.2779922, 75.7526321, 73.4137878, 71.2680206, 69.3164978, $
+              67.5542145,   65.9702301, 64.5487747, 63.2709846, 62.1167107, 61.0660896, 60.1005592, 59.2034454, $ 
+              58.360157,    57.5581551, 56.7868576, 56.0374718, 55.302887,  54.5775719, 53.8574829, 53.1399384, $
+              52.4233971,   51.7071419, 50.9909439, 50.2747459, 49.5585518, 48.8423538, 48.1261597, 47.4099655, $
+              46.6937675,   45.9775734, 45.2613792, 44.5451813, 43.8289604, 43.1125679, 42.395504,  41.6765518, $ 
+              40.9534645,   40.2226982, 39.4792557, 38.7165909, 37.9265213, 37.099102,  36.2225075, 35.2828789, $ 
+              34.2642632,   33.1487236, 31.9167614, 30.5481758, 29.0234509, 27.32551,   25.4415455, 23.3643932, $ 
+              21.0930786,   18.6323872, 15.9917059, 13.1835403, 10.2221146, 7.12228632, 3.89880991, 0.565921485,$
+              -2.86286402, -6.37483072, -9.95808029, -13.6015415, -17.2949867, -21.029068, -24.7953758, -28.5864811, $ 
+             -32.3959961, -36.2186012, -40.0500298, -43.8870316, -47.7272758, -51.5692215, -55.4119453, -59.2549667, $
+             -63.0980797, -66.9412079, -70.7843399, -74.6274719, -78.4706039, -82.313736, -86.156868, -90.]
+endif
+
+; Modif ThL: for ECRAD, we have 14 SW bands and 16 LW bands.
+NSW=14
+NLW=16
+;
+;--script only works on ciclad
+dir='/bdd/input4MIPs/VOLC/${vv}/'
+;;;not available for v4
+;;;filename=dir+'CMIP_1850_2014_extinction_550nm_${vv}.nc'
+;;;NETCDFREAD,filename,'altitude',altitude,dimaltitude
+;;;NETCDFREAD,filename,'ext550',ext550,dimext
+;
+filename=dir+'CMIP_ECRAD_radiation_v4_1850-2018.nc'
+;
+NETCDFREAD,filename,'altitude',altitude,dimaltitude
+NETCDFREAD,filename,'latitude',latitude,dimlatitude
+NETCDFREAD,filename,'wl1_sun',wl1_sun,solar_bands
+NETCDFREAD,filename,'wl2_sun',wl2_sun,solar_bands
+NETCDFREAD,filename,'wl1_earth',wl1_earth,terrestrial_bands
+NETCDFREAD,filename,'wl2_earth',wl2_earth,terrestrial_bands
+NETCDFREAD,filename,'month',month,dimmonth
+NETCDFREAD,filename,'ext_sun',ext_sun,dimext
+NETCDFREAD,filename,'omega_sun',omega_sun,dimext
+NETCDFREAD,filename,'g_sun',g_sun,dimext
+NETCDFREAD,filename,'ext_earth',ext_earth,dimext
+NETCDFREAD,filename,'omega_earth',omega_earth,dimext
+NETCDFREAD,filename,'g_earth',g_earth,dimext
+;
+print ,'wl1_sun=', wl1_sun
+print ,'wl2_sun=', wl2_sun
+print ,'wl1_ear=', wl1_earth
+print ,'wl2_ear=', wl2_earth
+;
+print ,'min max ext_sun=', min(ext_sun),   max(ext_sun)
+print ,'min max ome_sun=', min(omega_sun), max(omega_sun)
+;for nl=0, NSW-1 do begin
+;print ,'min max ome_sun nl=', nl, min(omega_sun(*,*,*,nl)), max(omega_sun(*,*,*,nl))
+;endfor
+print ,'min max ggg_sun=', min(g_sun),     max(g_sun)
+print ,'min max ext_ear=', min(ext_earth), max(ext_earth)
+;
+dimlat=dimlatitude(0)
+dimalt=dimaltitude(0)
+dimtime=dimmonth(0)
+;
+;--vertical resolution in Beiping's code
+dz=0.5
+;
+month_in_year=12
+;;month_in_year=1
+;
+;--compute optical depth of input data
+;--dimension ext_sun(l,k,j,nl)
+;ext_sun0=ext_sun(0,*,*,*)
+;print , 'size ext_sun0=', size(ext_sun0)
+;tau_sun=TOTAL(ext_sun0,2)*dz
+;print ,'min max tau_sun l=0 =', min(tau_sun),   max(tau_sun)
+;
+zz=fltarr(dimz)
+;
+;---exact altitudes of LMDZ -- L79
+;--only exists for VLR, LR, MR and zoom_128x89 at the moment
+;
+filename='/bdd/input4MIPs/ZALT/zalt_zonmean_${lmdz}_rev.nc'
+NETCDFREAD,filename,'GEOP',zz,dimzz
+NETCDFREAD,filename,'LAT',zzlat,dimzzlat
+NETCDFREAD,filename,'TIME_COUNTER',zztime,dimzztime
+dimzzilat=dimzzlat[0]
+dimzzitime=dimzztime[0]
+;--becareful zz comes with four dimensions
+;--lon lat k time 
+print, 'GEOP size=', size(zz)
+;--becareful zzlat comes with South Pole first
+print, 'LAT from zalt field=', zzlat
+if (dimzzilat ne dimlatlmdz) then begin
+  print , 'PB dimension latitude'
+endif
+;
+;
+;--reconstructing the vertical coordinate at interfaces (in unit km)
+;--and reverse lat axis for zzi  i==>ii
+;--and reverse zz axis for zzi   dimz-1-k==>k
+;--and forget about lon axis : index 0 in zz
+zzi=fltarr(dimzzilat,dimzzitime,dimz+1)
+for i=0, dimzzilat-1 do begin
+ ii=dimzzilat-1-i
+ for t=0, dimzzitime-1 do begin
+  zzi(ii,t,0)=0.0
+  for k=1, dimz-1 do begin
+   zzi(ii,t,k)=(zz(0,i,dimz-1-(k-1),t)+zz(0,i,dimz-1-k,t))/2.0
+  endfor
+  zzi(ii,t,dimz)=100.00
+ endfor
+endfor
+print, 'zzi=',zzi
+;
+lev=indgen(dimz)+1
+;
+dimzori=dimalt
+zziori=fltarr(dimzori+1)
+zziori(0)=altitude(0)-0.25
+for k=1, dimzori do begin
+ zziori(k)=altitude(k-1)+0.25
+endfor
+;
+;--550 nm properties
+;;;not available for v4
+;;;tau_550_lmdz=fltarr(dimlatlmdz,dimz,month_in_year)
+;
+;--SW properties
+tau_sun_lmdz=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
+ome_sun_lmdz=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
+ggg_sun_lmdz=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
+;
+tau_sun_lmdz_ave=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
+ome_sun_lmdz_ave=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
+ggg_sun_lmdz_ave=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
+;
+tau_sun_lmdz_tr=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
+ome_sun_lmdz_tr=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
+ggg_sun_lmdz_tr=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
+;
+;--LW properties
+tau_ear_lmdz=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
+ome_ear_lmdz=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
+ggg_ear_lmdz=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
+;
+tau_ear_lmdz_ave=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
+ome_ear_lmdz_ave=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
+ggg_ear_lmdz_ave=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
+;
+tau_ear_lmdz_tr=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
+ome_ear_lmdz_tr=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
+ggg_ear_lmdz_tr=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
+;
+tau_sun_lmdz_ave(*,*,*,0)=1.e-15
+ome_sun_lmdz_ave(*,*,*,0)=1.e-15
+ggg_sun_lmdz_ave(*,*,*,0)=0.0
+tau_ear_lmdz_ave(*,*,*,0)=1.e-15
+ome_ear_lmdz_ave(*,*,*,0)=1.e-15
+ggg_ear_lmdz_ave(*,*,*,0)=0.0
+;
+yearend=dimtime/month_in_year-1
+;
+for year=0, yearend do begin
+;
+;;;not available for v4
+;;;tau_550_lmdz(*,*,*)=1.e-15
+tau_sun_lmdz(*,*,*,*)=1.e-15
+ome_sun_lmdz(*,*,*,*)=1.e-15
+ggg_sun_lmdz(*,*,*,*)=0.0
+tau_ear_lmdz(*,*,*,*)=1.e-15
+ome_ear_lmdz(*,*,*,*)=1.e-15
+ggg_ear_lmdz(*,*,*,*)=0.0
+;
+chyr=strcompress(1850+year,/rem)
+;
+for mth=0,month_in_year-1 do begin
+;
+;--timestep
+l=mth+month_in_year*year
+print,'year mth l=',chyr, mth, l
+;
+;regridding
+for j=0, dimlatlmdz-1  do begin
+;
+;--finding latitude in beiping luo data
+jj=0
+for jluo=0,dimlat-2 do begin
+  if (latitudelmdz(j) gt (latitude(jluo)+latitude(jluo+1))/2. ) then begin
+    jj=jluo+1
+  endif
+endfor
+;
+for k=0, dimz-1 do begin
+for kori=0, dimzori-1 do begin
+;
+;fraction de la maille kori qui se trouve dans la maille k
+frac= max([0.0,min([zzi(j,mth,k+1),zziori(kori+1)])-max([zzi(j,mth,k),zziori(kori)])])/(zziori(kori+1)-zziori(kori))
+;
+;;;not available for v4
+;;;tau_550_lmdz(j,k,mth)=tau_550_lmdz(j,k,mth)+ext550(l,kori,jj)*dz*frac
+;
+for nl=0, NSW-1 do begin
+  tau_sun_lmdz(j,k,nl,mth)=tau_sun_lmdz(j,k,nl,mth)+ext_sun(l,kori,jj,nl)*dz*frac
+  ome_sun_lmdz(j,k,nl,mth)=ome_sun_lmdz(j,k,nl,mth)+ext_sun(l,kori,jj,nl)*omega_sun(l,kori,jj,nl)*dz*frac
+  ggg_sun_lmdz(j,k,nl,mth)=ggg_sun_lmdz(j,k,nl,mth)+ext_sun(l,kori,jj,nl)*omega_sun(l,kori,jj,nl)*g_sun(l,kori,jj,nl)*dz*frac
+;
+  tau_sun_lmdz_ave(j,k,nl,0)=tau_sun_lmdz_ave(j,k,nl,0)+ext_sun(l,kori,jj,nl)*dz*frac
+  ome_sun_lmdz_ave(j,k,nl,0)=ome_sun_lmdz_ave(j,k,nl,0)+ext_sun(l,kori,jj,nl)*omega_sun(l,kori,jj,nl)*dz*frac
+  ggg_sun_lmdz_ave(j,k,nl,0)=ggg_sun_lmdz_ave(j,k,nl,0)+ext_sun(l,kori,jj,nl)*omega_sun(l,kori,jj,nl)*g_sun(l,kori,jj,nl)*dz*frac
+endfor
+;
+; For LW, script differs from previous versions in that
+; (1.-omega_earth) factor becomes (omega_earth)
+;
+for nl=0, NLW-1 do begin
+  tau_ear_lmdz(j,k,nl,mth)=tau_ear_lmdz(j,k,nl,mth)+ext_earth(l,kori,jj,nl)*dz*frac
+  ome_ear_lmdz(j,k,nl,mth)=ome_ear_lmdz(j,k,nl,mth)+ext_earth(l,kori,jj,nl)*omega_earth(l,kori,jj,nl)*dz*frac
+  ggg_ear_lmdz(j,k,nl,mth)=ggg_ear_lmdz(j,k,nl,mth)+ext_earth(l,kori,jj,nl)*omega_earth(l,kori,jj,nl)*g_earth(l,kori,jj,nl)*dz*frac
+;
+  tau_ear_lmdz_ave(j,k,nl,0)=tau_ear_lmdz_ave(j,k,nl,0)+ext_earth(l,kori,jj,nl)*dz*frac
+  ome_ear_lmdz_ave(j,k,nl,0)=ome_ear_lmdz_ave(j,k,nl,0)+ext_earth(l,kori,jj,nl)*omega_earth(l,kori,jj,nl)*dz*frac
+  ggg_ear_lmdz_ave(j,k,nl,0)=ggg_ear_lmdz_ave(j,k,nl,0)+ext_earth(l,kori,jj,nl)*omega_earth(l,kori,jj,nl)*g_earth(l,kori,jj,nl)*dz*frac
+endfor
+;
+endfor
+;--end lat loop
+;
+endfor
+endfor 
+;--end k loops
+;
+endfor
+;--end month loop
+;
+;renormalizing intensive SW properties
+ggg_sun_lmdz(*,*,*,*)=ggg_sun_lmdz(*,*,*,*)/ome_sun_lmdz(*,*,*,*)
+ome_sun_lmdz(*,*,*,*)=ome_sun_lmdz(*,*,*,*)/tau_sun_lmdz(*,*,*,*)
+;
+;saving netcdf file
+;
+print ,'min max ext_sun lmdz=', min(tau_sun_lmdz), max(tau_sun_lmdz)
+print ,'min max ome_sun lmdz=', min(ome_sun_lmdz), max(ome_sun_lmdz)
+print ,'min max ggg_sun lmdz=', min(ggg_sun_lmdz), max(ggg_sun_lmdz)
+print ,'min max ext_ear lmdz=', min(tau_ear_lmdz), max(tau_ear_lmdz)
+print ,'min max ome_ear lmdz=', min(ome_ear_lmdz), max(ome_ear_lmdz)
+print ,'min max ggg_ear lmdz=', min(ggg_ear_lmdz), max(ggg_ear_lmdz)
+;
+print ,'min max ext_sun lmdz l=1=', min(tau_sun_lmdz(*,*,*,0)), max(tau_sun_lmdz(*,*,*,0))
+print ,'min max ome_sun lmdz l=1=', min(ome_sun_lmdz(*,*,*,0)), max(ome_sun_lmdz(*,*,*,0))
+print ,'min max ggg_sun lmdz l=1=', min(ggg_sun_lmdz(*,*,*,0)), max(ggg_sun_lmdz(*,*,*,0))
+print ,'min max ext_ear lmdz l=1=', min(tau_ear_lmdz(*,*,*,0)), max(tau_ear_lmdz(*,*,*,0))
+print ,'min max ome_ear lmdz l=1=', min(ome_ear_lmdz(*,*,*,0)), max(ome_ear_lmdz(*,*,*,0))
+print ,'min max ggg_ear lmdz l=1=', min(ggg_ear_lmdz(*,*,*,0)), max(ggg_ear_lmdz(*,*,*,0))
+
+;
+;--compute optical depth of output data
+;tau_sun=TOTAL(tau_sun_lmdz,2)
+;print ,'min max tau_sun_lmdz vert=', min(tau_sun),   max(tau_sun)
+;
+opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
+             wav:fltarr(NSW),time:fltarr(month_in_year),          $ 
+             tau_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
+             ome_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
+             ggg_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year) }
+;
+opticstruct.lat=latitudelmdz
+opticstruct.lev=lev
+opticstruct.wav=(wl1_sun(0:NSW-1)+wl2_sun(0:NSW-1))/2.
+opticstruct.time=float(indgen(month_in_year)+1)
+opticstruct.tau_sun=tau_sun_lmdz
+opticstruct.ome_sun=ome_sun_lmdz
+opticstruct.ggg_sun=ggg_sun_lmdz
+;
+attributes = {units:strarr(7),long_name:strarr(7)}
+attributes.units =     ['degrees_north','level','Âµm','month','-','-','-']
+attributes.long_name = ['latitude','level','wavelength','time','tau_sun','ome_sun','g_sun']
+;
+dimensions = {isdim:intarr(7), links:intarr(4,7)}
+       dimensions.isdim =  [1,1,1,1,0,0,0]  ; (1=dimension, 0=variable)
+       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
+                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
+                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
+;
+netcdfwrite,output+'tauswstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
+            attributes=attributes, dimensions=dimensions
+;
+; NOW FOR LW
+;
+opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
+             wav:fltarr(NLW),time:fltarr(month_in_year),          $
+             tau_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year),   $
+             ome_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year),   $
+             ggg_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year) }              
+;
+opticstruct.lat=latitudelmdz
+opticstruct.lev=lev
+opticstruct.wav=(wl1_earth+wl2_earth)/2.
+opticstruct.time=float(indgen(month_in_year)+1)
+opticstruct.tau_ear=tau_ear_lmdz
+opticstruct.ome_ear=ome_sun_lmdz
+opticstruct.ggg_ear=ggg_sun_lmdz
+
+;
+attributes = {units:strarr(7),long_name:strarr(7)}
+attributes.units = ['degrees_north','level','Âµm','month','-','-','-']
+attributes.long_name = ['latitude','level','wavelength','time','tau_ear','ome_ear','g_ear']
+;
+dimensions = {isdim:intarr(7), links:intarr(4,7)}
+       dimensions.isdim =  [1,1,1,1,0,0,0]  ; (1=dimension, 0=variable)
+       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
+                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
+                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
+;
+netcdfwrite,output+'taulwstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
+            attributes=attributes, dimensions=dimensions
+;
+;;;not available in v4
+;;;opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
+;;;             time:fltarr(month_in_year),                          $
+;;;             tau550:fltarr(dimlatlmdz,dimz,month_in_year) }
+;
+;;;opticstruct.lat=latitudelmdz
+;;;opticstruct.lev=lev
+;;;opticstruct.time=float(indgen(month_in_year)+1)
+;;;opticstruct.tau550=tau_550_lmdz
+;
+;;;attributes = {units:strarr(4),long_name:strarr(4)}
+;;;attributes.units = ['degrees_north','level','month','-']
+;;;attributes.long_name = ['latitude','level','time','tau550']
+;
+;;;dimensions = {isdim:intarr(4), links:intarr(3,4)}
+;;;       dimensions.isdim =  [1,1,1,0]  ; (1=dimension, 0=variable)
+;;;       dimensions.links = [ [-1,-1,-1],[-1,-1,-1],   $
+;;;                            [-1,-1,-1],[0,1,2]  ]
+;
+;;;netcdfwrite,output+'tau550strat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
+;;;            attributes=attributes, dimensions=dimensions
+;
+endfor
+;--end loop on years
+;
+; Now deal with average conditions
+;
+ggg_sun_lmdz_ave(*,*,*,0)=ggg_sun_lmdz_ave(*,*,*,0)/ome_sun_lmdz_ave(*,*,*,0)
+ome_sun_lmdz_ave(*,*,*,0)=ome_sun_lmdz_ave(*,*,*,0)/tau_sun_lmdz_ave(*,*,*,0)
+tau_sun_lmdz_ave(*,*,*,0)=tau_sun_lmdz_ave(*,*,*,0)/float(dimtime)
+;
+ggg_ear_lmdz_ave(*,*,*,0)=ggg_ear_lmdz_ave(*,*,*,0)/ome_ear_lmdz_ave(*,*,*,0)
+ome_ear_lmdz_ave(*,*,*,0)=ome_ear_lmdz_ave(*,*,*,0)/tau_ear_lmdz_ave(*,*,*,0)
+tau_ear_lmdz_ave(*,*,*,0)=tau_ear_lmdz_ave(*,*,*,0)/float(dimtime)
+;
+for mth=1, month_in_year-1 do begin
+  ggg_sun_lmdz_ave(*,*,*,mth)=ggg_sun_lmdz_ave(*,*,*,0)
+  ome_sun_lmdz_ave(*,*,*,mth)=ome_sun_lmdz_ave(*,*,*,0)
+  tau_sun_lmdz_ave(*,*,*,mth)=tau_sun_lmdz_ave(*,*,*,0)
+  ggg_ear_lmdz_ave(*,*,*,mth)=ggg_ear_lmdz_ave(*,*,*,0)
+  ome_ear_lmdz_ave(*,*,*,mth)=ome_ear_lmdz_ave(*,*,*,0)
+  tau_ear_lmdz_ave(*,*,*,mth)=tau_ear_lmdz_ave(*,*,*,0)
+endfor
+;
+opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
+             wav:fltarr(NSW),time:fltarr(month_in_year),          $ 
+             tau_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
+             ome_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
+             ggg_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year) }
+;
+opticstruct.lat=latitudelmdz
+opticstruct.lev=lev
+opticstruct.wav=(wl1_sun(0:NSW-1)+wl2_sun(0:NSW-1))/2.
+opticstruct.time=float(indgen(month_in_year))
+opticstruct.tau_sun=tau_sun_lmdz_ave
+opticstruct.ome_sun=ome_sun_lmdz_ave
+opticstruct.ggg_sun=ggg_sun_lmdz_ave
+;
+attributes = {units:strarr(7),long_name:strarr(7)}
+attributes.units =     ['degrees_north','level','Âµm','month','-','-','-']
+attributes.long_name = ['latitude','level','wavelength','time','tau_sun','ome_sun','g_sun']
+;
+dimensions = {isdim:intarr(7), links:intarr(4,7)}
+       dimensions.isdim =  [1,1,1,1,0,0,0]  ; (1=dimension, 0=variable)
+       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
+                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
+                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
+;
+netcdfwrite,output+'tauswstrat.2D.ave.nc',opticstruct,clobber=1, $
+            attributes=attributes, dimensions=dimensions
+;
+opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
+             wav:fltarr(NLW),time:fltarr(month_in_year),          $
+             tau_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year),   $
+             ome_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year),   $
+             ggg_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year) }
+;
+opticstruct.lat=latitudelmdz
+opticstruct.lev=lev
+opticstruct.wav=(wl1_earth+wl2_earth)/2.
+opticstruct.time=float(indgen(month_in_year))
+opticstruct.tau_ear=tau_ear_lmdz_ave
+opticstruct.ome_ear=tau_ome_lmdz_ave
+opticstruct.ggg_ear=tau_ggg_lmdz_ave
+;
+attributes = {units:strarr(7),long_name:strarr(7)}
+attributes.units = ['degrees_north','level','Âµm','month','-','-','-']
+attributes.long_name = ['latitude','level','wavelength','time','tau_ear','ome_ear','g_ear']
+;
+dimensions = {isdim:intarr(7), links:intarr(4,7)}
+       dimensions.isdim =  [1,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,-1],   $
+                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
+;
+netcdfwrite,output+'taulwstrat.2D.ave.nc',opticstruct,clobber=1, $
+            attributes=attributes, dimensions=dimensions
+;
+;--finally prepare the 10-year transition period 
+;--that is 2015-2023 for v3 or 2017-2025 for v4
+;--we mix the last year of the period with the climatology
+;
+for year=1850+yearend+1,1850+yearend+9 do begin
+;
+chyr=strcompress(year,/rem)
+print,'year =',chyr
+;--compute weights
+w1=float(1850+yearend+10-year)/10.
+w2=1.0-w1
+;--SW properties
+tau_sun_lmdz_tr=w1*tau_sun_lmdz+w2*tau_sun_lmdz_ave
+ome_sun_lmdz_tr=w1*tau_sun_lmdz*ome_sun_lmdz+w2*tau_sun_lmdz_ave*ome_sun_lmdz_ave
+ggg_sun_lmdz_tr=w1*tau_sun_lmdz*ome_sun_lmdz*ggg_sun_lmdz+ $ 
+                w2*tau_sun_lmdz_ave*ome_sun_lmdz_ave*ggg_sun_lmdz_ave
+ggg_sun_lmdz_tr=ggg_sun_lmdz_tr/ome_sun_lmdz_tr
+ome_sun_lmdz_tr=ome_sun_lmdz_tr/tau_sun_lmdz_tr
+;
+;--LW properties
+tau_ear_lmdz_tr=w1*tau_ear_lmdz+w2*tau_ear_lmdz_ave
+ome_ear_lmdz_tr=w1*tau_ear_lmdz*ome_ear_lmdz+w2*tau_ear_lmdz_ave*ome_ear_lmdz_ave
+ggg_ear_lmdz_tr=w1*tau_ear_lmdz*ome_ear_lmdz*ggg_ear_lmdz+ $ 
+                w2*tau_ear_lmdz_ave*ome_ear_lmdz_ave*ggg_ear_lmdz_ave
+ggg_ear_lmdz_tr=ggg_ear_lmdz_tr/ome_ear_lmdz_tr
+ome_ear_lmdz_tr=ome_ear_lmdz_tr/tau_ear_lmdz_tr
+;
+;--prepare SW output
+opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
+             wav:fltarr(NSW),time:fltarr(month_in_year),          $ 
+             tau_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
+             ome_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
+             ggg_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year) }
+;
+opticstruct.lat=latitudelmdz
+opticstruct.lev=lev
+opticstruct.wav=(wl1_sun(0:NSW-1)+wl2_sun(0:NSW-1))/2.
+opticstruct.time=float(indgen(month_in_year)+1)
+opticstruct.tau_sun=tau_sun_lmdz_tr
+opticstruct.ome_sun=ome_sun_lmdz_tr
+opticstruct.ggg_sun=ggg_sun_lmdz_tr
+;
+attributes = {units:strarr(7),long_name:strarr(7)}
+attributes.units =     ['degrees_north','level','Âµm','month','-','-','-']
+attributes.long_name = ['latitude','level','wavelength','time','tau_sun','ome_sun','g_sun']
+;
+dimensions = {isdim:intarr(7), links:intarr(4,7)}
+       dimensions.isdim =  [1,1,1,1,0,0,0]  ; (1=dimension, 0=variable)
+       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
+                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
+                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
+;
+netcdfwrite,output+'tauswstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
+            attributes=attributes, dimensions=dimensions
+;
+;--prepare LW output
+opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
+             wav:fltarr(NLW),time:fltarr(month_in_year),          $
+             tau_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year),   $
+             ome_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year),   $
+             ggg_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year) }
+;
+opticstruct.lat=latitudelmdz
+opticstruct.lev=lev
+opticstruct.wav=(wl1_earth+wl2_earth)/2.
+opticstruct.time=float(indgen(month_in_year)+1)
+opticstruct.tau_ear=tau_ear_lmdz_tr
+opticstruct.ome_ear=ome_ear_lmdz_tr
+opticstruct.ggg_ear=ggg_ear_lmdz_tr
+;
+attributes = {units:strarr(7),long_name:strarr(7)}
+attributes.units = ['degrees_north','level','Âµm','month','-','-','-']
+attributes.long_name = ['latitude','level','wavelength','time','tau_ear','ome_ear','g_ear']
+;
+dimensions = {isdim:intarr(7), links:intarr(4,7)}
+       dimensions.isdim =  [1,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,-1],   $
+                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
+;
+netcdfwrite,output+'taulwstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
+            attributes=attributes, dimensions=dimensions
+;
+endfor
+;--end loop on years
+;
+end
+EOF
+
+cat > volc_${lmdz}.job << EOF2
+#PBS -N process_volc_${lmdz}
+#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
+
+cd $dirpwd
+
+idl << eof
+.r netcdf 
+.r process_volc_${lmdz}
+regrid
+eof
+EOF2
+
+qsub volc_${lmdz}.job
