% layer_look_test5.m % % Look at bed time series from ROMS mudtoy % Plots sediment fraction, sediment mass, thickness, and critical shear % stress of each bed layer. Animates layers with timeseries of tau_b, % active layer thickness, suspended mass, and erosion flux. % csherwood@usgs.gov % 2007/09/05 % J.Paul Rinehimer - VIMS % 2007/06/06 % % Modified from layer_look.m % Benedicte Ferre - USGS % 01/16/07 %% % Set some parameters ncf = 'd:\bf\STRATAFORM\ROMS\bf_branch\coh_test\Data\Results\his_B6_bot.nc'; i = 2; j = 2; % Get data from netCDF files nco = netcdf(ncf) t = nco{'ocean_time'}(:); % dstart = nco{'dstart'}(:); nbeddim = nco('Nbed'); % number of bed sediment levels nzdim = nco('N'); % number of water-column levels nxdim = nco('xi_rho'); nydim = nco('eta_rho'); h = nco{'h'}(:); ero_flux = nco{'ero_flux'}(:); u = nco{'u'}(:,:,i,j); v = nco{'v'}(:,:,i,j); pcoh = nco{'cohesive_behavior'}(:); % Set parameters from netCDF files nt = length(t); nx = length(nxdim); ny = length(nydim); nz = length(nzdim); dt = nco{'dt'}(:); nbed = length(nbeddim); [nsed,nsand,nmud] = get_nseds(ncf); rhos = 2650 * ones([nsed 1]); % nsed=1; %% Get sediment data % Initialize sediment arrays 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); % Get cohesive fraction and mass from netcdf files if nmud~=0, for iSed = 1:nmud, iStr = sprintf('%02d',iSed); mud(iSed,:,:,:,:) = nco{['mud_' iStr]}(:); mudmass(iSed,:,:,:,:) = nco{['mudmass_' iStr]}(:); mudfrac(iSed,:,:,:,:) = nco{['mudfrac_' iStr]}(:); end end % Get noncohesive fraction and mass from netcdf if nsand~=0, for iSed = 1:nsand, iStr = sprintf('%02d',iSed); sand(iSed,:,:,:,:) = nco{['sand_' iStr]'}(:); sandmass(iSed,:,:,:,:) = nco{['sandmass_' iStr]}(:); sandfrac(iSed,:,:,:,:) = nco{['sandfrac_' iStr]}(:); end end %% % Get bed properties and bottom stress from netCDF files 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); taub =( abs(squeeze(bustr(:,i,j-1)))); % kludge alert!! % taumax = taub; % For no BBL model bustrcwmax = nco{'bustrcwmax'}(:,i,j); bvstrcwmax = nco{'bvstrcwmax'}(:,i,j); taumax = sqrt(bustrcwmax.^2+bvstrcwmax.^2); dep_net=nco{'dep_net'}(:); % Get vertical coordinate arrays %zr = zeros(nt,nz,ny,nx); %zw = zeros(nt,nz+1,ny,nx); zr = zeros(1,nz,ny,nx); zw = zeros(1,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 % This part takes way to long.... should modify get_z to return all time % values at once. % All the z_r, z_w are the same w.r.t time [z_r,z_w]=get_z(10,nco); zr(1,:,:,:)=z_r; zw(1,:,:,:)=z_w; % All dz are the same too, but we'll replicate to match other 4-D arrays, % in case grid is stretched later on. % dz(n,:,:,:)=diff(z_w); dz = repmat(diff(zw),[nt 1 1 1 ]); % ncclose all %% Calculate total suspended and bed mass % Combines cohesive/noncohesive to same data arrays % Initialize variables for combined arrays susmass = zeros(nsed,nt,ny,nx); % Depth-integrated suspended mass susmass_z = zeros(nsed,nt,nz,ny,nx); % Depth-variable suspended mass frac = zeros(nsed,nt,nbed,ny,nx); % Bed sediment fraction sedmass = zeros(nsed,nt,nbed,ny,nx); % Bed sediment mass sedmass_t = zeros(nsed,nt,ny,nx); % Depth-integrated bed sediment mass % Calculate totals in bed per class mudmass_t = sum(mudmass,3); % total mass of each mud class in bed sandmass_t = sum(sandmass,3); % total mass of each sand class in bed % Cohesive seds if nmud~=0 for iSed=1:nmud sdz = squeeze(mud(iSed,:,:,:,:)).*dz; susmass(iSed,:,:,:) = squeeze(sum(sdz,2)); susmass_z(iSed,:,:,:,:) = mud(iSed,:,:,:,:); sedmass(iSed,:,:,:,:) = mudmass(iSed,:,:,:,:); frac(iSed,:,:,:,:) = mudfrac(iSed,:,:,:,:); sedmass_t(iSed,:,:,:) = squeeze(mudmass_t(iSed,:,:,:,:)); end end % Noncohesives if nsand~=0, for iSand=1:nsand iSed = iSand+nmud; sdz = squeeze(sand(iSand,:,:,:,:)).*dz; susmass(iSed,:,:,:) = squeeze(sum(sdz,2)); susmass_z(iSed,:,:,:,:) = sand(iSand,:,:,:,:); sedmass(iSed,:,:,:,:) = sandmass(iSand,:,:,:,:); frac(iSed,:,:,:,:) = sandfrac(iSand,:,:,:,:); sedmass_t(iSed,:,:,:) = squeeze(sandmass_t(iSand,:,:,:,:)); end end susmasstot = squeeze(sum(susmass,1)); % total mass in suspension susmasstot_z = squeeze(sum(susmass_z,1)); % total mass in suspension per level sedmasstot = squeeze(sum(sedmass_t,1)); % total mass in bed thicki = sum(bed_thickness(1,:,:,:),2); % Initial sediment thickness zsed = cumsum(bed_thickness,2); % depth in sediment zsedb = flipdim(cumsum(flipdim(bed_thickness,2),2),2); % height above sediment bottom zsed0 = repmat(thicki,[nt,nbed,1,1])-zsedb; % Depth referenced to initial sediment thickness bmz = cumsum(squeeze(sum(sedmass,1)),2); % mass depth in sediment % this DOES equal mudmass+sandmass mass_calc = zeros(nsed,nt,nbed,ny,nx); for iSed=1:nsed mass_calc(iSed,:,:,:,:) = bed_thickness.*(1-bed_porosity).*rhos(iSed).*squeeze(frac(iSed,:,:,:,:)); % mass_calc(iSed,:,:,:,:) = squeeze(bed_thickness).*(1-squeeze(bed_porosity)).*rhos(iSed).*squeeze(frac(iSed,:,:,:,:)); end mass_tot = sum(mass_calc,3); total_mass = sedmass_t+susmass; % total sediment mass (should be 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 %% Setup animated figure figure(2); clf scrsz = get(0,'ScreenSize'); n=1; tdays = t./(24*3600); nsones = ones([nsed,1]); my_map = [0 0 1; 1 0 0 1 1 0; 0 1 0]; % Bottom Stress plot axh(6) = subplot(815); plot(tdays,taumax,'-r','linewidth',2); ylimit6=get(gca,'ylim'); text(0.1,ylimit6(2)+ylimit6(2)*0.15,'\tau_b (Pa)','Fontsize',12) hold on ph.tau = plot(tdays(n),taumax(n),'or'); hold off set(ph.tau,'MarkerFaceColor','r', ... 'xdatasource','tdays(n)', ... 'ydatasource','taumax(n)'); % Suspended mass plot axh(5) = subplot(816); area(tdays,squeeze(susmass(:,:,i,j)')) ylimit5=get(gca,'ylim'); text(0.1,ylimit5(2)+ylimit5(2)*0.2,'Sus. M (kg/m^2)','Fontsize',12) hold on ph.sm = plot(nsones.*tdays(n),cumsum(susmass(:,n,i,j)),'ko'); hold off set(ph.sm,'MarkerFaceColor','k', ... 'xdatasource','nsones.*tdays(n)', ... 'ydatasource','cumsum(susmass(:,n,i,j))') if nsed==4, legend({'Mud 4\mum','Mud 30\mum','Sand 62.5\mum','Sand 100\mum'}, ... 'Position',[0.2893 0.9509 0.433 0.03666],... 'Orientation','horizontal'); elseif nsed==1, legend({'Mud 4\mum'}, 'Position',[0.2893 0.9509 0.433 0.03666],... 'Orientation','horizontal'); elseif nsand>0 legend({'Sand 62.5\mum','Sand 100\mum'}, 'Position',[0.2893 0.9509 0.433 0.03666],... 'Orientation','horizontal'); elseif nmud>0 legend({'Mud 4\mum','Mud 30\mum'}, 'Position',[0.2893 0.9509 0.433 0.03666],... 'Orientation','horizontal'); end legend('boxoff') % Active Layer plot axh(7) = subplot(817); plot(tdays,1e6*squeeze(active_layer(:,i,j)),'-b','linewidth',2); ylimit7=get(gca,'ylim'); text(0.1,ylimit7(2)+ylimit7(2)*0.15,'AL Thick (\mum)','Fontsize',12) hold on ph.al = plot(tdays(1),1e6*squeeze(active_layer(1,i,j)),'ob'); hold off set(ph.al,'markerfacecolor','b',... 'xdatasource','tdays(n)',... 'ydatasource','1e6*squeeze(active_layer(n,i,j))') % net deposition/erosion axh(8) = subplot(818); plot(tdays,squeeze(dep_net(:,i,j))*1000,'-g','linewidth',2); ylimit8=get(gca,'ylim'); text(0.1,ylimit8(2)+ylimit8(2)*0.25,'Net deposition (kg/m^2)','Fontsize',12) hold on plot([tdays(1) tdays(end)],[0 0],'k') ph.ef = plot(tdays(1),squeeze(dep_net(1,i,j))*1000,'og'); hold off set(ph.ef,'markerfacecolor','g', ... 'xdatasource','tdays(n)', ... 'ydatasource','squeeze(dep_net(n,i,j))') xlabel('Time (days)','Fontsize',12) % Common axis properties set(axh(5:8),'xlim',[tdays(1) tdays(end)]); set(axh(5:7),'xticklabel',[]); % Figure Properties fprops.Colormap = my_map; fprops.Renderer = 'Painters'; fprops.color = 'white'; fprops.Position = [50 50 1000 700]; set(gcf,fprops) %% Animate the plot if(1) mkdir('frames') for n=1:10:nt % for n=300 % Layer/Fraction plot axh(1) = subplot(251); bfh = barh( squeeze(frac(:,n,:,i,j))' ,1,'stacked'); set(gca,'xlim',[0 1.15]) ylabel('Layer Number') xlabel('Fraction') % Layer/Mass plot axh(2) = subplot(252); bmh =barh( squeeze(sedmass(:,n,:,i,j))',1,'stacked'); set(gca,'xlim',[0 0.1]) xlabel('Mass (kg / m^2 )') ts = sprintf('Total = %5.2g kg/m^2',bmz(n,end,i,j)); text(.005,nbed+4,ts) % Layer/Thickness plot axh(3) = subplot(253); bth = barh( squeeze(1e3*bed_thickness(n,:,i,j))',1,'Facecolor',[0.5 0.5 0.5]); set(gca,'xlim',[0 0.2]) xlabel('Thickness (mm)') ts = sprintf('Total = %5.2g mm',1e3*zsed(n,end,i,j)); text(0.01,nbed+2.5,ts) ts = sprintf('AL thick = %5.2g mm',1e3*active_layer(n,i,j)); text(0.01,nbed+4,ts) % Critical Shear stress plot axh(4) = subplot(254); plot(bed_tau_crit(n,:,i,j),1e3.*zsed0(n,:,i,j),'d-r') tclims = [0 2]; ylims = [-1 2]; axis([tclims ylims]) hold on plot(tclims,[0 0],'k:') plot(taub(n)*[1 1],ylims,'k--') plot(tclims,1000 .* zsed0(n,1,i,j).*[1 1],'b--') xlabel('\tau_{crit}') text(-0.45,1.2,'Depth into Bed, mm','rotation',90) ts = sprintf('\\tau_b = %5.2g Pa',taub(n)); text(0.08,ylims(1)+0.15,ts) ts = sprintf('\\tau_c(1) = %5.2g Pa',bed_tau_crit(n,1,i,j)); text(0.08,ylims(1)+0.35,ts) ts = sprintf('\\tau_c(2) = %5.2g Pa',bed_tau_crit(n,2,i,j)); text(0.08,ylims(1)+0.55,ts) hold off % suspended sediment concentration axh(5) = subplot(255); for j=1:nsed, plot(squeeze(susmass_z(j,n,:,i,j)),h(i,j)+z_r(:,i,j),'LineWidth',2,'Color',my_map(j,:)) hold on end plot(squeeze(susmasstot_z(n,:,i,j)),h(i,j)+z_r(:,i,j),'LineWidth',2,'Color',[0 0 0]) hold off axis([0 1e-2 0 4]) xlabel('ssc kg/m^3') ylabel('Elevation, m') set(gca,'YAxisLocation','right') % Common axis properties set(axh(1:4),'ydir','reverse') set(axh(1:3),'ylim',[0 nbed+5]) refreshdata %set(ph.sm, 'xdata',tdays(n).*ones([nsed,1]),'ydata',cumsum(susmass(:,n,i,j))); %set(ph.tau,'xdata',tdays(n),'ydata',taumax(n)); %set(ph.al, 'xdata',tdays(n),'ydata',1e6*squeeze(active_layer(n,i,j))); if ~isempty(pcoh) if squeeze(pcoh(n,i,j))==1, col=[1 0 0]; % pure cohesive elseif squeeze(pcoh(n,i,j))==0, col=[0 0 1]; % pure non-cohesive else col=[0 1 1]; % mixed sediment end else col=[0 1 0]; end set(ph.ef, 'xdata',tdays(n),'ydata',squeeze(dep_net(n,i,j)),... 'MarkerFaceColor',col,'MarkerEdgeColor',col); drawnow shg %pause(.2) pfn = sprintf('f%05d.png',n) % anim_frameB('Mixed',n); eval(['print -dpng .\frames\',pfn]) %pause end end %% Plot mass conservation figure(3),clf smt = squeeze(susmasstot(:,i,j)); bmt = squeeze(sedmasstot(:,i,j)); tmt = smt+bmt; thrs = tdays*24; plot(thrs,smt,'r',thrs,bmt,'b',thrs,tmt,'--k') %,thrs([1 end]),ones(2,1).*tmt(1),'k:') xlabel('Hours') ylabel('Mass, kg m^{-2}') title('Mass Conservation?') legend('Suspended Mass','Bed Mass','Total Mass','Location', 'NorthWest') %% anim_make; %figure %pslice(repmat(tdays,1,20),squeeze(zw(:,:,i,j)),... % squeeze(susmasstot_z(:,:,i,j))*1000.) % ero_flux_GPu=ero_flux; % erosion_stress_GPu=erosion_stress; % maskCS_GPu=nco{'maskCS'}(:); % susmass_z_GPu=susmass_z; % susmasstot_z_GPu=susmasstot_z; % save results_mixed_GPu *_GPu