1 |
% |
---|
2 |
% gliderOptode_Generate_L1_Data.m |
---|
3 |
% |
---|
4 |
% Purpose: Generate Level 1 data mat files that include optode DO data. |
---|
5 |
% |
---|
6 |
% Requires: |
---|
7 |
% MATLAB folder - contains util |
---|
8 |
% |
---|
9 |
% Authors: Adapted from William Stark by Harvey Seim and Chris Calloway |
---|
10 |
% Marine Sciences Department |
---|
11 |
% UNC-Chapel Hill |
---|
12 |
% |
---|
13 |
% Created: March 2013 |
---|
14 |
% |
---|
15 |
% |
---|
16 |
|
---|
17 |
clear all; |
---|
18 |
|
---|
19 |
% add paths for required files... |
---|
20 |
addpath('MATLAB/util/'); |
---|
21 |
|
---|
22 |
% populate arrays for the deployment start and end dates... |
---|
23 |
% ex. strStart(2, 3) is start date for Ramses, Deployment 3 |
---|
24 |
strStart = {'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'; ... |
---|
25 |
'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'}; |
---|
26 |
strEnd = {'14-Feb-2012', '08-Mar-2012', '04-Apr-2012'; ... |
---|
27 |
'14-Feb-2012', '12-Mar-2012', '03-Apr-2012'}; |
---|
28 |
|
---|
29 |
% SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ... |
---|
30 |
for gliderIndex=1:2 |
---|
31 |
|
---|
32 |
% SET THE DEPLOYMENT NUMBER (1, 2 or 3) ... |
---|
33 |
for deploymentNumber=1:3 |
---|
34 |
|
---|
35 |
clearvars -except gliderIndex deploymentNumber strStart strEnd; |
---|
36 |
|
---|
37 |
% glider name string... |
---|
38 |
if (gliderIndex==1) |
---|
39 |
strGliderName = 'Pelagia'; |
---|
40 |
else |
---|
41 |
strGliderName = 'Ramses'; |
---|
42 |
end |
---|
43 |
|
---|
44 |
% deployment number string... |
---|
45 |
strDeploymentNumber = num2str(deploymentNumber); |
---|
46 |
|
---|
47 |
% deployment start date string... |
---|
48 |
strStartDate = strStart(gliderIndex, deploymentNumber); |
---|
49 |
|
---|
50 |
% deployment end date string... |
---|
51 |
strEndDate = strEnd(gliderIndex, deploymentNumber); |
---|
52 |
|
---|
53 |
% define the path to the glider ascii files... |
---|
54 |
datadir = strcat('GLIDER_DATA_LEVEL0/', strGliderName, ... |
---|
55 |
'_Deployment', strDeploymentNumber, '/'); |
---|
56 |
|
---|
57 |
% define optode foil calibration coefficients |
---|
58 |
% by glider and deployment |
---|
59 |
switch gliderIndex |
---|
60 |
case 1 % Pelagia |
---|
61 |
switch deploymentNumber |
---|
62 |
case {1, 2} % Deployment 1 or 2; |
---|
63 |
% before foil replacement; |
---|
64 |
% 21 June 2004 calibration date; |
---|
65 |
% batch 2204 |
---|
66 |
C0coef=[3.12899E+03,-1.05764E+02,2.06640E+00,-1.68983E-02]; |
---|
67 |
C1coef=[-1.67671E+02,4.80582E+00,-8.73323E-02,6.61507E-04]; |
---|
68 |
C2coef=[3.73685E+00,-8.78300E-02,1.47560E-03,-9.96701E-06]; |
---|
69 |
C3coef=[-3.96052E-02,7.46930E-04,-1,17804E-05,6.677619E-08]; |
---|
70 |
C4coef=[1.61999E-04,-2.37870E-06,3.63223E-08,-1.62194E-10]; |
---|
71 |
case 3 % Deployment 3; |
---|
72 |
% replaced foil on 12 March 2012; |
---|
73 |
% 18 August 2010 calibration date; |
---|
74 |
% batch 1023 |
---|
75 |
C0coef=[4.27019336E+03,-1.32723585E+02,2.15629751E+00,-1.40275831E-02]; |
---|
76 |
C1coef=[-2.29729690E+02,5.74242078E+00,-6.85357898E-02,1.88612346E-04]; |
---|
77 |
C2coef=[5.06401550E+00,-9.62084932E-02,5.22180779E-04,7.70889717E-06]; |
---|
78 |
C3coef=[-5.26332308E-02,7.15467419E-04,3.31185072E-06,-1.86124024E-07]; |
---|
79 |
C4coef=[2.10916841E-04,-1.84087896E-06,-4.28645540E-08,1.11120317E-09]; |
---|
80 |
end |
---|
81 |
case 2 % Ramses all deployments; |
---|
82 |
% 2 June 2010 calibration date; |
---|
83 |
% batch 5009 |
---|
84 |
C0coef=[4.53793e3 -1.62595e2 3.29574 -2.79285e-2]; |
---|
85 |
C1coef=[-2.50953e2 8.02322 -1.58398e-1 1.31141e-3]; |
---|
86 |
C2coef=[5.66417 -1.59647e-1 3.07910e-3 -2.46265e-5]; |
---|
87 |
C3coef=[-5.99449e-2 1.48326e-3 -2.82110e-5 2.15156e-7]; |
---|
88 |
C4coef=[2.43614e-4 -5.26759e-6 1.00064e-7 -7.14320e-10]; |
---|
89 |
end |
---|
90 |
|
---|
91 |
%*** READ IN EBD DATA **** |
---|
92 |
% declare arrays for accumulating data |
---|
93 |
oxyw_oxygen = []; |
---|
94 |
oxyw_saturation = []; |
---|
95 |
oxyw_temp = []; |
---|
96 |
oxyw_dphase = []; |
---|
97 |
% oxyw_bphase = []; |
---|
98 |
% oxyw_rphase = []; |
---|
99 |
% oxyw_bamp = []; |
---|
100 |
% oxyw_bpot = []; |
---|
101 |
% oxyw_ramp = []; |
---|
102 |
% oxyw_rawtemp = []; |
---|
103 |
% oxyw_time = []; |
---|
104 |
% oxyw_installed = []; |
---|
105 |
ptime_ebd = []; |
---|
106 |
|
---|
107 |
% try to load all *.ebdasc files at once... |
---|
108 |
[files, Dstruct] = wilddir(datadir, '.ebdasc'); |
---|
109 |
nfile = size(files, 1); |
---|
110 |
|
---|
111 |
for i=1:nfile-1 |
---|
112 |
% protect against empty ebd file |
---|
113 |
if(Dstruct(i).bytes>0) |
---|
114 |
data = read_gliderasc2([datadir, files(i,:)]); |
---|
115 |
|
---|
116 |
% if the number of values (in data.data) is less than the number |
---|
117 |
% of vars (in data.vars), this means that the data were not |
---|
118 |
% completely read in. To correct this, pad data.data with NaNs |
---|
119 |
% until its length equals that of data.vars... |
---|
120 |
if (length(data.data) < length(data.vars)) |
---|
121 |
data.data = padarray(data.data, ... |
---|
122 |
[0 length(data.vars)-length(data.data)], ... |
---|
123 |
NaN, 'post'); |
---|
124 |
end |
---|
125 |
|
---|
126 |
% concatenate variables with data... |
---|
127 |
if(~isempty(data.data)) |
---|
128 |
oxyw_oxygen = [oxyw_oxygen; ... |
---|
129 |
data.data(:,strmatch('sci_oxy3835_wphase_oxygen',... |
---|
130 |
data.vars, 'exact'))]; |
---|
131 |
oxyw_saturation = [oxyw_saturation; ... |
---|
132 |
data.data(:,strmatch('sci_oxy3835_wphase_saturation',... |
---|
133 |
data.vars, 'exact'))]; |
---|
134 |
oxyw_temp = [oxyw_temp; ... |
---|
135 |
data.data(:,strmatch('sci_oxy3835_wphase_temp',... |
---|
136 |
data.vars, 'exact'))]; |
---|
137 |
oxyw_dphase = [oxyw_dphase; ... |
---|
138 |
data.data(:,strmatch('sci_oxy3835_wphase_dphase',... |
---|
139 |
data.vars, 'exact'))]; |
---|
140 |
% oxyw_bphase = [oxyw_bphase; ... |
---|
141 |
% data.data(:,strmatch('sci_oxy3835_wphase_bphase',... |
---|
142 |
% data.vars, 'exact'))]; |
---|
143 |
% oxyw_rphase = [oxyw_rphase; ... |
---|
144 |
% data.data(:,strmatch('sci_oxy3835_wphase_rphase',... |
---|
145 |
% data.vars, 'exact'))]; |
---|
146 |
% oxyw_bamp = [oxyw_bamp; ... |
---|
147 |
% data.data(:,strmatch('sci_oxy3835_wphase_bamp',... |
---|
148 |
% data.vars, 'exact'))]; |
---|
149 |
% oxyw_bpot = [oxyw_bpot; ... |
---|
150 |
% data.data(:,strmatch('sci_oxy3835_wphase_bpot',... |
---|
151 |
% data.vars, 'exact'))]; |
---|
152 |
% oxyw_ramp = [oxyw_ramp; ... |
---|
153 |
% data.data(:,strmatch('sci_oxy3835_wphase_ramp',... |
---|
154 |
% data.vars, 'exact'))]; |
---|
155 |
% oxyw_rawtemp = [oxyw_rawtemp; ... |
---|
156 |
% data.data(:,strmatch('sci_oxy3835_wphase_rawtemp',... |
---|
157 |
% data.vars, 'exact'))]; |
---|
158 |
% oxyw_time = [oxyw_time; ... |
---|
159 |
% data.data(:,strmatch('sci_oxy3835_wphase_timestamp',... |
---|
160 |
% data.vars, 'exact'))]; |
---|
161 |
% oxyw_installed = [oxyw_installed; ... |
---|
162 |
% data.data(:,strmatch('sci_oxy3835_wphase_is_installed',... |
---|
163 |
% data.vars, 'exact'))]; |
---|
164 |
ptime_ebd = [ptime_ebd; ... |
---|
165 |
data.data(:,strmatch('sci_m_present_time',... |
---|
166 |
data.vars, 'exact'))]; |
---|
167 |
end |
---|
168 |
|
---|
169 |
data = []; |
---|
170 |
end |
---|
171 |
end |
---|
172 |
%*** END READ IN EBD DATA **** |
---|
173 |
|
---|
174 |
% remove nans from EBD data... |
---|
175 |
i = find(~isnan(oxyw_dphase)); |
---|
176 |
oxyw_dphase = oxyw_dphase(i); |
---|
177 |
ptime_ebd = ptime_ebd(i); |
---|
178 |
oxyw_oxygen = oxyw_oxygen(i); |
---|
179 |
oxyw_saturation = oxyw_saturation(i); |
---|
180 |
oxyw_temp = oxyw_temp(i); |
---|
181 |
% oxyw_bphase = oxyw_bphase(i); |
---|
182 |
% oxyw_rphase = oxyw_rphase(i); |
---|
183 |
% oxyw_bamp = oxyw_bamp(i); |
---|
184 |
% oxyw_bpot = oxyw_bpot(i); |
---|
185 |
% oxyw_ramp = oxyw_ramp(i); |
---|
186 |
% oxyw_rawtemp = oxyw_rawtemp(i); |
---|
187 |
% oxyw_time = oxyw_time(i); |
---|
188 |
% oxyw_installed = oxyw_installed(i); |
---|
189 |
|
---|
190 |
% apply the sort() function to ptime_oxy |
---|
191 |
% to make sure it increases monotonically... |
---|
192 |
[Y,I] = sort(ptime_ebd); |
---|
193 |
ptime_ebd = Y; |
---|
194 |
oxyw_oxygen = oxyw_oxygen(I); |
---|
195 |
oxyw_saturation = oxyw_saturation(I); |
---|
196 |
oxyw_temp = oxyw_temp(I); |
---|
197 |
oxyw_dphase = oxyw_dphase(I); |
---|
198 |
% oxyw_bphase = oxyw_bphase(I); |
---|
199 |
% oxyw_rphase = oxyw_rphase(I); |
---|
200 |
% oxyw_bamp = oxyw_bamp(I); |
---|
201 |
% oxyw_bpot = oxyw_bpot(I); |
---|
202 |
% oxyw_ramp = oxyw_ramp(I); |
---|
203 |
% oxyw_rawtemp = oxyw_rawtemp(I); |
---|
204 |
% oxyw_time = oxyw_time(I); |
---|
205 |
% oxyw_installed = oxyw_installed(I); |
---|
206 |
|
---|
207 |
% load CTD data from existing mat file... |
---|
208 |
strMatFilePath = strcat('GLIDER_DATA_LEVEL1/', strGliderName, ... |
---|
209 |
'_Deployment', strDeploymentNumber, '_CTD_L1.mat'); |
---|
210 |
load(strMatFilePath); |
---|
211 |
|
---|
212 |
% interpolate CTD to times of oxyw_dphase |
---|
213 |
tempi = interp1(ptime,temp,ptime_ebd); |
---|
214 |
salini = interp1(ptime,salinCorrected,ptime_ebd); |
---|
215 |
depthi = interp1(ptime,depth,ptime_ebd); |
---|
216 |
|
---|
217 |
% first, implement temperature-dependent correction to DO concentration |
---|
218 |
% that utilizes dphase as the input (manual page 30), coefficients from cal |
---|
219 |
% sheets (a 5x4 matrix of values) and will be glider (and foil) dependent. |
---|
220 |
C0 = C0coef(1) + (C0coef(2) .* tempi) + ... |
---|
221 |
(C0coef(3) .* (tempi.^2)) + (C0coef(4) .* (tempi.^3)); |
---|
222 |
C1 = C1coef(1) + (C1coef(2) .* tempi) + ... |
---|
223 |
(C1coef(3) .* (tempi.^2)) + (C1coef(4) .* (tempi.^3)); |
---|
224 |
C2 = C2coef(1) + (C2coef(2) .* tempi) + ... |
---|
225 |
(C2coef(3) .* (tempi.^2)) + (C2coef(4) .* (tempi.^3)); |
---|
226 |
C3 = C3coef(1) + (C3coef(2) .* tempi) + ... |
---|
227 |
(C3coef(3) .* (tempi.^2)) + (C3coef(4) .* (tempi.^3)); |
---|
228 |
C4 = C4coef(1) + (C4coef(2) .* tempi) + ... |
---|
229 |
(C4coef(3) .* (tempi.^2)) + (C4coef(4) .* (tempi.^3)); |
---|
230 |
|
---|
231 |
o2_tcorr = C0 + (C1 .* oxyw_dphase) + (C2 .* (oxyw_dphase.^2)) + ... |
---|
232 |
(C3 .* (oxyw_dphase.^3)) + (C4 .* (oxyw_dphase.^4)); |
---|
233 |
|
---|
234 |
% second, implement the salinity correction to DO concentration (page 31) |
---|
235 |
% ASSUMES DEFAULT SALINITY IN SENSOR WAS SET TO ZERO, OTHERWISE ANOHTER |
---|
236 |
% SCALING IS REQUIRED |
---|
237 |
B0=-6.24097e-3; |
---|
238 |
B1=-6.93498e-3; |
---|
239 |
B2=-6.90358e-3; |
---|
240 |
B3=-4.29155e-3; |
---|
241 |
C0=-3.11680e-7; |
---|
242 |
|
---|
243 |
Ts = log((298.15 - tempi) ./ (273.15 + tempi)); |
---|
244 |
|
---|
245 |
o2_tscorr = o2_tcorr .* ... |
---|
246 |
exp((salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + ... |
---|
247 |
(B3 .* (Ts.^3)))) + ... |
---|
248 |
(C0 .* (salini.^2))); |
---|
249 |
|
---|
250 |
% third, implement the pressure correction to DO concentration (page 32) |
---|
251 |
o2_tspcoor = o2_tscorr .* (1 + (0.04 .* depthi ./ 1000)); |
---|
252 |
|
---|
253 |
% use polynomial to calculate DO satuations using the measured temp and |
---|
254 |
% sal (manual page 30) |
---|
255 |
A0 = 2.00856; |
---|
256 |
A1 = 3.22400; |
---|
257 |
A2 = 3.99063; |
---|
258 |
A3 = 4.80299; |
---|
259 |
A4 = 9.78188e-1; |
---|
260 |
A5 = 1.71069; |
---|
261 |
B0 = -6.24097e-3; |
---|
262 |
B1 = -6.93498e-3; |
---|
263 |
B2 = -6.90358e-3; |
---|
264 |
B3 = -4.29155e-3; |
---|
265 |
C0 = -3.11680e-7; |
---|
266 |
|
---|
267 |
% need interpolated temp, salinity, pressure at times of DO obs |
---|
268 |
rslt = (A0 + (A1 .* Ts) + (A2 .* (Ts.^2)) + (A3 .* (Ts.^3)) + ... |
---|
269 |
(A4 .* (Ts.^4)) + (A5 .* (Ts.^5))) + ... |
---|
270 |
(salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + (B3 .*(Ts.^3)))) + ... |
---|
271 |
(C0 .* (salini.^2)); |
---|
272 |
|
---|
273 |
o2_sol = exp(rslt); |
---|
274 |
|
---|
275 |
o2_sat = (o2_tspcoor .* 2.2414) ./ o2_sol; |
---|
276 |
|
---|
277 |
% convert ptime into datenum style... |
---|
278 |
ptime_ebd_datenum = (ptime_ebd/3600/24) + datenum(1970, 1, 1, 0, 0, 0); |
---|
279 |
|
---|
280 |
% create configuration struct... |
---|
281 |
units = struct( ... |
---|
282 |
'oxyw_oxygen', '10e-6 mol/dm3', ... |
---|
283 |
'oxyw_saturation', 'percent', ... |
---|
284 |
'oxyw_temp', 'degrees C', ... |
---|
285 |
'oxyw_dphase', 'degrees', ... |
---|
286 |
... % 'oxyw_bphase', 'degrees', ... |
---|
287 |
... % 'oxyw_rphase', 'degrees', ... |
---|
288 |
... % 'oxyw_bamp', 'mA', ... |
---|
289 |
... % 'oxyw_bpot', 'mV', ... |
---|
290 |
... % 'oxyw_ramp', 'mA', ... |
---|
291 |
... % 'oxyw_rawtemp', 'degrees C', ... |
---|
292 |
... % 'oxyw_time', 'timestamp', ... |
---|
293 |
... % 'oxyw_installed', 'bool', ... |
---|
294 |
'ptime_ebd', 'seconds since 0000-01-01T00:00', ... |
---|
295 |
'ptime_ebd_datenum', 'days since 1970-01-01T00:00', ... |
---|
296 |
'tempi', 'degrees C', ... |
---|
297 |
'salini', 'psu', ... |
---|
298 |
'depthi', 'meters', ... |
---|
299 |
'o2_tcorr', '10e-6 mol/dm3', ... |
---|
300 |
'o2_tscorr', '10e-6 mol/dm3', ... |
---|
301 |
'o2_tspcoor', '10e-6 mol/dm3', ... |
---|
302 |
'o2_sol', 'cm3/liter at 1031 hPa', ... |
---|
303 |
'o2_sat', 'percent'); |
---|
304 |
|
---|
305 |
variable_description = struct( ... |
---|
306 |
'oxyw_oxygen', 'dissolved oxygen', ... |
---|
307 |
'oxyw_saturation', 'dissolved oxygen saturation', ... |
---|
308 |
'oxyw_temp', 'water temperature', ... |
---|
309 |
'oxyw_dphase', 'phase difference', ... |
---|
310 |
... % 'oxyw_bphase', 'blue phase', ... |
---|
311 |
... % 'oxyw_rphase', 'red phase', ... |
---|
312 |
... % 'oxyw_bamp', 'blue current bias', ... |
---|
313 |
... % 'oxyw_bpot', 'blue voltage bias', ... |
---|
314 |
... % 'oxyw_ramp', 'red current bias', ... |
---|
315 |
... % 'oxyw_rawtemp', 'raw water temperature', ... |
---|
316 |
... % 'oxyw_time', 'optode timestampe', ... |
---|
317 |
... % 'oxyw_installed', 'bool', ... |
---|
318 |
'ptime_ebd', 'science computer time', ... |
---|
319 |
'ptime_ebd_datenum', 'science computer date', ... |
---|
320 |
'tempi', 'interpolated CTD water temperature', ... |
---|
321 |
'salini', 'interpoloted CTD salinity', ... |
---|
322 |
'depthi', 'interpolated CTD depth', ... |
---|
323 |
'o2_tcorr', 'temperature corrected dissolved oxygen', ... |
---|
324 |
'o2_tscorr', 'temperature and salinity corrected dissolved oxygen', ... |
---|
325 |
'o2_tspcoor', 'temperature, salinity, and pressure corrected dissolved oxygen', ... |
---|
326 |
'o2_sol', 'corrected oxygen solubility', ... |
---|
327 |
'o2_sat', 'corrected oxygen saturation'); |
---|
328 |
|
---|
329 |
correction_coefficients = struct('C0coef', C0coef, ... |
---|
330 |
'C1coef', C1coef, ... |
---|
331 |
'C2coef', C2coef, ... |
---|
332 |
'C3coef', C3coef, ... |
---|
333 |
'C4coef', C4coef); |
---|
334 |
|
---|
335 |
config = struct('glider_name', strGliderName,... |
---|
336 |
'deployment_number', strDeploymentNumber,... |
---|
337 |
'start_date', strStartDate,... |
---|
338 |
'end_date', strEndDate,... |
---|
339 |
'correction_coefficients', correction_coefficients,... |
---|
340 |
'var_descriptions', variable_description,... |
---|
341 |
'var_units', units); |
---|
342 |
|
---|
343 |
% set Level 1 data mat file name... |
---|
344 |
strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_DO_L1.mat'); |
---|
345 |
|
---|
346 |
% save flight data to mat file... |
---|
347 |
save(strMatFileName,... |
---|
348 |
'config', ... |
---|
349 |
'oxyw_oxygen', ... |
---|
350 |
'oxyw_saturation', ... |
---|
351 |
'oxyw_temp', ... |
---|
352 |
'oxyw_dphase', ... |
---|
353 |
... % 'oxyw_bphase', ... |
---|
354 |
... % 'oxyw_rphase', ... |
---|
355 |
... % 'oxyw_bamp', ... |
---|
356 |
... % 'oxyw_bpot', ... |
---|
357 |
... % 'oxyw_ramp', ... |
---|
358 |
... % 'oxyw_rawtemp', ... |
---|
359 |
... % 'oxyw_time', ... |
---|
360 |
... % 'oxyw_installed', ... |
---|
361 |
'ptime_ebd', ... |
---|
362 |
'ptime_ebd_datenum', ... |
---|
363 |
'tempi', ... |
---|
364 |
'salini', ... |
---|
365 |
'depthi', ... |
---|
366 |
'o2_tcorr', ... |
---|
367 |
'o2_tscorr', ... |
---|
368 |
'o2_tspcoor', ... |
---|
369 |
'o2_sol', ... |
---|
370 |
'o2_sat'); |
---|
371 |
|
---|
372 |
end |
---|
373 |
end |
---|