% MIXED_LOOK - Look at results of test with MIXED_TOY

% Benedicte Ferre - USGS
% 01/31/2007

% Parameter to change manually, defines the number of sediment classes
% Should be read from the data file further
nmud  = 0;
nsand = 1;
nsed=nmud+nsand;

%% Load what we are interested in from the netcdf file
 fnc='ocean_mixed_toy_his.nc';
%fnc = 'ocean_test_chan_his_ZF.nc';
nco = netcdf(fnc)
nxi = length(nco('xi_rho'));
neta = length(nco('eta_rho'));
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('s_rho');      % number of water-column levels (?)
nz =length(nzdim);
nxdim = nco('xi_rho');
nx = length(nxdim);
nydim = nco('eta_rho');
ny = length(nydim);
h=nco{'h'}(:);
x_rho=nco{'x_rho'}(:);
y_rho=nco{'y_rho'}(:);
rmask=nco{'mask_rho'}(:);
h(~rmask)=nan;
%jdmat=cf_time(fnc,'ubar');
jdmat = datenum(datestr(t));

diam=nco{'grain_diameter'}(:);
slope=nco{'dmix_slope'}(:);
maskCS=nco{'maskCS'}(:);
ero_flux=nco{'ero_flux'}(:);
dep_net=nco{'dep_net'}(:);

rhos =  2650.*ones(1,nsed);

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


%% calculate
    mud = zeros(nmud,nt,nz,neta,nxi);
    mudmass = zeros(nmud,nt,nbed,neta,nxi);
    mudfrac = zeros(nmud,nt,nbed,neta,nxi);
if nsand~=0, 
    sand = zeros(nsand,nt,nz,neta,nxi);
    sandmass = zeros(nsand,nt,nbed,neta,nxi);
    sandfrac = zeros(nsand,nt,nbed,neta,nxi);
end
if nmud==0 && nsand==0,
    disp('No sediment included')
    break
end

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
%%
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
susmasstot = squeeze(sum(susmass,1));    % total mass in suspension
susmasstot_z = squeeze(sum(susmass_z,1));    % total mass in suspension per level


%%
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'}(:);
bustr = nco{'bustrcwmax'}(:);
bvstr = nco{'bvstrcwmax'}(:);
bustr2=bustr.^2;
bvstr2=bvstr.^2;
tau=sqrt(bustr2+bvstr2);
tau(:,ny,nx)=0;
Hwave=nco{'Hwave'}(:);
mudmass_z = zeros(nsed,nt,1,ny,nx);
sandmass_z = zeros(nsed,nt,1,ny,nx);
sedmass_z = zeros(nsed,nt,1,ny,nx);

if nmud~=0,
    mudmass_z = sum(mudmass,3);    % total mass of each mud class in bed
    for n=1:nmud
        sedmass(n,:,:,:,:)=mudmass(n,:,:,:,:);
        sedmass_z(n,:,1,:,:)=mudmass_z(n,:,:,:,:);
        frac(n,:,:,:,:)=mudfrac(n,:,:,:,:);
    end
end
if nsand~=0,
    sandmass_z= sum(sandmass,3);   % total mass of each sand class in bed
    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,:,:,:,:));
% end

sedmass_ztot = squeeze(sum(sedmass_z,1));
% all = sedmass_ztot + susmasstot;         % this should not change if mass is conserved;


sscbot = squeeze(sum(mud(:,:,1,:,:)));

%% SSC transect at i=2
z=2;
temp=cumsum(m_lldist(x_rho(z,:),y_rho(z,:)));
for n=1:nx-1,
    range(n,:)=temp(nx-n);
end
TransSSC=zeros(nt+1,nx+1);
TransSSC(1,1:2)=0;
TransSSC(2:end,1)=jdmat;
TransSSC(1,3:end)=range';
% for n=1:length(t)
%     TransSSC(n+1,2:end)=squeeze(sscbot(n,z,:))';
% end
% plot(h(z,:),TransSSC(10,2:end))
% prop=get(gca,'YLim');
% axis([0 100 0 prop(2)])
% title('ssc at the bottom along the 3rd transect','FontSize',14);
% xlabel('Depth (m)','FontSize',12)
% ylabel('ssc (kg m^3)','FontSize',12)
% 
% 
% %% Shear stress transect at i=4
% figure
% TransTau=zeros(nt+1,nx+1);
% TransTau(1,1:2)=0;
% TransTau(2:end,1)=jdmat;
% TransTau(1,3:end)=range';
% for i=1:length(t)
%     TransTau(i+1,2:end)=tau(i,z,:);
% end
% plot(h(z,:),TransTau(10,2:end))
% title('Shear stress along the 3rd transect','FontSize',14);
% xlabel('Depth (m)','FontSize',12)
% ylabel('\tau (N/m^{2})','FontSize',12)
% prop=get(gca,'YLim');
% axis([0 100 0 prop(2)])

a=[0:size(x_rho)-1]';
a=repmat(a,1,length(x_rho));
Y=[1,1;1,1];

figure
%% Plot the 1st transect
for nt= length(t)
    set(gcf,'color','white');
    subplot(7,10,1:40)
    pslice(repmat(x_rho(1,:),20,1),squeeze(z_r(:,1,:)),...
        squeeze(susmasstot_z(nt,:,z,:)),[0. 60.])
    hold on
	ylabel(gca, {'Depth (m)'})
    title(['Sediment transport: ' datestr(jdmat(nt)) ' UTC']);
    axis(gca,[0 max(max(x_rho))-50 -max(max(h)) 0]);
    box(gca,'on');
    hold(gca,'all');
    colorbar1 = colorbar('peer', gca,...
        'Position',[0.851 0.516 0.048 0.433],...
        'Box','on',...
        'XLim',[-0.5 1.5]);
	set(gca,...
      'Position',[0.13 0.519 0.68 0.433],...
      'TickDir','out');
    
    subplot(7,10,41:50)
    dep=squeeze(dep_net(nt,:,:))*1000.;
%     res=squeeze(ero_flux(nt,:,:))*1000.;
    maxi=max(max(max(squeeze(dep_net(:,:,:)))))*1000.;
    mini=min(min(min(squeeze(dep_net(:,:,:)))))*1000.;
    pslice(x_rho,a,dep,[mini maxi])
    axis(gca,[0 max(max(x_rho))-50 2.5 3]);
    ylabel(gca,{'net deposition (mm)'});
    box(gca,'on');
    hold(gca,'all');
    colorbar3 = colorbar('peer', gca,...
        'Position',[0.851 0.359 0.048 0.07],...
        'Box','on',...
        'XLim',[-0.5 1.5]);
    set(gca,...
        'Position',[0.13 0.364 0.68 0.07],...
        'YTickLabel',{' '});


    subplot(7,10,51:59)
    mask=squeeze(maskCS(nt,:,:));
    psliceB(x_rho,a,mask,[0 1])
    axis(gca,[0 max(max(x_rho))-50 2.5 3]);
    ylabel(gca,{'Ero. method'});
    box(gca,'on');
    hold(gca,'all');
    set(gca,'Position',[0.13 0.233 0.68 0.07],...
        'YTickLabel',{' '});
    subplot(7,10,60)
    area(Y)
    set(gca, 'Position',[0.851 0.233 0.048 0.07],...
        'XTickLabel',{' '},...
        'YTickLabel',{' '});
    text(2.135,1.483,'NCS')
    text(2.173,0.417,'CS')    
    
    
    subplot(7,10,61:70)
    d=squeeze(diam(nt,:,:)*1000);
    pslice(x_rho,a,d,[0.03 0.17])
    axis(gca,[0 max(max(x_rho))-50 2.5 3]);
    ylabel(gca,{'diameter (mm)'});
    box(gca,'on');
    hold(gca,'all');
    colorbar3 = colorbar('peer', gca,...
        'Position',[0.851 0.11 0.048 0.07],...
        'Box','on',...
        'XLim',[-0.5 1.5]);
    set(gca,...
        'Position',[0.13 0.11 0.68 0.07],...
        'YTickLabel',{' '});
    
    anim_frameB('MixedProfil',nt);
end
ncclose all

% figure
% %% Plot the last cell
% for nt=1:length(t)
%     pslice(x_rho,y_rho,squeeze(diam(nt,:,:))*1000.,[0.04 0.12])
%     pause
% end