% plot hypoxic area over time for SURA period
% [nice plots for the manuscript]

% chose among the different options

  tag = {               '058'     '053'     '043'}; 
  leg = {'observed' 'A20clim' 'B20clim' 'C20clim'}; % for legend
  lc =  {                 'b'       'g'       'r'}; % line color
  fn = 'comp_hyp_area_SURA_clim20_runs'; % output file name

% tag = {               '058'     '053'     '071'     '070'}; 
% leg = {'observed' 'A20clim' 'B20clim' 'A30clim' 'B30clim'}; % for legend
% lc =  {                 'b'       'm'       'g'       'c'}; % line color
% fn = 'comp_hyp_area_SURA_vert_res'; % output file name
% 
% tag = {               '071'    '081'    '083'}; 
% leg = {'observed' 'A30clim' 'A30HYC' 'A30IAS'}; % for legend
% lc =  {                 'b'       'm'       'g'}; % line color
% fn = 'comp_hyp_area_SURA_IR_bc'; % output file name

% tag = {               '070'    '080'    '082'}; 
% leg = {'observed' 'B30clim' 'B30HYC' 'B30IAS'}; % for legend
% lc =  {                 'b'       'm'       'g'}; % line color
% fn = 'comp_hyp_area_SURA_HD_bc'; % output file name


% tag = {               '058'     '053' }; 
% leg = {'observed' 'A20clim' 'B20clim' }; % for legend
% lc =  {                 'b'       'm' }; % line color
% fn = 'comp_hyp_area_SURA_w_envelope'; % output file name

w_envelope = 0; % true or false
env_option = {'IR' 'HD'};


% load Paul's uncertainty envelopes
%load rpc20p_hypoxicareaquantiles_katja
% variables in the files:
%   ir_areaquantiles       5x2189            87560  double              
%   ir_dates               1x2189            17512  double              
%   ir_largeenv            2x2189            35024  double              
%   ir_smallenv            2x2189            35024  double              
%   q                      1x5                  40  double              
%   sod_areaquantiles      5x1775            71000  double              
%   sod_dates              1x1775            14200  double              
%   sod_largeenv           2x1775            28400  double              
%   sod_smallenv           2x1775            28400  double      
% From Paul:
% Includes quantiles for both the SOD and IR scenario (I used SOD for my plots). I
% have also computed the resulting uncertainty envelopes which are really 
% just the respective quantiles minus the median.


print_figs=0;

% load data
for i=1:length(tag)
  eval(['load ./catZ/hyp_data_' tag{i}])
  x{i} = data;
  x{i}.datenum = datenum(x{i}.date(:,1:3));
end 
if w_envelope
  % Note Paul's data is in km^2; convert to m^2 by multiplying by 1e+6
  for i=1:length(tag)
    if env_option{i} == 'IR' 
      [C,inA,inB] = intersect(ir_dates,x{i}.datenum);
      x{i}.twentyfivepct = x{i}.hyp_area(inB)+ir_largeenv(1,inA)*1e+6;
      x{i}.seventyfivepct = x{i}.hyp_area(inB)+ir_largeenv(2,inA)*1e+6;
      x{i}.pcttime = ir_dates(inA);
    else
      [C,inA,inB] = intersect(sod_dates,x{i}.datenum);
      x{i}.twentyfivepct = x{i}.hyp_area(inB)+sod_largeenv(1,inA)*1e+6;
      x{i}.seventyfivepct = x{i}.hyp_area(inB)+sod_largeenv(2,inA)*1e+6;
      x{i}.pcttime = x{i}.datenum(inB);
    end
  end
end
% Note: to convert from m2 to km2, multiply by 1e-6
%      (1km2 = 10^3*2 = 10^6 m2)


% my own calculations of hypoxic area
obs.area = [ 11014  10124   9428  15299];
obs.yr   = [  2004   2005   2006   2007];
obs.time = datenum([obs.yr; 7*ones(1,4); 25*ones(1,4)]');

% ticks location for time axis
tickdates_SURA = [2004 01 01; 2004 07 01; 2005 01 01; 2005 07 01; 2006 01 01;
                  2006 07 01; 2007 01 01; 2007 07 01; 2008 01 01];            
        
figure
plot(obs.time,obs.area/1000,'k.','MarkerSize',16); hold on
for i=1:numel(tag)
  plot(x{i}.datenum,x{i}.hyp_area*1e-9,lc{i},'LineWidth',1);
end
for i=1:numel(tag)
  if w_envelope
    hp=plotarea(x{i}.pcttime,x{i}.twentyfivepct*1e-9,x{i}.seventyfivepct*1e-9,lc{i}); 
    set(hp,'FaceAlpha',0.3,'EdgeAlpha',0.3)
  end
end
legend(leg)
xlim([datenum(2004,1,1) datenum(2008,0,0)])
set(gca,'XTick',datenum(tickdates_SURA))
datetick('x',12,'keepticks')
set(gca,'FontSize',13)
yl = get(gca,'Ylim'); set(gca,'Ylim',[0 yl(2)])
ylabel('Hypoxic area (10^3 km^2)')
    
if print_figs
  set(gcf,'PaperPosition',[0.25 2.5 8 3])
  if w_envelope
    print(gcf,'-dpng','-r300',fn)
  else
    print(gcf,'-depsc',fn)     
  end
end



