function plot_Climatology(hMainFig, hPlotAxes)
% Function to plot mean climatology and standard deviation
%
% Input
% -----
% hTsgGUI ............ Handel to the main user interface
% hPlotAxes .......... Handels to the 3 graphic axes
%
% Output
% ------
% none
%
% $Id$

% Get data after read_Climatology because it could be change tsg.levitus.type
% -------------------------------------------------------------------------
tsg = getappdata( hMainFig, 'tsg_data' );

% get climatology depth: 0 or 10m  -> indice 1 or 2
% -------------------------------------------------
depth = tsg.preference.levitus_depth_value;

% if reading error, tsg.levitus.data is empty, no action
% ------------------------------------------------------
if isempty(tsg.levitus.data)
  return
end

% get last selected climatology structure
% ---------------------------------------
s = get(findobj('Tag', 'TAG_UIMENU_CLIMATO_MAIN'), 'Userdata');

% select time dimension for climatology from saved structure s
% ------------------------------------------------------------
time_dim = s.time;
lat=tsg.LATX;
lon=tsg.LONX;

% shift to [-180 180] longitude range if needed
% ---------------------------------------------
if ~isempty(tsg.indcross)
  lon = mod(lon+180,360)-180;
end

dayd = tsg.DAYD;

% valid parameters in climatology file are:
% -----------------------------------------
valid_parameters = {'SSTP','SSPS','DOX1'};

% Get Climatology
%           LATX(80)  = -0.5 et LATX(81)  = 0.5
%           LONX(180) = -0.5 et LONX(181) = 0.5
% ---------------------------------------------
for i = 1:2
  
  % get current parameter
  % ----------------------
  para = tsg.plot.parameter{i};
  
  % for climatology, plot SSJT as SSTP (SST)
  % ----------------------
  if strcmp( para, 'SSJT' )
    para = 'SSTP';
  end
  
  % goto next loop if parameter is not a member of climatology file
  % ----------------------
  if isempty(find(strcmp(para, valid_parameters),1))
    continue
  end
  
  % for a valid parameter, plot climatology
  % --------------------------------------
  if para
    
    % prepare the interpolation
    % ---------------------------
    latc = tsg.levitus.data.WOA_LATX;
    lonc = tsg.levitus.data.WOA_LONX;
    clim = tsg.levitus.data.(['WOA_MEAN_' para]);
    % remove sigleton dimension
    clim = squeeze(clim(time_dim,depth,:,:));
    % interpolation for 2-D gridded
    meanc = interp2(lonc,latc,clim,lon,lat);
    % same for standard deviation
    sclim = tsg.levitus.data.(['WOA_STD_' para]);
    sclim = squeeze(sclim(time_dim,depth,:,:));
    stdc = interp2(lonc,latc,sclim,lon,lat);
    
    % Plot mean salinity climatology
    % ------------------------------
    if verLessThan('matlab', '9.1.0')  % R2016b
      % Select the axes
      % ---------------
      axes( hPlotAxes(i) );
      line(dayd, meanc, ...
        'Tag', ['TAG_LINE_CLIMATO_MEAN_' para], 'Linestyle', '-', 'Color','k');
      
      % Plot with 3 standard deviation
      % ------------------------------
      line(dayd,  meanc + 3 * stdc, ...
        'Tag', ['TAG_LINE_CLIMATO_STDDEV_PLUS_' para], 'Linestyle', '-', 'Color','r');
      line(dayd,  meanc - 3 * stdc, ...
        'Tag', ['TAG_LINE_CLIMATO_STDDEV_MINUS_' para], 'Linestyle', '-', 'Color','r');
      
      % Plot marker if climatological value has no neighbour
      % --------------------------------------------------
      indnon=find(isfinite([meanc',nan]) & isnan(circshift([meanc',nan],[0 1]))...
        & isnan(circshift([meanc',nan],[0 -1])));
      if ~isempty(indnon)
        line(dayd(indnon), meanc(indnon), ...
          'Tag', ['TAG_LINE_CLIMATO_MEAN_' para], 'Linestyle', 'none', 'Color','k','Marker','d');
        line(dayd(indnon),  meanc(indnon) + 3 * stdc(indnon), ...
          'Tag', ['TAG_LINE_CLIMATO_STDDEV_PLUS_' para], 'Linestyle', 'none', 'Color','r','Marker','d');
        line(dayd(indnon),  meanc(indnon) - 3 * stdc(indnon), ...
          'Tag', ['TAG_LINE_CLIMATO_STDDEV_MINUS_' para], 'Linestyle', 'none', 'Color','r','Marker','d');
      end
    else
      line(hPlotAxes(i), dayd, meanc, ...
        'Tag', ['TAG_LINE_CLIMATO_MEAN_' para], 'Linestyle', '-', 'Color','k');
      
      % Plot with 3 standard deviation
      % ------------------------------
      line(hPlotAxes(i), dayd,  meanc + 3 * stdc, ...
        'Tag', ['TAG_LINE_CLIMATO_STDDEV_PLUS_' para], 'Linestyle', '-', 'Color','r');
      line(hPlotAxes(i), dayd,  meanc - 3 * stdc, ...
        'Tag', ['TAG_LINE_CLIMATO_STDDEV_MINUS_' para], 'Linestyle', '-', 'Color','r');
      
      % Plot marker if climatological value has no neighbour
      % --------------------------------------------------
      indnon=find(isfinite([meanc',nan]) & isnan(circshift([meanc',nan],[0 1]))...
        & isnan(circshift([meanc',nan],[0 -1])));
      if ~isempty(indnon)
        line(hPlotAxes(i), dayd(indnon), meanc(indnon), ...
          'Tag', ['TAG_LINE_CLIMATO_MEAN_' para], 'Linestyle', 'none', 'Color','k','Marker','d');
        line(hPlotAxes(i), dayd(indnon),  meanc(indnon) + 3 * stdc(indnon), ...
          'Tag', ['TAG_LINE_CLIMATO_STDDEV_PLUS_' para], 'Linestyle', 'none', 'Color','r','Marker','d');
        line(hPlotAxes(i), dayd(indnon),  meanc(indnon) - 3 * stdc(indnon), ...
          'Tag', ['TAG_LINE_CLIMATO_STDDEV_MINUS_' para], 'Linestyle', 'none', 'Color','r','Marker','d');
      end
    end
    
  end
  
end

% save tsg structure
% ------------------
setappdata( hMainFig, 'tsg_data', tsg );

% Update the map if already displayed
% -----------------------------------
hMapFig = get(get(hPlotAxes(4),'parent'),'parent');

if strcmp( get(hMapFig,'visible'), 'on') == 1
  erase_Line( hPlotAxes, 4 );
  plot_map( hMainFig, hPlotAxes);
end