Index: TOOLS/CMIP6_FORCING/SOLAR/ECRAD/CMIP6_14bands_ECRAD_spectralsolarscaling.m
===================================================================
--- TOOLS/CMIP6_FORCING/SOLAR/ECRAD/CMIP6_14bands_ECRAD_spectralsolarscaling.m	(revision 6266)
+++ TOOLS/CMIP6_FORCING/SOLAR/ECRAD/CMIP6_14bands_ECRAD_spectralsolarscaling.m	(revision 6266)
@@ -0,0 +1,322 @@
+% Updated 05/02/2018, ThL
+% * changed paths to retrieve new version of input data;
+% * version considered is now 3.2.
+
+% Updated 03/05/2022, ThL
+% * now adapted to deal with 14 bands of ECRAD. 
+
+% Updated 22/09/2022, ThL
+% * added computation of spectral_solar_scaling vector for ECRAD (over 14 wavebands).
+
+
+% MATLAB run script to integrate spectral solar irradiance (SSI) data
+% recommended and provided by SPARC/SOLARIS-HEPPA for use in CMIP6 to
+% spectral bands of any desired resolution (according to
+% radiation/photolysis scheme of employed (chemistry) climate model)
+%
+% Author: Tim Kruschke (tkruschke[at]geomar.de)
+% Last modified: 02 May 2016 (compliance with solarforcing_*_*_3.1.nc)
+%
+% Description:
+%    1. reading data from netcdf-file as provided by SPARC/SOLARIS-HEPPA
+%    2. reading definitions of spectral bands (user-specified by creating an
+%       ascii-textfile according to specific formatting, as described below)
+%    3. calling MATLAB-function "regrid_ssi" (provided with this run script as
+%       part of a joint zip-archive) to integrate SSI over specified bands
+%    4. writing new netcdf-file with same structure and attributes as original
+%       (input) file, but only with (unchanged) total solar irradiance
+%       (TSI), SSI integrated over spectral bands specified by user, and
+%       irradiance fraction (of TSI) in respective band (plus
+%       necessary/related dimensions/variables time, wlen, wlenbinsize)
+%
+% Formatting of ascii-textfile, containing definitions of spectral bands:
+%    NO HEADER LINES, NO UNITS, NO CHARACTER DELIMITERS
+%    It's as simple as this:
+%    - one line per spectral band
+%    - each line contains two floating point numbers: lower and upper edge
+%      of respective spectral band (in nanometers!!!), separated by
+%      blank space(s) or tabulator(s)
+%    See example_spectral_bands.dat (provided with this run script as
+%    part of a joint zip-archive) as an example.
+%--------------------------------------------------------------------------
+
+%  execute with run CMIP6.m
+
+%% User-defined paths and filenames
+
+% path to directory of input data; adapt to your file structure
+dir_in=['.'];
+
+% path to directory of output data; adapt to your file structure if desired
+% to differ from input directory
+dir_out=dir_in;
+
+% filename of ascii-textfile (to be placed in dir_in) containing
+% definitions of spectral bands (see explanations in script header);
+% adapt to name of your band definition file
+%filename_spectral_bands='example_spectral_bands.dat';
+filename_spectral_bands='ECRAD_14bands_nm.dat';
+%filename_spectral_bands='example_spectral_bands_2_SSI_FRAC=1.dat';
+
+% filename of original data file, provided by SPARC/SOLARIS-HEPPA (to be
+% placed in dir_in); adapt to respective file that is to be considered
+%filename_in='solarforcing_ref_mon_3.1.nc';
+filename_in='/bdd/input4MIPs/CMIP6/CMIP/SOLARIS-HEPPA/SOLARIS-HEPPA-3-2/atmos/day/multiple/gz/v20170103/solarforcing-ref-day_input4MIPs_solar_CMIP_SOLARIS-HEPPA-3-2_gn_18500101-22991231.nc'
+
+% lines 53-54 only used to match input and output filenames
+% ATTENTION: works only on unix platforms; if run on other platforms,
+% delete/comment following two lines and enter complete filename_out in
+% l. 59
+[~,filename_in_base]=unix(['basename ' dir_in filesep filename_in ' .nc']);
+filename_in_base=strcat(filename_in_base);
+
+% filename of new file to be created, containing SSI data in spectral
+% resolution as desired by the user; adapt to whatever name you like if
+% not desired to match input filename
+filename_out=[filename_in_base '_ECRAD_14bands.nc'];
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+%% Usually, nothing should be edited below this line!!!!!!!!!!!!!!!!!!!!!!!
+
+ncid_in=netcdf.open([dir_in filesep filename_in],'NOWRITE');
+ncid_out=netcdf.create(strcat([dir_out filesep filename_out]),'CLOBBER'); % existing files are overwritten
+
+% adapting title of original file
+nc_title=netcdf.getAtt(ncid_in,-1,'title');
+nc_title=['Spectral solar irradiance (SSI) of ' nc_title ', converted to custom spectral resolution'];
+netcdf.putAtt(ncid_out,-1,'title',nc_title);
+
+%% time stamp
+netcdf.putAtt(ncid_out,-1,'creation_date',datestr(now)); % creates a time stamp, though not compliant with ISO 8601
+%c% You can use the following lines (commented out) instead to generate a
+%c% time stamp compliant to ISO 8601; change last digits according to
+%c% your time zone in this case (currently +02:00 for CEST)
+%     create_date_vec=datevec(now);
+%     create_date_str=cell(1,6);
+%     for id=1:6;
+%         if create_date_vec(id)<10;
+%             create_date_str{1,id}=['0' num2str(round(create_date_vec(id)))];
+%         else
+%             create_date_str{1,id}=num2str(round(create_date_vec(id)));
+%         end
+%     end
+%     netcdf.putAtt(ncid_out,-1,'creation_date',[create_date_str{1,1} '-' create_date_str{1,2} '-' ...
+%         create_date_str{1,3} 'T' create_date_str{1,4} ':' create_date_str{1,5} ':' create_date_str{1,6} '+02:00'])
+%%
+
+% copy global attributes to new file
+netcdf.copyAtt(ncid_in,-1,'institution',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'institution_id',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'activity_id',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'comment',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'time_coverage_start',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'time_coverage_end',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'frequency',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'source',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'source_id',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'realm',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'further_info_url',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'metadata_url',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'contributor_name',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'references',ncid_out,-1)
+netcdf.copyAtt(ncid_in,-1,'Conventions',ncid_out,-1)
+netcdf.putAtt(ncid_out,-1,'history',['Original SSI data from file ' filename_in ' integrated ' ...
+    'over spectral bands as defined in file ' filename_spectral_bands ' ' ...
+    'by MATLAB run script CMIP6_SSI_convert_resolution_run.m, using function regrid_ssi.m']);
+
+% reading SSI-related variables from CMIP6 input file
+wvl_in_ID=netcdf.inqVarID(ncid_in,'wlen');
+wvl_edges_in_ID=netcdf.inqVarID(ncid_in,'wlen_bnds');
+bw_in_ID=netcdf.inqVarID(ncid_in,'wlenbinsize');
+time_in_ID=netcdf.inqVarID(ncid_in,'time');
+time_bnds_in_ID=netcdf.inqVarID(ncid_in,'time_bnds');
+year_in_ID=netcdf.inqVarID(ncid_in,'calyear');
+month_in_ID=netcdf.inqVarID(ncid_in,'calmonth');
+day_in_ID=netcdf.inqVarID(ncid_in,'calday');
+tsi_in_ID=netcdf.inqVarID(ncid_in,'tsi');
+ssi_in_ID=netcdf.inqVarID(ncid_in,'ssi');
+
+wvl_in=netcdf.getVar(ncid_in,wvl_in_ID);
+wvl_edges_in=(netcdf.getVar(ncid_in,wvl_edges_in_ID))';
+bw_in=netcdf.getVar(ncid_in,bw_in_ID);
+ssi_in=netcdf.getVar(ncid_in,ssi_in_ID);
+tsi_in=netcdf.getVar(ncid_in,tsi_in_ID);
+time_in=netcdf.getVar(ncid_in,time_in_ID);
+time_bnds_in=netcdf.getVar(ncid_in,time_bnds_in_ID);
+year_in=netcdf.getVar(ncid_in,year_in_ID);
+month_in=netcdf.getVar(ncid_in,month_in_ID);
+day_in=netcdf.getVar(ncid_in,day_in_ID);
+
+% read user specified definitions of spectral bands
+fid=fopen([dir_in filesep filename_spectral_bands]);
+tmp=textscan(fid,'%f %f');
+fclose(fid);
+wvl_edges_out=tmp{1,1};
+wvl_edges_out(:,2)=tmp{1,2};
+clear tmp
+
+% calculation of output central wave lengths
+wvl_out=sum(wvl_edges_out,2)/2;
+
+% calculation of output bin width
+bw_out=wvl_edges_out(:,2)-wvl_edges_out(:,1);
+
+% check for consistency of output bin definitions
+if min(bw_out)<0;
+    error(['Error: For at least one spectral band, the specified upper edge \n' ...
+        '(right column) in your input file ' filename_in '\n ' ...
+        'is smaller than the respective lower edge (left column).'])
+end
+
+% check for spectral coverage of input data with respect to desired output
+% bins
+if min(wvl_edges_out(:,1))<min(wvl_edges_in(:,1)) || max(wvl_edges_out(:,2))>max(wvl_edges_in(:,2));
+    warning(['You are extrapolating into spectral ranges not covered by the input spectrum! ' ...
+        'This routine supports this by linearly extrapolating based on change between last two bins. ' ...
+        'However, this might be physically inappropriate!']);
+end
+
+% Integration of SSI data for each specified spectral band
+% (regrid_ssi.m calculates weighted means of SSI per nm, this is multiplied by binwidth)
+ssi_out=bsxfun(@times,regrid_ssi(ssi_in,wvl_edges_in,wvl_edges_out),bw_out);
+
+% Calculation of solar irradiance fraction (compared to total solar
+% irradiance) in each band
+ssi_frac_out=bsxfun(@rdivide,ssi_out,tsi_in');
+
+% Calculation the average value of tsi
+avtsi=mean(tsi_in)
+disp (avtsi)
+
+%Calculate the anomaly of tsi
+antsi(:,1)=tsi_in(:,1)-avtsi
+
+disp (antsi)
+
+format longEng
+disp (avtsi)
+
+% Now dealing with creation of the "spectral_solar_scaling" vector
+
+% Number of g-points per j-waveband (j in [1:14])
+q = [6,12,8,8,10,10,2,10,8,6,6,8,6,12];
+
+% Cumulative sum...
+qsum = cumsum(q);
+% ... used to yield indices of min and max g-points (in [1:112]) tallying with each j-band ([1:14])
+Nmax = qsum;
+Nmin = [1];
+for j=2:14, Nmin(j) = qsum(j-1)+1; end;
+
+% importing SSI_i^ECRAD tabulated values
+SSI_ECRAD=csvread('q_i.txt');
+
+specsolscal = zeros(14,length(time_in));
+
+if ( size(specsolscal) ~= size(ssi_out) )
+  disp('Mismatch in dimensions!');
+  return
+end
+
+for j=1:14,
+  % According to literal derivation
+  % r_j = SSI_j^CMIP6 / SSI_j^ECRAD
+  specsolscal(j,:) = ssi_out(j,:) / sum(SSI_ECRAD(Nmin(j):Nmax(j)));
+end
+
+% writing data (including dimensions) to output file
+  % dimension definitions
+time_dimID=netcdf.defDim(ncid_out,'time',0);
+wvl_dimID=netcdf.defDim(ncid_out,'wlen',length(wvl_out));
+bounds_dimID=netcdf.defDim(ncid_out,'nbd',2);
+
+  % variable definitions
+time_out_ID=netcdf.defVar(ncid_out,'time','double',time_dimID);
+time_bnds_out_ID=netcdf.defVar(ncid_out,'time_bnds','double',[bounds_dimID time_dimID]);
+wvl_out_ID=netcdf.defVar(ncid_out,'wlen','double',wvl_dimID);
+wvl_edges_out_ID=netcdf.defVar(ncid_out,'wlen_bnds','double',[wvl_dimID bounds_dimID]);
+bw_out_ID=netcdf.defVar(ncid_out,'wlenbinsize','double',wvl_dimID);
+year_out_ID=netcdf.defVar(ncid_out,'calyear','double',time_dimID);
+month_out_ID=netcdf.defVar(ncid_out,'calmonth','double',time_dimID);
+day_out_ID=netcdf.defVar(ncid_out,'calday','double',time_dimID);
+tsi_out_ID=netcdf.defVar(ncid_out,'tsi','float',time_dimID);
+antsi_out_ID=netcdf.defVar(ncid_out,'antsi','float',time_dimID);
+ssi_out_ID=netcdf.defVar(ncid_out,'ssi','float',[wvl_dimID time_dimID]);
+ssi_frac_out_ID=netcdf.defVar(ncid_out,'ssi_frac','float',[wvl_dimID time_dimID]);
+spectral_solar_scaling_out_ID=netcdf.defVar(ncid_out,'spectral_solar_scaling','float',[wvl_dimID time_dimID]);
+
+  % copy variable attributes from input to output file
+netcdf.copyAtt(ncid_in,time_in_ID,'axis',ncid_out,time_out_ID)
+netcdf.copyAtt(ncid_in,time_in_ID,'long_name',ncid_out,time_out_ID)
+netcdf.copyAtt(ncid_in,time_in_ID,'standard_name',ncid_out,time_out_ID)
+netcdf.copyAtt(ncid_in,time_in_ID,'units',ncid_out,time_out_ID)
+netcdf.copyAtt(ncid_in,time_in_ID,'calendar',ncid_out,time_out_ID)
+netcdf.copyAtt(ncid_in,time_in_ID,'_CoordinateAxisType',ncid_out,time_out_ID)
+netcdf.copyAtt(ncid_in,time_in_ID,'bounds',ncid_out,time_out_ID)
+
+netcdf.copyAtt(ncid_in,time_bnds_in_ID,'long_name',ncid_out,time_bnds_out_ID)
+netcdf.copyAtt(ncid_in,time_bnds_in_ID,'units',ncid_out,time_bnds_out_ID)
+
+netcdf.copyAtt(ncid_in,wvl_in_ID,'standard_name',ncid_out,wvl_out_ID)
+netcdf.copyAtt(ncid_in,wvl_in_ID,'long_name',ncid_out,wvl_out_ID)
+netcdf.copyAtt(ncid_in,wvl_in_ID,'units',ncid_out,wvl_out_ID)
+netcdf.copyAtt(ncid_in,wvl_in_ID,'bounds',ncid_out,wvl_out_ID)
+
+netcdf.copyAtt(ncid_out,wvl_edges_in_ID,'long_name',ncid_out,wvl_edges_out_ID)
+netcdf.copyAtt(ncid_out,wvl_edges_in_ID,'units',ncid_out,wvl_edges_out_ID)
+
+netcdf.copyAtt(ncid_in,bw_in_ID,'long_name',ncid_out,bw_out_ID)
+netcdf.copyAtt(ncid_in,bw_in_ID,'units',ncid_out,bw_out_ID)
+
+netcdf.copyAtt(ncid_in,year_in_ID,'long_name',ncid_out,year_out_ID)
+
+netcdf.copyAtt(ncid_in,month_in_ID,'long_name',ncid_out,month_out_ID)
+
+netcdf.copyAtt(ncid_in,day_in_ID,'long_name',ncid_out,day_out_ID)
+
+netcdf.copyAtt(ncid_in,tsi_in_ID,'standard_name',ncid_out,tsi_out_ID)
+netcdf.copyAtt(ncid_in,tsi_in_ID,'long_name',ncid_out,tsi_out_ID)
+netcdf.copyAtt(ncid_in,tsi_in_ID,'units',ncid_out,tsi_out_ID)
+netcdf.copyAtt(ncid_in,tsi_in_ID,'cell_methods',ncid_out,tsi_out_ID)
+
+netcdf.putAtt(ncid_out,antsi_out_ID,'standard_name','total solar_irradiance anomaly')
+netcdf.putAtt(ncid_out,antsi_out_ID,'long_name','reconstructed total solar irradiance anomaly at 1 AU')
+netcdf.putAtt(ncid_out,antsi_out_ID,'units','W m^-2')
+
+netcdf.putAtt(ncid_out,ssi_out_ID,'long_name','reconstructed solar irradiance at 1 AU within respective bands')
+netcdf.putAtt(ncid_out,ssi_out_ID,'units','W m^-2')
+
+netcdf.putAtt(ncid_out,ssi_frac_out_ID,'long_name','reconstructed fraction of total solar irradiance at 1 AU within respective bands')
+netcdf.putAtt(ncid_out,ssi_frac_out_ID,'units','1')
+
+netcdf.putAtt(ncid_out,spectral_solar_scaling_out_ID,'long_name','spectral solar scaling necessary for ECRAD code')
+netcdf.putAtt(ncid_out,spectral_solar_scaling_out_ID,'units','1')
+
+netcdf.close(ncid_in);
+netcdf.endDef(ncid_out);
+
+netcdf.putVar(ncid_out,time_out_ID,0,length(time_in),time_in)
+netcdf.putVar(ncid_out,time_bnds_out_ID,time_bnds_in)
+netcdf.putVar(ncid_out,wvl_out_ID,wvl_out)
+netcdf.putVar(ncid_out,wvl_edges_out_ID,wvl_edges_out)
+netcdf.putVar(ncid_out,bw_out_ID,bw_out)
+netcdf.putVar(ncid_out,year_out_ID,year_in)
+netcdf.putVar(ncid_out,month_out_ID,month_in)
+netcdf.putVar(ncid_out,day_out_ID,day_in)
+netcdf.putVar(ncid_out,tsi_out_ID,tsi_in)
+netcdf.putVar(ncid_out,antsi_out_ID,antsi)
+netcdf.putVar(ncid_out,ssi_out_ID,ssi_out)
+netcdf.putVar(ncid_out,ssi_frac_out_ID,ssi_frac_out)
+netcdf.putVar(ncid_out,spectral_solar_scaling_out_ID,specsolscal)
+
+netcdf.sync(ncid_out);
+netcdf.close(ncid_out)
+
Index: TOOLS/CMIP6_FORCING/SOLAR/ECRAD/ECRAD_14bands_nm.dat
===================================================================
--- TOOLS/CMIP6_FORCING/SOLAR/ECRAD/ECRAD_14bands_nm.dat	(revision 6266)
+++ TOOLS/CMIP6_FORCING/SOLAR/ECRAD/ECRAD_14bands_nm.dat	(revision 6266)
@@ -0,0 +1,17 @@
+3076.92307692308 3846.15384615385
+2500.0           3076.92307692308
+2150.5376344086	 2500.0
+1941.74757281553 2150.5376344086
+1626.0162601626	 1941.74757281553
+1298.7012987013	 1626.0162601626
+1242.23602484472 1298.7012987013
+778.210116731518 1242.23602484472
+625.0            778.210116731518
+441.501103752759 625.0
+344.827586206897 441.501103752759
+263.157894736842 344.827586206897
+200.0            263.157894736842
+3846.15384615385 12195.1219512195
+
+
+
Index: TOOLS/CMIP6_FORCING/SOLAR/ECRAD/README
===================================================================
--- TOOLS/CMIP6_FORCING/SOLAR/ECRAD/README	(revision 6266)
+++ TOOLS/CMIP6_FORCING/SOLAR/ECRAD/README	(revision 6266)
@@ -0,0 +1,14 @@
+=======================
+Script CMIP6_14bands_ECRAD_spectralsolarscaling.m was adapted to comply with the ECRAD configuration; it is based on the previous versions of the matlab script.
+This script currently takes as input the solarforcing-ref-day_input4MIPs_solar_CMIP_SOLARIS-HEPPA-3-2_gn_18500101-22991231.nc file as stored on the /bdd/,
+and generates a solarforcing-ref-day_input4MIPs_solar_CMIP_SOLARIS-HEPPA-3-2_gn_18500101-22991231_ECRAD_14bands.nc file over the 14 ECRAD wavebands.
+This output file is stored in the working directory.
+
+The computation relies on the ECRAD_14bands_nm.dat file.
+
+In addition to the usual tsi and ssi variables, the generated .nc file also includes a "spectral_solar_scaling" variable, which used by the ECRAD code.
+The computation of this variable uses the q_i.txt file.
+
+========================
+Th. Lurton 3rd Nov. 2022
+========================
Index: TOOLS/CMIP6_FORCING/SOLAR/ECRAD/q_i.txt
===================================================================
--- TOOLS/CMIP6_FORCING/SOLAR/ECRAD/q_i.txt	(revision 6266)
+++ TOOLS/CMIP6_FORCING/SOLAR/ECRAD/q_i.txt	(revision 6266)
@@ -0,0 +1,112 @@
+   3.65113
+   3.22777
+   2.69139
+  1.803999
+  0.680689
+0.05459027
+3.109565915
+3.023690161
+2.927294682
+2.658536836
+2.363372287
+3.788830093
+1.292562484
+0.9716177386
+0.07898652562
+0.1085012327
+0.02948349337
+0.01264108828
+3.925607787
+3.68861524
+3.355607991
+3.070245675
+5.083395247
+3.326048488
+1.192390541
+0.08781259959
+3.229298958
+ 3.2398307
+3.161025723
+2.981079742
+5.058397827
+3.400453832
+1.258349523
+0.09926272694
+   9.34081
+    8.9372
+   8.19346
+   7.39196
+   6.12127
+   5.23956
+   4.24941
+   3.20013
+  2.394979
+0.55783362
+15.07162699
+15.18318075
+14.74718102
+13.59462684
+12.26629121
+10.67700998
+8.717317104
+6.563135924
+4.943989969
+1.167176383
+22.8828322
+1.410780502
+  104.6244
+    94.796
+   40.8294
+   35.1801
+   28.6947
+   21.5751
+   14.6388
+   1.59111
+   2.35878
+ 1.4541238
+67.48451956
+60.11852508
+46.27861136
+31.19057333
+10.59461472
+1.553080993
+0.8278299959
+0.1392118059
+   42.6858
+    45.772
+   90.9534
+   88.2384
+  75.01118
+ 4.5315347
+   29.0079
+   28.4088
+   33.3382
+   21.8203
+ 15.982813
+0.93700512
+14.56973103
+11.90183812
+9.047010234
+5.774621646
+3.949449876
+1.634929398
+3.077229609
+0.1974403841
+0.3546402293
+0.3309576722
+0.6749451021
+0.5701053859
+1.094298089
+0.05499232719
+    1.3288
+   2.14018
+   1.97612
+      1.79
+   2.74219
+  1.861516
+ 0.8344899
+   0.12731
+ 0.0446988
+ 0.0307441
+ 0.0116728
+0.00165573
