
%run_tag = {'010' '011' '012' '013' '014' '015' '016' '017' '018'};
%run_tag = {'020' '021'};
%run_tag = {'030'}; modlist = { 'b2' 'b4' 'b6' 'b8'};
%run_tag = { '031'}; modlist = { 'b1' 'b2'};
%run_tag = { '033'};
%run_tag={'035'};
%run_tag = {'036' '037' '038' '039'};
%run_tag = {'038'};
%run_tag={'040'};
%run_tag={'041'};
%run_tag={'052' '054'};
%dir_root = '/share/3/output/fennelk/mch2mod/out_';
%run_tag = {'050'};
%dir_root = '/share/3/output/fennelk/mch2mod/out_';
%run_tag = {'050b'};
%dir_root = '/share/2/output/jiatang/mch2mod/out_';
%run_tag={'056' '057'};
%dir_root = '/share/5/output/fennelk/mch2mod/out_';

%run_tag={'043' '053' '070' '071'};
%dir_root = '/share/5/output/fennelk/mch2mod/out_';

%run_tag={'058'};
%dir_root = '/share/3/output/fennelk/mch2mod/out_';

run_tag = { '082' '083'};
dir_root = '/share/5/output/fennelk/mch2mod/out_';

% need to set value for ox_crit
% standard 2ml/l or 62.5 mmol/m3
ox_crit=63; % mmol/m3
fn_str = '';
% higher: 3ml/l or 93.75 mmol/m3
%ox_crit = 94;
%fn_str = '3mlpl_';


for rt=1:length(run_tag)
  run_tag{rt}

% location of oxygen data
  dirname=[dir_root run_tag{rt}];
  %filenameroot=['his_mch_' run_tag{rt}];
  filenameroot='his_mch_';

  f = roms_getfilenames(dirname,filenameroot);
  
  if exist('modlist')
    for modno =1:length(modlist)
      dirname2 = [dir_root run_tag{rt} modlist{modno}]
      f2 = roms_getfilenames(dirname2,filenameroot);
      f = [f f2];
    end
  elseif ~isempty( dir([dirname 'b/' filenameroot '*']) )
    dirname2 = [dir_root run_tag{rt} 'b']
    f2 = roms_getfilenames(dirname2,filenameroot);
    f = [f f2];
  end

  ctl = rnt_timectl(f,'ocean_time',[1970,1,0]);
  n_files = length(ctl.file);

  % load grid
  grid=grd('mch');
  % create inverted land mask
  inv_land_mask = ones(size(grid.mask_rho))-grid.mask_rho;

  % make 3D array that contains volume of each grid box
  area_mask = 1./grid.pm./grid.pn; % area in m^2 
  delta_z = diff(grid.z_w); % vertical thicknesses in m
  vol_mask = zeros(size(delta_z));
  for k=1:grid.N
    vol_mask(k,:,:) = squeeze(delta_z(k,:,:)).*area_mask; % volume in m^3
  end 

  hyp_data = []; hyparea = []; hypvol = [];
  % loop over time
  ii=1;
  for nf=1:n_files
    nc=netcdf(ctl.file{nf});
    for t = ctl.segm(nf)+1:ctl.segm(nf+1)
      %ii
      oxygen = nc{ 'oxygen'}(ctl.ind(t),:,:,:);
      % create oxygen mask (1 where hypoxic, zero elsewhere and on land!)
      ox_mask_3D = zeros(size(vol_mask));
      for k=1:grid.N
        ox = squeeze(oxygen(k,:,:));
        ind = find(ox<=ox_crit);
        ox_mask = zeros(size(ox));
        ox_mask(ind)=1;
        ox_mask = ox_mask-inv_land_mask;
        ox_mask_3D(k,:,:) = ox_mask;
      end
      data.hyp_2d_area{ii} = sum(ox_mask_3D)/grid.N;
      % calculate hypoxic area
      ox_area_mask = change(sum(ox_mask_3D),'>',0,1);
      hyparea(ii) = area_mask(:)'*ox_area_mask(:); % in m^2
      % calculate hypoxic volume
      hypvol(ii)=vol_mask(:)'*ox_mask_3D(:); % in m^3
      ii=ii+1; 
    end
    close(nc)
  end

  data.dirname = dirname;
  data.filenameroot = filenameroot;
  data.ox_crit = ox_crit;
  data.date = ctl.date;
  data.hyp_area = hyparea;
  data.hyp_vol = hypvol;
  data.run_tag = run_tag{rt};

  outputname = ['hyp_data_' fn_str run_tag{rt}];
  eval(['save ' outputname ' data']);

end




