% % correctThermalLag.m % % This function receives as input parameter a CTD profile, and applies a % thermal lag correction on it. % % The correction applied uses the recursive scheme described in % Morison, J., R. Andersen, N. Larson, E. D'Asaro, and T. Boyd, 1994: % The Correction for Thermal-Lag Effects in Sea-Bird CTD Data. % Journal of Atmospheric and Oceanic Technology, vol. 11, pages 1151�1164. % % Syntax: % correctedProfileData = correctThermalLag(basicProfileData) % correctedProfileData = correctThermalLag(basicProfileData, correctionParams) % [correctedProfileData, correctionParams] = correctThermalLag(basicProfileData) % % Inputs: % basicProfileData - A profile structure* % correctionParams - The set of parameters to be used in the correction* % gliderVelocity - Glider velocity from the dbd files (m_speed) % % Outputs: % correctedProfileData - A profile structure* % correctionParams - The set of parameters used in the correction* % % Profile structure: % A struct that contains several fields, all of them column vectors with the same length: % - ptime: Present time instant at which this row was collected % - depth: Depth (pressure in decibars) measured by the CTD % - temp: Temperature measured by the CTD % - cond: Conductivity measured by the CTD % - pitch: Pitch angle of the glider (optional) % The output profile has the same information of the input plus two fields with the % corrected profile properties: % - condOutCell: corrected conductivity, removing the effects of the % temperature difference between the outer and inner parts of the % conductivity cell. % - tempInCell: corrected temperature, that is, the temperature of % the water mass lying inside the conductivity cell. % % From this information, the user can choose which one of the two % corrections to use in order to compute salinity: % - Combine 'condOutCell' with 'temp' (Expected values outside of the % conductivity cell). % - Combine 'cond' with 'tempInCell' (Expected values inside of the % conductivity cell). % % Correction parameters: % A vector of four elements, consisting of alpha_offset, alpha_slope, % tau_offset and tau_slope. These parameters are used to compute alpha and tau, % the amplitude and time constant respectively, which are inversely % proportional to the flow speed. % % Example: % correctedProfileData = correctThermalLag(basicProfileData) corrects the % profile information contained in the input 'basicProfileData', % using Morison parameters % % correctedProfileData = correctThermalLag(basicProfileData, correctionParams) % corrects the profile information contained in 'basicProfileData', % using correctionParams as the parameters to be used % during the correction. % % [correctedProfileData, correctionParams] = correctThermalLag(basicProfileData) % corrects the profile information contained in 'basicProfileData', and % provides in the output the correction parameters used for the correction. % % Other m-files required: none % Subfunctions: none % MAT-files required: none % % See also: ADJUSTTHERMALLAGPARAMS % % Author: Bartolome Garau % Work address: Parc Bit, Naorte, Bloc A 2ºp. pta. 3; Palma de Mallorca SPAIN. E-07121 % Author e-mail: tgarau@socib.es % Website: http://www.socib.es % Created: 17-Feb-2011 % %////////////////////////////////////////////////////////////////////////// function [correctedProfileData, varargout] = correctThermalLag(basicProfileData, varargin) % Extract the information contained in the profile time = basicProfileData.ptime; depth = basicProfileData.depth; temp = basicProfileData.temp; cond = basicProfileData.cond; if isfield(basicProfileData, 'pitch'), pitch = basicProfileData.pitch * pi / 180; else pitch = (26 * pi / 180) * ones(size(depth)); end; % Copy the original fields into the output struct correctedProfileData = basicProfileData; if numel(time) <= 1 correctedProfileData.condOutCell = correctedProfileData.cond; correctedProfileData.tempInCell = correctedProfileData.temp; return; end; if nargin == 1, % Only basicProfileData is provided as input % These values are proposed in Morison94 alpha_offset = 0.0135; alpha_slope = 0.0264; tau_offset = 7.1499; tau_slope = 2.7858; elseif nargin == 2, % offset and slope for alpha and tau are also given paramVector = varargin{1}; alpha_offset = paramVector(1); alpha_slope = paramVector(2); tau_offset = paramVector(3); tau_slope = paramVector(4); % ADDED 6/13/2012 (William Stark, whs@unc.edu): Allow for glider velocity vector % to be passed in. elseif nargin == 3, % Glider velocity vector is also provided % copy in offset and slope for alpha and tau as before... paramVector = varargin{1}; alpha_offset = paramVector(1); alpha_slope = paramVector(2); tau_offset = paramVector(3); tau_slope = paramVector(4); % copy in glider velocity... gliderVelocity = varargin{2}; else disp('Incorrect number of parameters. Type ''help correctThermalLag'' for more information'); return; end; % Some initial precomputations... deltaTime = abs(diff(time)); deltaTemp = diff(temp) ; samplingFreq = 1 ./ deltaTime ; % Calculate the surge speed from the depth rate and pitch deltaDepth = abs(diff(depth)); % does not matter if downcast or upcast depthRate = deltaDepth ./ deltaTime; pitchSinus = sin(pitch); surgeSpeed = depthRate ./ pitchSinus(1:end-1); % The relative coefficient between the flow speed inside and outside % of the conductivity cell. This is still uncertain (ask Gordon for origin). % Here are three choices for first three orders polynomial. speedFactorPols = [0.00, 0.00, 0.40; % 0th order degree 0.00, 0.03, 0.45; % 1st order degree 1.58, 1.15, 0.70]; % 2nd order degree selectedDegree = 1; % First order approximation, second row of the matrix speedFactor = polyval(speedFactorPols(selectedDegree+1, :), surgeSpeed); % ADDED 6/13/2012 (William Stark, whs@unc.edu): if/else statement added % so that if a glider velocity was passed in it is used as the flow % speed... otherwise the flow speed is calculated using the surge speed % and speed factor as before. if (exist('gliderVelocity')==1), % if glider velocity was passed in % use gliderVelocity as the flow speed flowSpeed = gliderVelocity; % calculation of flow speed using surge speed and speed factor % results in a flow speed vector that is one element shorter. For % this reason we must trim the length of this 'passed-in' flow % speed vector by one... flowSpeed(end) = []; else; % calculate the flow speed using surge speed and speed factor... % ADDED 6/6/12 (William Stark, whs@unc.edu): set eps=0.01 TO AVOID BLOWUPS WHEN % CALLING THIS FUNCTION DIRECTLY... %flowSpeed = abs(speedFactor .* surgeSpeed) + eps; % Avoid division by zero flowSpeed = abs(speedFactor .* surgeSpeed) + 0.01; % Avoid division by zero end; % The alpha and tau parameters, as suggested in the reference paper, % depend on the flow with tne next expressions alpha = alpha_offset + alpha_slope ./ flowSpeed ; tau = tau_offset + tau_slope ./ sqrt(flowSpeed); % Relation between a and b coefficients with respect to alpha and tau coefa = 4 .* samplingFreq .* alpha .* tau ./ (1 + 4 .* samplingFreq .* tau); coefb = 1 - 2 .* coefa ./ alpha; % Sensitivity of conductivity with respect to temperature, % approximation suggested by SeaBird: SBE Data Processing User�s Manual % at Section 6: Data Processing Modules, Cell Thermal Mass % Software Release 7.16a and later. Date: 01/18/08 % dCdT = 0.1 .* (1 + 0.006 .* (temp - 20)); dCdT = 0.088 + 0.0006 * temp; % Recursive processing of the corrections condCorrection = zeros(size(cond)); tempCorrection = zeros(size(temp)); for depthLevel = 1:length(depth)-1, % Compute corrections for next depth level condCorrection(depthLevel+1) = ... - coefb(depthLevel) .* condCorrection(depthLevel) + ... coefa(depthLevel) .* dCdT(depthLevel) .* deltaTemp(depthLevel); tempCorrection(depthLevel+1) = ... - coefb(depthLevel) .* tempCorrection(depthLevel) + ... coefa(depthLevel) .* deltaTemp(depthLevel); end; % Apply corrections and save them as fields in the output struct correctedProfileData.condOutCell = cond + condCorrection; correctedProfileData.tempInCell = temp - tempCorrection; if nargout == 2, % If required, output alpha and tau offsets and slopes varargout{1} = [alpha_offset, alpha_slope, tau_offset, tau_slope]; end; if nargout == 3, % If required, output also alpha and tau parameters varargout{2} = [alpha, tau]; end; clear gliderVelocity; return;