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