
load oxygenplot_alldata



z = [0 , 10, 20, 30, 50, 75, 100, 125, 150, 200, 500, 1000];
level = {'_z0m' '_z10m' '_z20m' '_z30m' '_z50m' '_z75m' '_z100m' '_z125m' '_z150m' '_z200m' '_z500m' '_z1000m'};
for i = 1 : 12
  levelt = level{i};
   oxx=ox;
   latt=lat;
   lonn=lon;
   monthss=months;
   
   
  bad1=find(depth1 <z(i)-2 | depth1 >z(i)+2 | depth1<0);
  
  oxx(bad1)=[];
  latt(bad1)=[];
  lonn(bad1)=[];
  monthss(bad1)=[];
  
  
  bad2 =find(depth2 ~= z(i));
  oxx(bad2)=[];
  latt(bad2)=[];
  lonn(bad2)=[];
  monthss(bad2)=[];
  
  
delta = 1;
     latmin = 30;
     latmax = 50;
     lonmin = -80;
     lonmax = -60;
     imax = ceil((lonmax-lonmin)/delta);
     jmax = ceil((latmax-latmin)/delta);
     bimax = imax;
     bjmax = jmax;

     lattrans = (latt - latmin)/(latmax - latmin);
     latind = ceil(jmax*lattrans);
     lontrans = (lonn - lonmin)/(lonmax - lonmin);
     lonind = ceil(imax*lontrans);
  
     bin = cell(12, bimax, bjmax);

 oxmonthmean = nan(12, bimax, bjmax);
oxmonthstd = nan(12, bimax, bjmax);
oxmonthcount = nan(12, bimax, bjmax); 
oxseasonmean = nan (4,bimax, bjmax);
oxseasonstd = nan (4,bimax, bjmax);
oxseasoncount = nan (4,bimax, bjmax); 
oxseasoncountnew = nan (4, bimax, bjmax);
    for m = 1:12
       
 
        for bi = 1:bimax;
            for bj = 1:bjmax;
                bin{m, bi, bj} = find(latind == bi & lonind == bj & monthss == m & oxx >= 0 & oxx< 600);
                %             fprintf('%d %d numel: %d\n', bi, bj, numel(bin{m, bi, bj}))
               
                
                oxmonthmean(m, bi, bj) = mean(oxx(bin{m, bi, bj}));
                oxmonthstd(m, bi, bj) = std(oxx(bin{m, bi, bj}));
                oxmonthcount(m, bi, bj) = length(bin{m, bi, bj});
               
            end
        end
    end
    

  for bi = 1 : bimax;
            for bj = 1 : bjmax
                oxseasonmean(1, bi, bj )=mean([oxx(bin{3,bi,bj});oxx(bin{4,bi,bj});oxx(bin{5,bi,bj})]);
                oxseasonstd(1,bi,bj)=std([oxx(bin{3,bi,bj});oxx(bin{4,bi,bj});oxx(bin{5,bi,bj})]);
                oxseasoncount(1,bi,bj)=oxmonthcount(3,bi,bj)+oxmonthcount(4,bi,bj)+oxmonthcount(5,bi,bj);
                oxseasoncountnew(1, bi,bj)= change(oxseasoncount(1, bi, bj), '==', 0, nan );
                oxseasonmean(2, bi, bj )=mean([oxx(bin{6,bi,bj});oxx(bin{7,bi,bj});oxx(bin{8,bi,bj})]);
                oxseasonstd(2,bi,bj)=std([oxx(bin{6,bi,bj});oxx(bin{7,bi,bj});oxx(bin{8,bi,bj})]);
                oxseasoncount(2,bi,bj)=oxmonthcount(6,bi,bj)+oxmonthcount(7,bi,bj)+oxmonthcount(8,bi,bj);
                 oxseasoncountnew(2, bi,bj)= change(oxseasoncount(2, bi, bj), '==', 0, nan );
                oxseasonmean(3, bi, bj)=mean([oxx(bin{9,bi,bj});oxx(bin{10,bi,bj});oxx(bin{11,bi,bj})]);
                oxseasonstd(3,bi,bj)=std([oxx(bin{9,bi,bj});oxx(bin{10,bi,bj});oxx(bin{11,bi,bj})]);
                oxseasoncount(3,bi,bj)=oxmonthcount(9,bi,bj)+oxmonthcount(10,bi,bj)+oxmonthcount(11,bi,bj);
                 oxseasoncountnew(3, bi,bj)= change(oxseasoncount(3, bi, bj), '==', 0, nan );
                oxseasonmean(4, bi, bj)=mean([oxx(bin{1,bi,bj});oxx(bin{2,bi,bj});oxx(bin{12,bi,bj})]);
                oxseasonstd(4,bi,bj)=std([oxx(bin{1,bi,bj});oxx(bin{2,bi,bj});oxx(bin{12,bi,bj})]);
                oxseasoncount(4,bi,bj)=oxmonthcount(1,bi,bj)+oxmonthcount(2,bi,bj)+oxmonthcount(12,bi,bj);
                 oxseasoncountnew(4, bi,bj)= change(oxseasoncount(4, bi, bj), '==', 0, nan );
            end
  end   

season = {'Spring' 'Summer' 'Fall' 'Winter'};

longrid = linspace(lonmin+delta/2,lonmax-delta/2,20); %names variable to be used for plotting
latgrid = linspace(latmin+delta/2,latmax-delta/2,20); %names variable to be used for plotting
 
  for s = 1:4  
     figure
    %subplot(1,4,s), 
    pcolorjw(longrid, latgrid, squeeze(oxseasonstd(s,:,:)))
    
    %subplots surf. mean O2 values of each month
    shading flat
    pnc
    colorbar
    caxis([0,200])
    han = colorbar;
    set(han , 'pos',[0.3 0.55 0.02 0.35])
    
    title(['Seasonal std value [O_2] (mmol/m^3) in ' ' ', season{s}])
    xlabel('Longitude (^oE)')
    ylabel('Latitude (^oN)')
    base = ['oxygenplot/seasonal/standard deviation/1degree/oxyg_1deg_' '' season{s} ];
    print(gcf, '-dpng', [base levelt]);
  end
end 
 