% Updated 05/02/2018, ThL % * changed paths to retrieve new version of input data; % * version considered is now 3.2. % 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='example_spectral_bands_6_SSI_FRAC=1.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='/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' % 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 % resolutiuon as desired by the user; adapt to whatever name you like if % not desired to match input filename %filename_out=[filename_in_base '_SSI_regridded.nc']; filename_out=[filename_in_base '_SSI_regridded_6_TSI_anomaly.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))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) % 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]); % 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.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.sync(ncid_out); netcdf.close(ncid_out);