% % gliderOptode_Generate_L1_Data.m % % Purpose: Generate Level 1 data mat files that include optode DO data. % % Requires: % MATLAB folder - contains util % % Authors: Adapted from William Stark by Harvey Seim and Chris Calloway % Marine Sciences Department % UNC-Chapel Hill % % Created: March 2013 % %/////////////////////////////////////////////////////////////////////////////// clear all; % add paths for required files... addpath('MATLAB/util/'); % populate arrays for the deployment start and end dates... % ex. strStart(2, 3) is start date for Ramses, Deployment 3 strStart = {'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'; ... '26-Jan-2012', '16-Feb-2012', '16-Mar-2012'}; strEnd = {'14-Feb-2012', '08-Mar-2012', '04-Apr-2012'; ... '14-Feb-2012', '12-Mar-2012', '03-Apr-2012'}; % SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ... for gliderIndex=1:2 % SET THE DEPLOYMENT NUMBER (1, 2 or 3) ... for deploymentNumber=1:3 clearvars -except gliderIndex deploymentNumber strStart strEnd; % glider name string... if (gliderIndex==1) strGliderName = 'Pelagia'; else strGliderName = 'Ramses'; end % deployment number string... strDeploymentNumber = num2str(deploymentNumber); % deployment start date string... strStartDate = strStart(gliderIndex, deploymentNumber); % deployment end date string... strEndDate = strEnd(gliderIndex, deploymentNumber); % define the path to the glider ascii files... datadir = strcat('GLIDER_DATA_LEVEL0/', strGliderName, ... '_Deployment', strDeploymentNumber, '/'); % define optode foil calibration coefficients % and oxygen bounds values by glider and deployment switch gliderIndex case 1 % Pelagia switch deploymentNumber case {1, 2} % Deployment 1 or 2; % before foil replacement; % 21 June 2004 calibration date; % batch 2204 C0coef=[3.12899E+03,-1.05764E+02,2.06640E+00,-1.68983E-02]; C1coef=[-1.67671E+02,4.80582E+00,-8.73323E-02,6.61507E-04]; C2coef=[3.73685E+00,-8.78300E-02,1.47560E-03,-9.96701E-06]; C3coef=[-3.96052E-02,7.46930E-04,-1.17804E-05,6.677619E-08]; C4coef=[1.61999E-04,-2.37870E-06,3.63223E-08,-1.62194E-10]; lo_oxy = 50; hi_oxy = 2000; lo_sat = 50; hi_sat = 2000; case 3 % Deployment 3; % replaced foil on 12 March 2012; % 18 August 2010 calibration date; % batch 1023 C0coef=[4.27019336E+03,-1.32723585E+02,2.15629751E+00,-1.40275831E-02]; C1coef=[-2.29729690E+02,5.74242078E+00,-6.85357898E-02,1.88612346E-04]; C2coef=[5.06401550E+00,-9.62084932E-02,5.22180779E-04,7.70889717E-06]; C3coef=[-5.26332308E-02,7.15467419E-04,3.31185072E-06,-1.86124024E-07]; C4coef=[2.10916841E-04,-1.84087896E-06,-4.28645540E-08,1.11120317E-09]; lo_oxy = 100; hi_oxy = 300; lo_sat = 20; hi_sat = 140; end case 2 % Ramses all deployments; % 2 June 2010 calibration date; % batch 5009 C0coef=[4.53793e3 -1.62595e2 3.29574 -2.79285e-2]; C1coef=[-2.50953e2 8.02322 -1.58398e-1 1.31141e-3]; C2coef=[5.66417 -1.59647e-1 3.07910e-3 -2.46265e-5]; C3coef=[-5.99449e-2 1.48326e-3 -2.82110e-5 2.15156e-7]; C4coef=[2.43614e-4 -5.26759e-6 1.00064e-7 -7.14320e-10]; lo_oxy = 10; hi_oxy = 250; lo_sat = 20; hi_sat = 100; end disp(['Processing DO for ', strGliderName, ' Deployment ', strDeploymentNumber]); disp(['lo_oxy filter = ', num2str(lo_oxy)]); disp(['hi_oxy filter = ', num2str(hi_oxy)]); disp(['lo_sat filter = ', num2str(lo_sat)]); disp(['hi_sat filter = ', num2str(hi_sat)]); %*** READ IN EBD DATA **** % declare arrays for accumulating data oxyw_oxygen = []; oxyw_saturation = []; oxyw_temp = []; oxyw_dphase = []; % oxyw_bphase = []; % oxyw_rphase = []; % oxyw_bamp = []; % oxyw_bpot = []; % oxyw_ramp = []; % oxyw_rawtemp = []; % oxyw_time = []; % oxyw_installed = []; ptime_ebd = []; % try to load all *.ebdasc files at once... [files, Dstruct] = wilddir(datadir, '.ebdasc'); nfile = size(files, 1); for i=1:nfile-1 % protect against empty ebd file if(Dstruct(i).bytes>0) data = read_gliderasc2([datadir, files(i,:)]); % if the number of values (in data.data) is less than the number % of vars (in data.vars), this means that the data were not % completely read in. To correct this, pad data.data with NaNs % until its length equals that of data.vars... if (length(data.data) < length(data.vars)) data.data = padarray(data.data, ... [0 length(data.vars)-length(data.data)], ... NaN, 'post'); end % concatenate variables with data... if(~isempty(data.data)) if (~isempty(strmatch('sci_oxy3835_wphase_oxygen',data.vars,'exact'))) oxyw_oxygen = [oxyw_oxygen; ... data.data(:,strmatch('sci_oxy3835_wphase_oxygen',... data.vars, 'exact'))]; oxyw_saturation = [oxyw_saturation; ... data.data(:,strmatch('sci_oxy3835_wphase_saturation',... data.vars, 'exact'))]; oxyw_temp = [oxyw_temp; ... data.data(:,strmatch('sci_oxy3835_wphase_temp',... data.vars, 'exact'))]; oxyw_dphase = [oxyw_dphase; ... data.data(:,strmatch('sci_oxy3835_wphase_dphase',... data.vars, 'exact'))]; % oxyw_bphase = [oxyw_bphase; ... % data.data(:,strmatch('sci_oxy3835_wphase_bphase',... % data.vars, 'exact'))]; % oxyw_rphase = [oxyw_rphase; ... % data.data(:,strmatch('sci_oxy3835_wphase_rphase',... % data.vars, 'exact'))]; % oxyw_bamp = [oxyw_bamp; ... % data.data(:,strmatch('sci_oxy3835_wphase_bamp',... % data.vars, 'exact'))]; % oxyw_bpot = [oxyw_bpot; ... % data.data(:,strmatch('sci_oxy3835_wphase_bpot',... % data.vars, 'exact'))]; % oxyw_ramp = [oxyw_ramp; ... % data.data(:,strmatch('sci_oxy3835_wphase_ramp',... % data.vars, 'exact'))]; % oxyw_rawtemp = [oxyw_rawtemp; ... % data.data(:,strmatch('sci_oxy3835_wphase_rawtemp',... % data.vars, 'exact'))]; % oxyw_time = [oxyw_time; ... % data.data(:,strmatch('sci_oxy3835_wphase_timestamp',... % data.vars, 'exact'))]; % oxyw_installed = [oxyw_installed; ... % data.data(:,strmatch('sci_oxy3835_wphase_is_installed',... % data.vars, 'exact'))]; ptime_ebd = [ptime_ebd; ... data.data(:,strmatch('sci_m_present_time',... data.vars, 'exact'))]; end end data = []; end end %*** END READ IN EBD DATA **** % remove nans from EBD data... i = find(~isnan(oxyw_dphase)); oxyw_dphase = oxyw_dphase(i); ptime_ebd = ptime_ebd(i); oxyw_oxygen = oxyw_oxygen(i); oxyw_saturation = oxyw_saturation(i); oxyw_temp = oxyw_temp(i); % oxyw_bphase = oxyw_bphase(i); % oxyw_rphase = oxyw_rphase(i); % oxyw_bamp = oxyw_bamp(i); % oxyw_bpot = oxyw_bpot(i); % oxyw_ramp = oxyw_ramp(i); % oxyw_rawtemp = oxyw_rawtemp(i); % oxyw_time = oxyw_time(i); % oxyw_installed = oxyw_installed(i); % remove unreasonably low oxygen values from EBD data... i = find(gt(oxyw_oxygen, lo_oxy)); oxyw_oxygen = oxyw_oxygen(i); ptime_ebd = ptime_ebd(i); oxyw_saturation = oxyw_saturation(i); oxyw_dphase = oxyw_dphase(i); oxyw_temp = oxyw_temp(i); % oxyw_bphase = oxyw_bphase(i); % oxyw_rphase = oxyw_rphase(i); % oxyw_bamp = oxyw_bamp(i); % oxyw_bpot = oxyw_bpot(i); % oxyw_ramp = oxyw_ramp(i); % oxyw_rawtemp = oxyw_rawtemp(i); % oxyw_time = oxyw_time(i); % oxyw_installed = oxyw_installed(i); % remove unreasonably high oxygen values from EBD data... i = find(lt(oxyw_oxygen, hi_oxy)); oxyw_oxygen = oxyw_oxygen(i); ptime_ebd = ptime_ebd(i); oxyw_saturation = oxyw_saturation(i); oxyw_dphase = oxyw_dphase(i); oxyw_temp = oxyw_temp(i); % oxyw_bphase = oxyw_bphase(i); % oxyw_rphase = oxyw_rphase(i); % oxyw_bamp = oxyw_bamp(i); % oxyw_bpot = oxyw_bpot(i); % oxyw_ramp = oxyw_ramp(i); % oxyw_rawtemp = oxyw_rawtemp(i); % oxyw_time = oxyw_time(i); % oxyw_installed = oxyw_installed(i); % remove unreasonably low saturation values from EBD data... i = find(gt(oxyw_saturation, lo_sat)); oxyw_saturation = oxyw_saturation(i); ptime_ebd = ptime_ebd(i); oxyw_oxygen = oxyw_oxygen(i); oxyw_dphase = oxyw_dphase(i); oxyw_temp = oxyw_temp(i); % oxyw_bphase = oxyw_bphase(i); % oxyw_rphase = oxyw_rphase(i); % oxyw_bamp = oxyw_bamp(i); % oxyw_bpot = oxyw_bpot(i); % oxyw_ramp = oxyw_ramp(i); % oxyw_rawtemp = oxyw_rawtemp(i); % oxyw_time = oxyw_time(i); % oxyw_installed = oxyw_installed(i); % remove unreasonably high saturation values from EBD data... i = find(lt(oxyw_saturation, hi_sat)); oxyw_saturation = oxyw_saturation(i); ptime_ebd = ptime_ebd(i); oxyw_oxygen = oxyw_oxygen(i); oxyw_dphase = oxyw_dphase(i); oxyw_temp = oxyw_temp(i); % oxyw_bphase = oxyw_bphase(i); % oxyw_rphase = oxyw_rphase(i); % oxyw_bamp = oxyw_bamp(i); % oxyw_bpot = oxyw_bpot(i); % oxyw_ramp = oxyw_ramp(i); % oxyw_rawtemp = oxyw_rawtemp(i); % oxyw_time = oxyw_time(i); % oxyw_installed = oxyw_installed(i); % apply the sort() function to ptime_oxy % to make sure it increases monotonically... [Y,I] = sort(ptime_ebd); ptime_ebd = Y; oxyw_oxygen = oxyw_oxygen(I); oxyw_saturation = oxyw_saturation(I); oxyw_temp = oxyw_temp(I); oxyw_dphase = oxyw_dphase(I); % oxyw_bphase = oxyw_bphase(I); % oxyw_rphase = oxyw_rphase(I); % oxyw_bamp = oxyw_bamp(I); % oxyw_bpot = oxyw_bpot(I); % oxyw_ramp = oxyw_ramp(I); % oxyw_rawtemp = oxyw_rawtemp(I); % oxyw_time = oxyw_time(I); % oxyw_installed = oxyw_installed(I); % load CTD data from existing mat file... strMatFilePath = strcat('GLIDER_DATA_LEVEL1/', strGliderName, ... '_Deployment', strDeploymentNumber, '_CTD_L1.mat'); load(strMatFilePath); % interpolate CTD to times of oxyw_dphase tempi = interp1(ptime,temp,ptime_ebd); salini = interp1(ptime,salinCorrected,ptime_ebd); depthi = interp1(ptime,depth,ptime_ebd); % first, implement temperature-dependent correction to DO concentration % that utilizes dphase as the input (manual page 30), coefficients from cal % sheets (a 5x4 matrix of values) and will be glider (and foil) dependent. C0 = C0coef(1) + (C0coef(2) .* tempi) + ... (C0coef(3) .* (tempi.^2)) + (C0coef(4) .* (tempi.^3)); C1 = C1coef(1) + (C1coef(2) .* tempi) + ... (C1coef(3) .* (tempi.^2)) + (C1coef(4) .* (tempi.^3)); C2 = C2coef(1) + (C2coef(2) .* tempi) + ... (C2coef(3) .* (tempi.^2)) + (C2coef(4) .* (tempi.^3)); C3 = C3coef(1) + (C3coef(2) .* tempi) + ... (C3coef(3) .* (tempi.^2)) + (C3coef(4) .* (tempi.^3)); C4 = C4coef(1) + (C4coef(2) .* tempi) + ... (C4coef(3) .* (tempi.^2)) + (C4coef(4) .* (tempi.^3)); o2_tcorr = C0 + (C1 .* oxyw_dphase) + (C2 .* (oxyw_dphase.^2)) + ... (C3 .* (oxyw_dphase.^3)) + (C4 .* (oxyw_dphase.^4)); % second, implement the salinity correction to DO concentration (page 31) % ASSUMES DEFAULT SALINITY IN SENSOR WAS SET TO ZERO, OTHERWISE ANOHTER % SCALING IS REQUIRED B0=-6.24097e-3; B1=-6.93498e-3; B2=-6.90358e-3; B3=-4.29155e-3; C0=-3.11680e-7; Ts = log((298.15 - tempi) ./ (273.15 + tempi)); o2_tscorr = o2_tcorr .* ... exp((salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + ... (B3 .* (Ts.^3)))) + ... (C0 .* (salini.^2))); % third, implement the pressure correction to DO concentration (page 32) o2_tspcorr = o2_tscorr .* (1 + (0.04 .* depthi ./ 1000)); % use polynomial to calculate DO satuations using the measured temp and % sal (manual page 30) A0 = 2.00856; A1 = 3.22400; A2 = 3.99063; A3 = 4.80299; A4 = 9.78188e-1; A5 = 1.71069; B0 = -6.24097e-3; B1 = -6.93498e-3; B2 = -6.90358e-3; B3 = -4.29155e-3; C0 = -3.11680e-7; % need interpolated temp, salinity, pressure at times of DO obs rslt = (A0 + (A1 .* Ts) + (A2 .* (Ts.^2)) + (A3 .* (Ts.^3)) + ... (A4 .* (Ts.^4)) + (A5 .* (Ts.^5))) + ... (salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + (B3 .*(Ts.^3)))) + ... (C0 .* (salini.^2)); o2_sol = exp(rslt); o2_sat = (o2_tspcorr .* 2.2414) ./ o2_sol; % convert ptime into datenum style... ptime_ebd_datenum = (ptime_ebd/3600/24) + datenum(1970, 1, 1, 0, 0, 0); % create configuration struct... units = struct( ... 'oxyw_oxygen', '10e-6 mol/dm3', ... 'oxyw_saturation', 'percent', ... 'oxyw_temp', 'degrees C', ... 'oxyw_dphase', 'degrees', ... ... % 'oxyw_bphase', 'degrees', ... ... % 'oxyw_rphase', 'degrees', ... ... % 'oxyw_bamp', 'mA', ... ... % 'oxyw_bpot', 'mV', ... ... % 'oxyw_ramp', 'mA', ... ... % 'oxyw_rawtemp', 'degrees C', ... ... % 'oxyw_time', 'timestamp', ... ... % 'oxyw_installed', 'bool', ... 'ptime_ebd', 'seconds since 0000-01-01T00:00', ... 'ptime_ebd_datenum', 'days since 1970-01-01T00:00', ... 'tempi', 'degrees C', ... 'salini', 'psu', ... 'depthi', 'meters', ... 'o2_tcorr', '10e-6 mol/dm3', ... 'o2_tscorr', '10e-6 mol/dm3', ... 'o2_tspcorr', '10e-6 mol/dm3', ... 'o2_sol', 'cm3/liter at 1031 hPa', ... 'o2_sat', 'percent'); variable_description = struct( ... 'oxyw_oxygen', 'dissolved oxygen', ... 'oxyw_saturation', 'dissolved oxygen saturation', ... 'oxyw_temp', 'water temperature', ... 'oxyw_dphase', 'phase difference', ... ... % 'oxyw_bphase', 'blue phase', ... ... % 'oxyw_rphase', 'red phase', ... ... % 'oxyw_bamp', 'blue current bias', ... ... % 'oxyw_bpot', 'blue voltage bias', ... ... % 'oxyw_ramp', 'red current bias', ... ... % 'oxyw_rawtemp', 'raw water temperature', ... ... % 'oxyw_time', 'optode timestampe', ... ... % 'oxyw_installed', 'bool', ... 'ptime_ebd', 'science computer time', ... 'ptime_ebd_datenum', 'science computer date', ... 'tempi', 'interpolated CTD water temperature', ... 'salini', 'interpoloted CTD salinity', ... 'depthi', 'interpolated CTD depth', ... 'o2_tcorr', 'temperature corrected dissolved oxygen', ... 'o2_tscorr', 'temperature and salinity corrected dissolved oxygen', ... 'o2_tspcorr', 'temperature, salinity, and pressure corrected dissolved oxygen', ... 'o2_sol', 'corrected oxygen solubility', ... 'o2_sat', 'corrected oxygen saturation'); correction_coefficients = struct('C0coef', C0coef, ... 'C1coef', C1coef, ... 'C2coef', C2coef, ... 'C3coef', C3coef, ... 'C4coef', C4coef); config = struct('glider_name', strGliderName,... 'deployment_number', strDeploymentNumber,... 'start_date', strStartDate,... 'end_date', strEndDate,... 'correction_coefficients', correction_coefficients,... 'var_descriptions', variable_description,... 'var_units', units); % set Level 1 data mat file name... strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_DO_L1.mat'); % save flight data to mat file... save(strMatFileName,... 'config', ... 'oxyw_oxygen', ... 'oxyw_saturation', ... 'oxyw_temp', ... 'oxyw_dphase', ... ... % 'oxyw_bphase', ... ... % 'oxyw_rphase', ... ... % 'oxyw_bamp', ... ... % 'oxyw_bpot', ... ... % 'oxyw_ramp', ... ... % 'oxyw_rawtemp', ... ... % 'oxyw_time', ... ... % 'oxyw_installed', ... 'ptime_ebd', ... 'ptime_ebd_datenum', ... 'tempi', ... 'salini', ... 'depthi', ... 'o2_tcorr', ... 'o2_tscorr', ... 'o2_tspcorr', ... 'o2_sol', ... 'o2_sat'); end end