function [error] = corTsgGradient(hMainFig, PARA, dateMin, dateMax)
% Correct the TSG salinity time series with the Water sample.
% The correction is a linear interpolation of the sample/tsg difference
% between 2 dates as a function of salinity or temperature (NOT time)
% Only adviced when the sample/tsg difference shifts suddenly when
% crossing a strong salinity gradient, to avoid a correction discontinuity
%
% Input
% hMainFig ..... Handle to the main GUI
% PARA ..........Cell array
% PARA{1} contains the characters (SSP, SSJT, SSTP)
% PARA{2} contains either the cahracters (SSPS, SSJT, SSTP)
% or (SSPS_CAL, SSJT_CAL, SSTP_CAL)
% dateMin ...... the correction is applied between dateMin and date Max
% dateMax ...... the correction is applied between dateMin and date Max
%
% Output
% Error ........ 1 everything OK
% ........ -1 dateMax <= date Min
%
% 2009/12/01 G. Alory
% This method is adapted from G. Reverdin, and used when Nuka Arctica
% crosses the salinity front at the southern tip of Groenland
%
% $Id$
% Get application data
% --------------------
tsg = getappdata( hMainFig, 'tsg_data');
SAMPLE = tsg.plot.sample;
% -------------------------------------------------------------------------
% Get from the checkbox the QC code on which the correction will be applied
% -------------------------------------------------------------------------
% get list of keys from hashtable tsg.qc.hash, defined inside
% tsg_initialisation.m
% -----------------------------------------------------------
qc_list = keys(tsg.qc.hash);
% TODO: define size of keptCode
% -----------------------------
%keptCode = zeros(numel(qc_list), 1);
% iterate (loop) on each key store inside hastable
% ------------------------------------------------
keptCode = [];
nKeptCode = 0;
for key = qc_list
% get handle of checkbox
% ----------------------
hCb = findobj(hMainFig, 'tag', ['TAG_CHECK_CORRECTION_' char(key)]);
if get( hCb, 'value' )
nKeptCode = nKeptCode + 1;
keptCode(nKeptCode) = tsg.qc.hash.(key).code;
end
end
% Get PROBABLY_GOOD, PROBABLY_BAD and VALUE_CHANGED codes
% -------------------------------------------------------
PROBABLY_GOOD = tsg.qc.hash.PROBABLY_GOOD.code;
PROBABLY_BAD = tsg.qc.hash.PROBABLY_BAD.code;
VALUE_CHANGED = tsg.qc.hash.VALUE_CHANGED.code;
% Intialisation
% 01/09/2009 : intialisation to NaN for real and 0 for byte (QC)
% BE CAREFUL:
% netcdf toolbox failed with assertion when we write NaN to ncbyte variable
% -------------------------------------------------------------------------
if isempty( tsg.([PARA{1} '_ADJUSTED_ERROR']) )
tsg.([PARA{1} '_ADJUSTED']) = NaN*ones(size(tsg.(PARA{1})));
tsg.([PARA{1} '_ADJUSTED_QC']) = zeros(size(tsg.([PARA{1} '_QC'])));
tsg.([PARA{1} '_ADJUSTED_ERROR']) = NaN*ones(size(tsg.(PARA{1})));
end
if dateMax > dateMin
% Find samples within TIME_WINDOWS with Good, probably Good, QC
% -------------------------------------------------------------
ind = find( tsg.DAYD_EXT >= dateMin & tsg.DAYD_EXT <= dateMax &...
tsg.([SAMPLE '_EXT_QC']) <= PROBABLY_GOOD);
if ~isempty(ind)
% detect NaN in sample.SSPS_DIF due to bad QC code for tsg.SSPS
% -------------------------------------------------------------
ind2 = find(~isnan(tsg.EXT_DIF(ind)));
if ~isempty(ind2) && length(ind2) >= 2
% Locate front with big shift in TSG/sample difference
% -------------------------------------------------------------
% tsg_diff=diff(tsg.EXT_SMOOTH(ind(ind2)));
% corr_diff=diff(tsg.EXT_DIF(ind(ind2)));
% The correction is applied to the TSG between dateMin and dateMax using
% a linear interpolation only on measurements with keptCode Quality Codes
% ------------------------------------------------------------------------
dtTsgQCkept=[];
for icode = 1 : length( keptCode )
dtTsg = find( tsg.DAYD >= tsg.DAYD_EXT(ind(ind2(1))) & tsg.DAYD <= tsg.DAYD_EXT(ind(ind2(end))) &...
tsg.([PARA{1} '_QC']) == keptCode( icode ));
if ~isempty( dtTsg )
dtTsgQCkept=[dtTsgQCkept; dtTsg];
% Compute the correction + the error
% ----------------------------------
tsg_diff=diff(tsg.EXT_SMOOTH(ind(ind2)));
for iint=1:(length(ind2)-1)
% in each interval between 2 samples, the correction is a linear interpolation
% of the sample/tsg difference as a function of Tsg SSPS or SSTP,
% the error an interpolation of 0.5*abs(sample/tsg diff.)
dtTsg2 = find( tsg.DAYD(dtTsg) >= tsg.DAYD_EXT(ind(ind2(iint))) &...
tsg.DAYD(dtTsg) <= tsg.DAYD_EXT(ind(ind2(iint+1))) &...
tsg.([PARA{1} '_QC'])(dtTsg) == keptCode( icode ));
if ~isempty(dtTsg2)
tsg_rel=(tsg.(PARA{2})(dtTsg(dtTsg2))-tsg.EXT_SMOOTH(ind(ind2(iint))))/tsg_diff(iint);
tsg_rel(tsg_rel<0)=0;
tsg_rel(tsg_rel>1)=1;
tsg.([PARA{1} '_ADJUSTED'])(dtTsg(dtTsg2)) =...
(1-tsg_rel)*tsg.EXT_DIF(ind(ind2(iint))) + tsg_rel*tsg.EXT_DIF(ind(ind2(iint+1)));
tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsg(dtTsg2)) =...
(1-tsg_rel)*abs(tsg.EXT_DIF(ind(ind2(iint)))/2) + tsg_rel*abs(tsg.EXT_DIF(ind(ind2(iint+1)))/2);
% Compute the corrected value : original value + correction
% ---------------------------------------------------------
tsg.([PARA{1} '_ADJUSTED'])(dtTsg(dtTsg2)) =...
tsg.(PARA{2})(dtTsg(dtTsg2)) + tsg.([PARA{1} '_ADJUSTED'])(dtTsg(dtTsg2));
end
end
% Transfer the QC
% ---------------
tsg.([PARA{1} '_ADJUSTED_QC'])(dtTsg) = tsg.([PARA{1} '_QC'])(dtTsg);
end
end
% The error minimum cannot be lower than the sensor accuracy
% ----------------------------------------------------------
accuracy= tsg_accuracy(hMainFig, PARA{1},dtTsgQCkept );
error_too_small=find(abs(tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsgQCkept)) < accuracy);
if ~isempty(error_too_small)
if length(ind2) > 2
tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsgQCkept(error_too_small))=accuracy(error_too_small);
else
tsg.([PARA{1} '_ADJUSTED_ERROR'])(dtTsgQCkept(error_too_small))=-accuracy(error_too_small);
end
end
% Case with 1 sample only: we can't apply this correction
% --------------------------------------------------------------------------
else
msgbox( {'Function corTsgGradient'; ' '; ...
'At least 2 samples are needed';...
'The soft cannot make the interpolation'},...
'Correction Method', 'warn');
error = -1;
return
end
end
% Update tsg application data
% ---------------------------
setappdata( hMainFig, 'tsg_data', tsg);
% everything OK
% -------------
error = 1;
else
% DateMax <= DateMin
% ------------------
error = -1;
end