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$


% Read surface climatology (annual, seasonal or monthly)
% ------------------------------------------------------
read_Climatology(hMainFig);

% 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 form saved structure s
% ------------------------------------------------------------
time_dim = s.time;

% round positive latitude and Longitude toward zero
% -------------------------------------------------
ind = find(tsg.LATX > 0);
lat(ind) = fix(tsg.LATX(ind)) + 0.5;

ind = find(tsg.LONX > 0);
lon(ind) = fix(tsg.LONX(ind)) + 0.5;

% rounds negative latitude and Longitudeto the nearest lowest integers
% ---------------------------------------------------------------------
ind = find(tsg.LATX <= 0);
lat(ind) = floor(tsg.LATX(ind)) + 0.5;

ind = find(tsg.LONX <= 0);
lon(ind) = floor(tsg.LONX(ind)) + 0.5;

% shift to [-180 180] longitude range if needed
% ---------------------------------------------
if ~isempty(tsg.indcross)
    lon=mod(lon+180,360)-180;
end
    
% Calculates differences between adjacent elements of X.
% 0 if adajacent latitude or longitude are equal
% - 1 or -1 otherwise
% ------------------------------------------------------------
lat_diff = [diff( lat )'; 0];
lon_diff = [diff( lon )'; 0];

% Select latitude and longitude
% -----------------------------
ind  = find(abs(lat_diff) == 1 | abs(lon_diff == 1));
lat2 = lat( ind );
lon2 = lon( ind );
dayd = tsg.DAYD( ind );

% Get Climatology
%           LATX(80)  = -0.5 et LATX(81)  = 0.5
%           LONX(180) = -0.5 et LONX(181) = 0.5
% ---------------------------------------------
for i = 1:2
  
  para = tsg.plot.parameter{i};
  if strcmp( para, 'SSJT' )
    para = 'SSTP';
  end

  mean = zeros(size(ind));
  std  = zeros(size(ind));
  for ii=1:length(ind)
    ilat     = find(tsg.levitus.data.WOA_LATX == lat2(ii));
    ilon     = find(tsg.levitus.data.WOA_LONX == lon2(ii));
    mean(ii) = tsg.levitus.data.(['WOA_MEAN_' para])(time_dim,depth,ilat,ilon);
    std(ii)  = tsg.levitus.data.(['WOA_STD_'  para])(time_dim,depth,ilat,ilon);
  end

  % Select the axes
  % ---------------
  axes( hPlotAxes(i) );

  % Plot mean salinity climatology
  % ------------------------------
  line(dayd, mean, ...
    'Tag', ['TAG_LINE_CLIMATO_MEAN_' para], 'Linestyle', '-', 'Color','k');

  % Plot with 3 standard deviation
  % ------------------------------
  line(dayd,  mean + 3 * std, ...
    'Tag', ['TAG_LINE_CLIMATO_STDDEV_PLUS_' para], 'Linestyle', '-', 'Color','r');
  line(dayd,  mean - 3 * std, ...
    'Tag', ['TAG_LINE_CLIMATO_STDDEV_MINUS_' para], 'Linestyle', '-', 'Color','r');
end


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

end
