diff --git a/tsg_util/plot_Validation.m b/tsg_util/plot_Validation.m
index 0c5380f460e9334ebc7ce031984f0a45d93c2193..63ce1321d8b0f2c9ffbfaf46758ff78179d4ad20 100644
--- a/tsg_util/plot_Validation.m
+++ b/tsg_util/plot_Validation.m
@@ -19,9 +19,9 @@ switch nPlot
     
     if ~isempty( tsg.([SAMPLE '_EXT']) )
       
-      % Plot circle for WS data
+      % Plot squares for WS data
       % -----------------------
-      ind   = 1: length( tsg.([SAMPLE '_EXT_TYPE']));
+      ind   = 1: size( tsg.([SAMPLE '_EXT_TYPE']),1);
       indWS = strmatch( 'WS', tsg.([SAMPLE '_EXT_TYPE']), 'exact');
       if ~isempty(ind)
         plot_Tsg( hMainFig, hPlotAxes, 1,...
@@ -30,7 +30,7 @@ switch nPlot
           [SAMPLE '_EXT_1'],'','none','square',5);
       end
       
-      % Plot squares for CTD, ARGO, etc. data
+      % Plot circles for CTD, ARGO, etc. data
       % -------------------------------------
       indEXT = setxor(ind, indWS);
       if ~isempty(indEXT)
diff --git a/tsg_util/shipVelocity.m b/tsg_util/shipVelocity.m
index 9ff773b93deffb54f645c0f6338d57baa1662757..e6b2014bbdcf1d7e80255dcd9e1479eb83ac5af7 100644
--- a/tsg_util/shipVelocity.m
+++ b/tsg_util/shipVelocity.m
@@ -14,55 +14,75 @@ tsg = getappdata( hMainFig, 'tsg_data');
 tsgNames     = fieldnames(tsg);
 nbFieldNames = length( tsgNames );
 
+% Record number
+% -------------
+tsgLength = length(tsg.DAYD);
+
 % Test for bad velocity.
 % Suppress the bad records
 % ------------------------
 nBadVelocity   = 0;
 indBadVelocity = 1;
+indStart       = 1;
+velocity = nan*ones(tsgLength,1);
+
 while ~isempty( indBadVelocity )
 
-  % Record number
-  % -------------
-  tsgLength = length(tsg.DAYD);
-  
-  % Spherical earth distance in km - Velocity in knots
-  % --------------------------------------------------
-  range    = m_lldist(tsg.LONX,tsg.LATX);
-  ind      = size(tsg.DAYD);
-  velocity = zeros(size(ind));
-  for i = 1:tsgLength-1
-    velocity(i) = range(i) / ((tsg.DAYD(i+1)-tsg.DAYD(i)) * 24 * 1.854);
-  end
-  velocity = [velocity';0];
-
-  % Find velocity > 50 knots
-  % ------------------------
-  indBadVelocity = find( velocity > 50 );
-  
-  if ~isempty( indBadVelocity )
-    
-    % Keep only the first indice. One bad position lead to 2 bad velocities
-    % ---------------------------------------------------------------------
-    indBadVelocity = indBadVelocity(1) + 1;
-
-    % Store the number of bad Velocity
-    % --------------------------------
-    nBadVelocity = nBadVelocity + 1;
-
-    % Delete bad record
-    % -----------------
-    for i = 1 : nbFieldNames
-
-      para = char( tsgNames{i} );
-
-      % Find array of length DAYD in the TSG structure
-      % ----------------------------------------------
-      if length( tsg.(para) ) == tsgLength
-        tsg.(para)(indBadVelocity,:) = [];
-      end
+    % Spherical earth distance in km - Velocity in knots
+    % --------------------------------------------------
+    range    = m_lldist(tsg.LONX(indStart:tsgLength),tsg.LATX(indStart:tsgLength));
+    velocity(indStart:tsgLength) = [range;0] ./ ((circshift(tsg.DAYD(indStart:tsgLength),-1)...
+        -tsg.DAYD(indStart:tsgLength)) * 24 * 1.854);
+
+    % Find velocity > 50 knots
+    % ------------------------
+    indBadVelocity = find( velocity(indStart:tsgLength) > 50 );
+
+    if ~isempty( indBadVelocity )
+
+        nextBad=circshift(indBadVelocity,-1)-indBadVelocity;
+
+        if (nextBad(1)==1)
+            
+            % Case when one bad position lead to 2 bad velocities
+            % ---------------------------------------------------
+            indBadVelocity = indBadVelocity(1) + 1;
+
+            % Delete bad record
+            % -----------------
+            for i = 1 : nbFieldNames
+
+                para = char( tsgNames{i} );
+
+                % Find array of length DAYD in the TSG structure
+                % ----------------------------------------------
+                if length( tsg.(para) ) == tsgLength
+                    tsg.(para)(indStart+indBadVelocity-1,:) = [];
+                end
+
+            end
+
+            tsgLength=tsgLength-1;
+            indStart=indStart+indBadVelocity-2;
+
+        else
+            
+            % Case when sudden shift in position: apply NaN to velocity
+            % ---------------------------------------------------------
+            indBadVelocity = indBadVelocity(1);
+            velocity(indStart+indBadVelocity-1)=nan;
+            indStart=indStart+indBadVelocity;
+
+        end
+
+        % Store the number of bad Velocity
+        % --------------------------------
+        nBadVelocity = nBadVelocity + 1;
+
+
+
     end
-  end
-  
+
 end
 
 % Keep the number of data not in increasing date
@@ -72,7 +92,7 @@ tsg.report.badvelocity = nBadVelocity;
 % Keep ship velocity from positions if sog not available
 % ------------------------------------------------------
 if isempty(tsg.SPDC)
-  tsg.SPDC = velocity;
+  tsg.SPDC = velocity(1:tsgLength);
 end
 
 % Save tsg structure