source: TOOLS/CMIP6_FORCING/SOLAR/CMIP6_piControl_2bands.m @ 3595

Last change on this file since 3595 was 3595, checked in by tlurton, 4 years ago

20/02/2018 Th. Lurton
Updated paths to final location of input files on Ciclad: /prodigfs/project/input4MIPs/.

File size: 12.4 KB
Line 
1% Updated 05/02/2018, ThL
2% * changed paths to retrieve new version of input data;
3% * version considered is now 3.2.
4
5% MATLAB run script to integrate spectral solar irradiance (SSI) data
6% recommended and provided by SPARC/SOLARIS-HEPPA for use in CMIP6 to
7% spectral bands of any desired resolution (according to
8% radiation/photolysis scheme of employed (chemistry) climate model)
9%
10% Author: Tim Kruschke (tkruschke[at]geomar.de)
11% Last modified: 02 May 2016 (compliance with solarforcing_*_*_3.1.nc)
12%
13% Description:
14%    1. reading data from netcdf-file as provided by SPARC/SOLARIS-HEPPA
15%    2. reading definitions of spectral bands (user-specified by creating an
16%       ascii-textfile according to specific formatting, as described below)
17%    3. calling MATLAB-function "regrid_ssi" (provided with this run script as
18%       part of a joint zip-archive) to integrate SSI over specified bands
19%    4. writing new netcdf-file with same structure and attributes as original
20%       (input) file, but only with (unchanged) total solar irradiance
21%       (TSI), SSI integrated over spectral bands specified by user, and
22%       irradiance fraction (of TSI) in respective band (plus
23%       necessary/related dimensions/variables time, wlen, wlenbinsize)
24%
25% Formatting of ascii-textfile, containing definitions of spectral bands:
26%    NO HEADER LINES, NO UNITS, NO CHARACTER DELIMITERS
27%    It's as simple as this:
28%    - one line per spectral band
29%    - each line contains two floating point numbers: lower and upper edge
30%      of respective spectral band (in nanometers!!!), separated by
31%      blank space(s) or tabulator(s)
32%    See example_spectral_bands.dat (provided with this run script as
33%    part of a joint zip-archive) as an example.
34%--------------------------------------------------------------------------
35
36%% User-defined paths and filenames
37
38% path to directory of input data; adapt to your file structure
39dir_in=['.'];
40
41% path to directory of output data; adapt to your file structure if desired
42% to differ from input directory
43dir_out=dir_in;
44
45% filename of ascii-textfile (to be placed in dir_in) containing
46% definitions of spectral bands (see explanations in script header);
47% adapt to name of your band definition file
48%filename_spectral_bands='example_spectral_bands.dat';
49%filename_spectral_bands='example_spectral_bands_6_SSI_FRAC=1.dat';
50filename_spectral_bands='example_spectral_bands_2_SSI_FRAC=1.dat';
51
52% filename of original data file, provided by SPARC/SOLARIS-HEPPA (to be
53% placed in dir_in); adapt to respective file that is to be considered
54%filename_in='solarforcing_ref_mon_3.1.nc';
55%filename_in='solarforcing_ref_day_3.1.nc';
56filename_in='/prodigfs/project/input4MIPs/CMIP6/CMIP/SOLARIS-HEPPA/SOLARIS-HEPPA-3-2/atmos/fx/multiple/gz/v20170103/solarforcing-picontrol-fx_input4MIPs_solar_CMIP_SOLARIS-HEPPA-3-2_gn_18500101-18730128.nc'
57
58% lines 53-54 only used to match input and output filenames
59% ATTENTION: works only on unix platforms; if run on other platforms,
60% delete/comment following two lines and enter complete filename_out in
61% l. 59
62[~,filename_in_base]=unix(['basename ' dir_in filesep filename_in ' .nc']);
63filename_in_base=strcat(filename_in_base);
64
65% filename of new file to be created, containing SSI data in spectral
66% resolutiuon as desired by the user; adapt to whatever name you like if
67% not desired to match input filename
68%filename_out=[filename_in_base '_SSI_regridded.nc'];
69filename_out=[filename_in_base '_SSI_regridded_2_TSI_anomaly.nc'];
70
71%% Usually, nothing should be edited below this line!!!!!!!!!!!!!!!!!!!!!!!
72
73ncid_in=netcdf.open([dir_in filesep filename_in],'NOWRITE');
74ncid_out=netcdf.create(strcat([dir_out filesep filename_out]),'CLOBBER'); % existing files are overwritten
75
76% adapting title of original file
77nc_title=netcdf.getAtt(ncid_in,-1,'title');
78nc_title=['Spectral solar irradiance (SSI) of ' nc_title ', converted to custom spectral resolution'];
79netcdf.putAtt(ncid_out,-1,'title',nc_title);
80
81%% time stamp
82netcdf.putAtt(ncid_out,-1,'creation_date',datestr(now)); % creates a time stamp, though not compliant with ISO 8601
83%c% You can use the following lines (commented out) instead to generate a
84%c% time stamp compliant to ISO 8601; change last digits according to
85%c% your time zone in this case (currently +02:00 for CEST)
86%     create_date_vec=datevec(now);
87%     create_date_str=cell(1,6);
88%     for id=1:6;
89%         if create_date_vec(id)<10;
90%             create_date_str{1,id}=['0' num2str(round(create_date_vec(id)))];
91%         else
92%             create_date_str{1,id}=num2str(round(create_date_vec(id)));
93%         end
94%     end
95%     netcdf.putAtt(ncid_out,-1,'creation_date',[create_date_str{1,1} '-' create_date_str{1,2} '-' ...
96%         create_date_str{1,3} 'T' create_date_str{1,4} ':' create_date_str{1,5} ':' create_date_str{1,6} '+02:00'])
97%%
98
99% copy global attributes to new file
100netcdf.copyAtt(ncid_in,-1,'institution',ncid_out,-1)
101netcdf.copyAtt(ncid_in,-1,'institution_id',ncid_out,-1)
102netcdf.copyAtt(ncid_in,-1,'activity_id',ncid_out,-1)
103netcdf.copyAtt(ncid_in,-1,'comment',ncid_out,-1)
104%netcdf.copyAtt(ncid_in,-1,'time_coverage_start',ncid_out,-1)
105%netcdf.copyAtt(ncid_in,-1,'time_coverage_end',ncid_out,-1)
106netcdf.copyAtt(ncid_in,-1,'frequency',ncid_out,-1)
107netcdf.copyAtt(ncid_in,-1,'source',ncid_out,-1)
108netcdf.copyAtt(ncid_in,-1,'source_id',ncid_out,-1)
109netcdf.copyAtt(ncid_in,-1,'realm',ncid_out,-1)
110netcdf.copyAtt(ncid_in,-1,'further_info_url',ncid_out,-1)
111netcdf.copyAtt(ncid_in,-1,'metadata_url',ncid_out,-1)
112netcdf.copyAtt(ncid_in,-1,'contributor_name',ncid_out,-1)
113netcdf.copyAtt(ncid_in,-1,'references',ncid_out,-1)
114netcdf.copyAtt(ncid_in,-1,'Conventions',ncid_out,-1)
115netcdf.putAtt(ncid_out,-1,'history',['Original SSI data from file ' filename_in ' integrated ' ...
116    'over spectral bands as defined in file ' filename_spectral_bands ' ' ...
117    'by MATLAB run script CMIP6_SSI_convert_resolution_run.m, using function regrid_ssi.m']);
118
119% reading SSI-related variables from CMIP6 input file
120wvl_in_ID=netcdf.inqVarID(ncid_in,'wlen');
121wvl_edges_in_ID=netcdf.inqVarID(ncid_in,'wlen_bnds');
122bw_in_ID=netcdf.inqVarID(ncid_in,'wlenbinsize');
123time_in_ID=netcdf.inqVarID(ncid_in,'time');
124time_bnds_in_ID=netcdf.inqVarID(ncid_in,'time_bnds');
125year_in_ID=netcdf.inqVarID(ncid_in,'calyear');
126month_in_ID=netcdf.inqVarID(ncid_in,'calmonth');
127day_in_ID=netcdf.inqVarID(ncid_in,'calday');
128tsi_in_ID=netcdf.inqVarID(ncid_in,'tsi');
129ssi_in_ID=netcdf.inqVarID(ncid_in,'ssi');
130
131wvl_in=netcdf.getVar(ncid_in,wvl_in_ID);
132wvl_edges_in=(netcdf.getVar(ncid_in,wvl_edges_in_ID))';
133bw_in=netcdf.getVar(ncid_in,bw_in_ID);
134ssi_in=netcdf.getVar(ncid_in,ssi_in_ID);
135tsi_in=netcdf.getVar(ncid_in,tsi_in_ID);
136time_in=netcdf.getVar(ncid_in,time_in_ID);
137time_bnds_in=netcdf.getVar(ncid_in,time_bnds_in_ID);
138year_in=netcdf.getVar(ncid_in,year_in_ID);
139month_in=netcdf.getVar(ncid_in,month_in_ID);
140day_in=netcdf.getVar(ncid_in,day_in_ID);
141
142% read user specified definitions of spectral bands
143fid=fopen([dir_in filesep filename_spectral_bands]);
144tmp=textscan(fid,'%f %f');
145fclose(fid);
146wvl_edges_out=tmp{1,1};
147wvl_edges_out(:,2)=tmp{1,2};
148clear tmp
149
150% calculation of output central wave lengths
151wvl_out=sum(wvl_edges_out,2)/2;
152
153% calculation of output bin width
154bw_out=wvl_edges_out(:,2)-wvl_edges_out(:,1);
155
156% check for consistency of output bin definitions
157if min(bw_out)<0;
158    error(['Error: For at least one spectral band, the specified upper edge \n' ...
159        '(right column) in your input file ' filename_in '\n ' ...
160        'is smaller than the respective lower edge (left column).'])
161end
162
163% check for spectral coverage of input data with respect to desired output
164% bins
165if min(wvl_edges_out(:,1))<min(wvl_edges_in(:,1)) || max(wvl_edges_out(:,2))>max(wvl_edges_in(:,2));
166    warning(['You are extrapolating into spectral ranges not covered by the input spectrum! ' ...
167        'This routine supports this by linearly extrapolating based on change between last two bins. ' ...
168        'However, this might be physically inappropriate!']);
169end
170
171% Integration of SSI data for each specified spectral band
172% (regrid_ssi.m calculates weighted means of SSI per nm, this is multiplied by binwidth)
173ssi_out=bsxfun(@times,regrid_ssi(ssi_in,wvl_edges_in,wvl_edges_out),bw_out);
174
175% Calculation of solar irradiance fraction (compared to total solar
176% irradiance) in each band
177ssi_frac_out=bsxfun(@rdivide,ssi_out,tsi_in');
178
179% Calculation the average value of tsi
180avtsi=mean(tsi_in)
181disp (avtsi)
182
183%Calculate the anomaly of tsi
184antsi(:,1)=tsi_in(:,1)-avtsi
185
186disp (antsi)
187
188format longEng
189disp (avtsi)
190
191% writing data (including dimensions) to output file
192  % dimension definitions
193time_dimID=netcdf.defDim(ncid_out,'time',0);
194wvl_dimID=netcdf.defDim(ncid_out,'wlen',length(wvl_out));
195bounds_dimID=netcdf.defDim(ncid_out,'nbd',2);
196
197  % variable definitions
198time_out_ID=netcdf.defVar(ncid_out,'time','double',time_dimID);
199time_bnds_out_ID=netcdf.defVar(ncid_out,'time_bnds','double',[bounds_dimID time_dimID]);
200wvl_out_ID=netcdf.defVar(ncid_out,'wlen','double',wvl_dimID);
201wvl_edges_out_ID=netcdf.defVar(ncid_out,'wlen_bnds','double',[wvl_dimID bounds_dimID]);
202bw_out_ID=netcdf.defVar(ncid_out,'wlenbinsize','double',wvl_dimID);
203year_out_ID=netcdf.defVar(ncid_out,'calyear','double',time_dimID);
204month_out_ID=netcdf.defVar(ncid_out,'calmonth','double',time_dimID);
205day_out_ID=netcdf.defVar(ncid_out,'calday','double',time_dimID);
206tsi_out_ID=netcdf.defVar(ncid_out,'tsi','float',time_dimID);
207antsi_out_ID=netcdf.defVar(ncid_out,'antsi','float',time_dimID);
208ssi_out_ID=netcdf.defVar(ncid_out,'ssi','float',[wvl_dimID time_dimID]);
209ssi_frac_out_ID=netcdf.defVar(ncid_out,'ssi_frac','float',[wvl_dimID time_dimID]);
210
211  % copy variable attributes from input to output file
212netcdf.copyAtt(ncid_in,time_in_ID,'axis',ncid_out,time_out_ID)
213netcdf.copyAtt(ncid_in,time_in_ID,'long_name',ncid_out,time_out_ID)
214netcdf.copyAtt(ncid_in,time_in_ID,'standard_name',ncid_out,time_out_ID)
215netcdf.copyAtt(ncid_in,time_in_ID,'units',ncid_out,time_out_ID)
216netcdf.copyAtt(ncid_in,time_in_ID,'calendar',ncid_out,time_out_ID)
217netcdf.copyAtt(ncid_in,time_in_ID,'_CoordinateAxisType',ncid_out,time_out_ID)
218netcdf.copyAtt(ncid_in,time_in_ID,'bounds',ncid_out,time_out_ID)
219
220netcdf.copyAtt(ncid_in,time_bnds_in_ID,'long_name',ncid_out,time_bnds_out_ID)
221netcdf.copyAtt(ncid_in,time_bnds_in_ID,'units',ncid_out,time_bnds_out_ID)
222
223netcdf.copyAtt(ncid_in,wvl_in_ID,'standard_name',ncid_out,wvl_out_ID)
224netcdf.copyAtt(ncid_in,wvl_in_ID,'long_name',ncid_out,wvl_out_ID)
225netcdf.copyAtt(ncid_in,wvl_in_ID,'units',ncid_out,wvl_out_ID)
226netcdf.copyAtt(ncid_in,wvl_in_ID,'bounds',ncid_out,wvl_out_ID)
227
228netcdf.copyAtt(ncid_out,wvl_edges_in_ID,'long_name',ncid_out,wvl_edges_out_ID)
229netcdf.copyAtt(ncid_out,wvl_edges_in_ID,'units',ncid_out,wvl_edges_out_ID)
230
231netcdf.copyAtt(ncid_in,bw_in_ID,'long_name',ncid_out,bw_out_ID)
232netcdf.copyAtt(ncid_in,bw_in_ID,'units',ncid_out,bw_out_ID)
233
234netcdf.copyAtt(ncid_in,year_in_ID,'long_name',ncid_out,year_out_ID)
235
236netcdf.copyAtt(ncid_in,month_in_ID,'long_name',ncid_out,month_out_ID)
237
238netcdf.copyAtt(ncid_in,day_in_ID,'long_name',ncid_out,day_out_ID)
239
240netcdf.copyAtt(ncid_in,tsi_in_ID,'standard_name',ncid_out,tsi_out_ID)
241netcdf.copyAtt(ncid_in,tsi_in_ID,'long_name',ncid_out,tsi_out_ID)
242netcdf.copyAtt(ncid_in,tsi_in_ID,'units',ncid_out,tsi_out_ID)
243netcdf.copyAtt(ncid_in,tsi_in_ID,'cell_methods',ncid_out,tsi_out_ID)
244
245netcdf.putAtt(ncid_out,antsi_out_ID,'standard_name','total solar_irradiance anomaly')
246netcdf.putAtt(ncid_out,antsi_out_ID,'long_name','reconstructed total solar irradiance anomaly at 1 AU')
247netcdf.putAtt(ncid_out,antsi_out_ID,'units','W m^-2')
248
249netcdf.putAtt(ncid_out,ssi_out_ID,'long_name','reconstructed solar irradiance at 1 AU within respective bands')
250netcdf.putAtt(ncid_out,ssi_out_ID,'units','W m^-2')
251
252netcdf.putAtt(ncid_out,ssi_frac_out_ID,'long_name','reconstructed fraction of total solar irradiance at 1 AU within respective bands')
253netcdf.putAtt(ncid_out,ssi_frac_out_ID,'units','1')
254
255netcdf.close(ncid_in);
256netcdf.endDef(ncid_out);
257
258netcdf.putVar(ncid_out,time_out_ID,0,length(time_in),time_in)
259netcdf.putVar(ncid_out,time_bnds_out_ID,time_bnds_in)
260netcdf.putVar(ncid_out,wvl_out_ID,wvl_out)
261netcdf.putVar(ncid_out,wvl_edges_out_ID,wvl_edges_out)
262netcdf.putVar(ncid_out,bw_out_ID,bw_out)
263netcdf.putVar(ncid_out,year_out_ID,year_in)
264netcdf.putVar(ncid_out,month_out_ID,month_in)
265netcdf.putVar(ncid_out,day_out_ID,day_in)
266netcdf.putVar(ncid_out,tsi_out_ID,tsi_in)
267netcdf.putVar(ncid_out,antsi_out_ID,antsi)
268netcdf.putVar(ncid_out,ssi_out_ID,ssi_out)
269netcdf.putVar(ncid_out,ssi_frac_out_ID,ssi_frac_out)
270
271netcdf.sync(ncid_out);
272netcdf.close(ncid_out);
Note: See TracBrowser for help on using the repository browser.