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