% LAYER_LOOK 
%
% Look time series at one location in each sediment layer
%
% Benedicte Ferre - USGS
% 01/16/07
%

%% Change by the user
ncf = 'D:\bf\STRATAFORM\ROMS\roms_sed_crs_branch\MyApp\ocean_his.nc';

nmud=1;
nsand=1;
nsed = nmud+nsand; % how do I figure out how many there are?
d50 = [4 63];     % [4 22] um
tcr = [0.05 0.1]; % [0.05 0.1] N/m2
ws = [0.1 2];     % [0.1 2] mm/s
i = 2;
j = 5;

%% write on the screen
for cl = 1 : size(d50)
    sprintf('Class %d',cl)
    sprintf('Size %d %s',d50(cl),'\mum')
    sprintf('ws = %g %s',ws(cl),'mm/s')
    sprintf('tcr = %g %s',tcr(cl),'N/m2')
end
nco = netcdf(ncf)

%%
rhos =  2650.*ones(1,nsed);
t=nco{'ocean_time'}(:);
nt = length(t);
dstart = nco{'dstart'}(:);
dt = nco{'dt'}(:);
nbeddim = nco('Nbed'); % number of bed sediment levels
nbed = length(nbeddim);
nzdim = nco('N');      % number of water-column levels
nz =length(nzdim);
nxdim = nco('xi_rho');
nx = length(nxdim);
nydim = nco('eta_rho');
ny = length(nydim);

mud = zeros(nmud,nt,nz,ny,nx);
mudmass = zeros(nmud,nt,nbed,ny,nx);
mudfrac = zeros(nmud,nt,nbed,ny,nx);
sand = zeros(nsand,nt,nz,ny,nx);
sandmass = zeros(nsand,nt,nbed,ny,nx);
sandfrac = zeros(nsand,nt,nbed,ny,nx);
frac = zeros(nsed,nt,nbed,ny,nx);
sedmass = zeros(nsed,nt,nbed,ny,nx);
h=nco{'h'}(:);
ero_flux=nco{'ero_flux'}(:);

if nmud~=0, 
    for ns = 1:nmud,
       ts = sprintf('mud(%d,:,:,:,:) = nco{''mud_%02d''}(:);',ns,ns);
       eval(ts);
       ts = sprintf('mudmass(%d,:,:,:,:) = nco{''mudmass_%02d''}(:);',ns,ns);
       eval(ts);
       ts = sprintf('mudfrac(%d,:,:,:,:) = nco{''mudfrac_%02d''}(:);',ns,ns);
       eval(ts);
    end
end
if nsand~=0, 
    for ns = 1:nsand,
       ts = sprintf('sand(%d,:,:,:,:) = nco{''sand_%02d''}(:);',ns,ns);
       eval(ts);
       ts = sprintf('sandmass(%d,:,:,:,:) = nco{''sandmass_%02d''}(:);',ns,ns);
       eval(ts);
       ts = sprintf('sandfrac(%d,:,:,:,:) = nco{''sandfrac_%02d''}(:);',ns,ns);
       eval(ts);
    end
end

%%
bed_thickness = nco{'bed_thickness'}(:);
bed_porosity = nco{'bed_porosity'}(:);
bed_age = nco{'bed_age'}(:);
active_layer = nco{'active_layer_thickness'}(:);
erosion_stress = nco{'erosion_stress'}(:);
bed_biodiff = nco{'bed_biodiff'}(:);
bed_tau_crit = nco{'bed_tau_crit'}(:);
bustr = nco{'bustr'}(:);
bvstr = nco{'bvstr'}(:);
bustr2=bustr(:,1:end-1,:).^2;
bvstr2=bvstr(:,:,1:end-1).^2;
taub=sqrt(bustr2+bvstr2);
bustrcwmax = nco{'bustrcwmax'}(:);
bvstrcwmax = nco{'bvstrcwmax'}(:);
taumax=sqrt(bustrcwmax.^2+bvstrcwmax.^2);

zr = zeros(nt,nz,ny,nx);
zw = zeros(nt,nz+1,ny,nx);
dz = zeros(nt,nz,ny,nx);
for n=1:nt,
   [z_r,z_w]=get_z(n,nco);
   zr(n,:,:,:)=z_r;
   zw(n,:,:,:)=z_w;
   dz(n,:,:,:)=diff(z_w);
end

%%
susmass=zeros(nsed,nt,ny,nx);
susmass_z=zeros(nsed,nt,nz,ny,nx);

if nmud~=0,
    for n=1:nmud
        s = squeeze(mud(n,:,:,:,:)).*dz;
        susmass(n,:,:,:)=squeeze( sum(s,2) );
        susmass_z(n,:,:,:,:)= mud(n,:,:,:,:);
    end
end
s=[];
if nsand~=0,
    for n=nmud+1:nmud+nsand
        s = squeeze(sand(n-nmud,:,:,:,:)).*dz;
        susmass(n,:,:,:)=squeeze( sum(s,2) );
        susmass_z(n,:,:,:,:)= sand(n-nmud,:,:,:,:);
    end
end
taub = ( abs(squeeze(bustr(:,i,j-1))));  % kludge alert!!
% taub=squeeze(taub(:,i,j));
susmasstot = squeeze(sum(susmass,1));    % total mass in suspension
susmasstot_z = squeeze(sum(susmass_z,1));    % total mass in suspension per level

mudmass_z = sum(mudmass,3);    % total mass of each mud class in bed
sandmass_z= sum(sandmass,3);   % total mass of each sand class in bed
sedmass_z = zeros(nsed,nt,1,ny,nx);
if nmud~=0,
    for n=1:nmud
        sedmass(n,:,:,:,:)=mudmass(n,:,:,:,:);
        sedmass_z(n,:,1,:,:)=mudmass_z(n,:,:,:,:);
        frac(n,:,:,:,:)=mudfrac(n,:,:,:,:);
    end
end
if nsand~=0,
    for n=nmud+1:nmud+nsand
        sedmass(n,:,:,:,:)=sandmass(n-nmud,:,:,:,:);
        sedmass_z(n,:,1,:,:)=sandmass_z(n-nmud,:,:,:,:);
        frac(n,:,:,:,:)=sandfrac(n-nmud,:,:,:,:);
    end
end
sedmass_z=squeeze(sedmass_z);

zsed = cumsum(bed_thickness,2);           % depth in sediment
bmz =  cumsum(squeeze(sum(mudmass,1))+squeeze(sum(sandmass,1)),2); % mass depth in sediment

% this DOES equal mudmass+sandmass
mass_calc = zeros(nsed,nt,nbed,ny,nx);
for n=1:nsed
mass_calc(n,:,:,:,:) = bed_thickness.*(1-bed_porosity).*rhos(n).*squeeze(frac(n,:,:,:,:));
% mass_calc(n,:,:,:,:) = squeeze(bed_thickness).*(1-squeeze(bed_porosity)).*rhos(n).*squeeze(frac(n,:,:,:,:));
end

sedmass_ztot = squeeze(sum(sedmass_z,1));
% all = sedmass_ztot + susmasstot;         % this should not change if mass is conserved;
%%
% figure(2)
% clf
% plot(t./(24*3600),squeeze(bed_tau_crit(:,1,j,1)))

%%
if(0)
   
figure; clf

ff = zeros(2*nbed,1);
zz = zeros(2*nbed,1);
for n=1:nt
   ff(1:2:end-1,1)=squeeze(frac(1,n,:,i,j))';
   ff(2:2:end,1)=squeeze(frac(1,n,:,i,j))';
   zz(3:2:end-1)=squeeze(bmz(n,1:end-1,i,j));
   zz(2:2:end)=  squeeze(bmz(n,1:end,i,j));
   plot(ff,zz);
   %hold on
   set(gca,'ydir','reverse')
   axis([ 0 1.1 0 20])
   pause
end

end
%%
tdays = t./(24*3600);
figure('Color',[1 1 1]);
clf
scrsz = get(0,'ScreenSize');
% load('D:\bf\STRATAFORM\ROMS\Sandford_sed\MyMap','my_map')
my_map = [0 0 1;    % Blue
          1 0 0;    % Red
          1 1 0;    % Yellow
          0 1 0];   % Green
 
subplot(815)
harea = area(tdays,squeeze(susmass(:,:,i,j)'));
for ind = 1:length(harea)
    set(harea(ind),'FaceColor',my_map(ind,:))
end
set(gca,'xlim',[tdays(1) tdays(end)],'XTickLabel',[]);
ylabel('Sus. Mass (kg/m^2)')

tsM=sprintf('Mud %dum ',d50(1:nmud));
tsS=sprintf('Sand %dum ',d50(nmud+1:nsed));
legend({tsM,tsS},'Position',[0.2893 0.9509 0.433 0.03666],...
  'Orientation','horizontal');
legend('boxoff')

subplot(816)
plot(tdays,taumax(:,i,j),'-r','linewidth',2);
% plot(tdays,taub,'-r','linewidth',2);
set(gca,'xlim',[tdays(1) tdays(end)],'XTickLabel',[]);
ylabel('\tau_b (Pa)')

subplot(817)
plot(tdays,1e6*squeeze(active_layer(:,i,j)),'-b','linewidth',2);
set(gca,'xlim',[tdays(1) tdays(end)],'XTickLabel',[]);
ylabel('AL Thick (\mum)')

subplot(818)
% area(tdays,squeeze(sedmass_z(:,:,i,j)'))
plot(tdays,squeeze(ero_flux(:,i,j)),'-g','linewidth',2);
set(gca,'xlim',[tdays(1) tdays(end)]);
% ylabel('Bed Mass (kg/m^2)')
ylabel('Ero flux (kg/m^2)')
xlabel('Time (days)')

%
for n=1:nt
   subplot(241)
   hbar1 = barh( squeeze(frac(:,n,:,i,j))' ,1,'stacked');
   for ind = 1:length(hbar1)
       set(hbar1(ind),'FaceColor',my_map(ind,:))
   end
   set(gca,'ydir','reverse')
   set(gca,'xlim',[0 1])
   ylabel('Layer Number')
   xlabel('Fraction')
   
   subplot(242)
   hbar2 = barh( squeeze(sedmass(:,n,:,i,j))',1,'stacked');
   for ind = 1:length(hbar2)
       set(hbar2(ind),'FaceColor',my_map(ind,:))
   end
   set(gca,'ydir','reverse')
   xlabel('Mass (kg / m^2 )')
   maxi=ceil(max(max(max(max(max(sum(sedmass(:,:,:,i,j))))))));
   axis([0 maxi 0 12])
   ts = sprintf('Total = %g kg/m^2',bmz(n,end,i,j));
   text(.1,11,ts)
%    set(gcf,'Colormap',my_map)
   shg
   
   subplot(243)
   barh( squeeze(1000*bed_thickness(n,:,i,j))',1,'Facecolor',[0.5 0.5 0.5])
   set(gca,'ydir','reverse')
   xlabel('Thickness (mm)')
   maxi=ceil(max(max(max(1000*bed_thickness(:,:,i,j)))));
   axis([0 maxi 0 12])
   ts = sprintf('Total = %g mm',1000*zsed(n,end,i,j));
   text(0.01,11,ts)
   ts = sprintf('AL thickness = %g mm',1000*active_layer(n,i,j));
   text(0.01,11.5,ts)

% Bene for Total mass depth vs. Fraction   
   subplot(244)
   cla
   zz2   = zeros(nbed,1);
   zz2I  = zeros(nbed,1);
   taucr = zeros(nbed,1);
   poro  = zeros(nbed,1);
   ff2=squeeze(frac(:,n,:,i,j))';
   zz2(:,1)=squeeze(bmz(n,:,i,j))';
   temp=0.;
   for m = nbed:-1:1 %1 : nbed
       taucr(nbed-m+1,:)=squeeze(bed_tau_crit(n,m,i,j));
       poro(nbed-m+1,:)=squeeze(bed_porosity(n,m,i,j));
       if m ~= 1, 
           val = zz2(m-1,1);
       else
           val = 0.;
       end
       temp1 = ff2(m,1);
       x(:,1) = [ 0;0;temp1;temp1 ];
       for ns = 2 : nsed
           temp2 = temp1+ff2(m,ns);
           x(:,ns) = [ temp1;temp1;temp2;temp2 ];
           temp1 = temp2;
       end
       y(:,1) = [ temp;zz2(nbed,1)-val;zz2(nbed,1)-val;temp ] ;
       for ns = 2 : nsed
           y(:,ns) = y(:,1) ;
       end
       temp = zz2(nbed,1)-val;
       zz2I(nbed-m+1,:)=temp; 
       pl1=patch(x(:,1),y(:,1),my_map(1,:));
       hold on
       for ns = 2 : nsed
           pl1=patch(x(:,ns),y(:,ns),my_map(ns,:));
       end
   end
   ax1 = gca;
   ylabel('Total mass depth (kg/m2)')
   set(ax1,'YAxisLocation','right')
   maxi=ceil(max(max(bmz(:,:,i,j))));
   xlabel('Fraction')
   set(gcf,'Colormap',my_map)
   prems=squeeze(bmz(1,end,i,j))';
   plot([0 1],[prems prems],'--k')
   axis([0 1 0 maxi])
   hold on
   
% Trick Matlab
   PosAx1 = get(ax1,'Position');
   minbp = min(min(min(min(bed_porosity))));
   maxbp = max(max(max(max(bed_porosity))));
   ax3=axes('Position',[PosAx1(1) PosAx1(2) PosAx1(3) 0.375],...
           'YColor',[1 1 1],...
           'XAxisLocation','top',...
           'Color','none',...
           'XColor',[0 1 1],...
           'YTickLabel','','YTick',zeros(1,0));
   axis([ minbp maxbp 0 maxi])
   text(maxbp + (maxbp-minbp)/10,maxi,...
       'Porosity','Color',[0 1 1])
   
   maxbt = max(max(max(max(bed_tau_crit))));
   ax2=axes('Position',PosAx1,...
           'XAxisLocation','top',...
           'Color','none',...
           'XColor',[0 1 0]);
   axis([0 maxbt 0 maxi])
   text(maxbt+maxbt/10.,maxi,...
       'Critical shear stress','Color',[0 1 0])   
    
% Plot the critical shear stress profile
   ax4=axes('Position',get(ax1,'Position'),...
           'XAxisLocation','top',...
           'Color','none',...
           'XColor',[0 1 0]);
   hl1=plot(taucr,zz2I,'Color',[0 1 0],'linewidth',2);
   axis([0 max(max(max(max(bed_tau_crit)))) 0 maxi])
   set(ax4,'visible','off')
   box on;
   
% Plot the porosity profile
   ax5=axes('Position',get(ax1,'Position'),...
           'XAxisLocation','top',...
           'Color','none');
   hl2=plot(poro,zz2I,'Color',[0 1 1],'linewidth',2);
   axis([min(min(min(min(bed_porosity)))) ...
       max(max(max(max(bed_porosity)))) 0 maxi])
   set(ax5,'visible','off')
   box on;   

   subplot(816)
   cla
%    plot(tdays,taumax(:,i,j),'-r','linewidth',2);
   plot(tdays,taub,'-r','linewidth',2);
   set(gca,'xlim',[tdays(1) tdays(end)],'XTickLabel',[]);
   ylabel('\tau_b (Pa)')
   hold on
%    hh=plot(tdays(n),taumax(n,i,j),'or');
   hh=plot(tdays(n),taub(n,1),'or');
   set(hh,'markerfacecolor',[1 0 0])
   set(gca,'xlim',[tdays(1) tdays(end)],'XTickLabel',[]);
   
   subplot(817)
   cla
   plot(tdays,1e6*squeeze(active_layer(:,i,j)),'-b','linewidth',2);
   hold on
   hh=plot(tdays(n),1e6*squeeze(active_layer(n,i,j)),'ob');
   set(hh,'markerfacecolor',[0 0 1])
   set(gca,'xlim',[tdays(1) tdays(end)],'XTickLabel',[]);
   ylabel('AL(\mum)')

   subplot(818)
   cla
   plot(tdays,squeeze(ero_flux(:,i,j)),'-g','linewidth',2);
   ylabel('Ero flux (kg/m^2)')
   hold on
   hh=plot(tdays(n),squeeze(ero_flux(:,i,j)),'og');
   set(hh,'markerfacecolor',[0 1 0])
   set(gca,'xlim',[tdays(1) tdays(end)],'XTickLabel',[]);

   set(gcf, 'Colormap', my_map)
   set(gcf, 'Renderer', 'Painters');
   set(gcf, 'color', 'white',...
   'Position',[1 70 1280 880]);
   drawnow
   shg
%    if(n==1), pause, end
   pause
   pfn = sprintf(' ./frames/f%05d.png',n);
   %anim_frameB('Mixed',n);
end