diff --git a/Edgeerosion.m b/Edgeerosion.m new file mode 100644 index 0000000..a73b672 --- /dev/null +++ b/Edgeerosion.m @@ -0,0 +1,68 @@ +function [deltaY2,Pedge,Y2OX,EdgeERY2]=EdgeerosionXXX(P,z,aw,maxedgeheight,fox,dt,dx,MASK,A,Y2OX); + + + +%MASK 1 are the mudflat +%MASK 2 is where the are no waves +%fox=0; + +%the non-periodic cell are walls +MASK(1,:)=0; +MASK(end,:)=0; + +%you get the wave power from the 4 sourronding cells +Pedge=[P(:,1)*0 P(:,1:end-1)]+[P(1,:)*0; P(1:end-1,:)]+[P(:,2:end) P(:,end)*0]+[P(2:end,:); P(end,:)*0]; +Pedge(MASK==1)=0;%the wave power in the mudflat becomes zero! +%find the marsh cells with some wave power around it +edg=find(A==1 & MASK==0 & Pedge>0); +%rng('shuffle'); +r=rand(length(edg),1);% you only need rng at the beginnign of the loop +a=find(r0 %yes, there are adjacent cells + dz=([z(edger(i))-z(p(a))]); + dz=mean(max(dz,0));%dz=min(max(dz,0)); %dz=max(max(dz,0)); + dz=min(dz,maxedgeheight); + %dz=2;%maxedgeheight; + + %this is how much the bed is lowered, includes everything!!! + deltaY2(edger(i))=deltaY2(edger(i))+dz; + + %This goes into resuspension. %This is what is conserved!!! + EdgeERY2(edger(i))=EdgeERY2(edger(i))+dz.*(1-fox); + + %Keep track of how much you oxydeized!!!! + Y2OX=Y2OX+dz.*fox; + + + end + +end + + + \ No newline at end of file diff --git a/REFERENCEMARSHr07.mat b/REFERENCEMARSHr07.mat new file mode 100644 index 0000000..79a5f9b Binary files /dev/null and b/REFERENCEMARSHr07.mat differ diff --git a/SeaWaves.m b/SeaWaves.m new file mode 100644 index 0000000..34e7cc2 --- /dev/null +++ b/SeaWaves.m @@ -0,0 +1,135 @@ +function [Um,TP,HS,F,kwave,PW]=SeaWaves(h,angle,hwSea_lim,range,wind,MASK,ndir,dx); + +%angle=0; +Um=0*h;TP=0*h;HS=0*h; + +extrafetch=5000;%[m} +Lbasin=1000/dx; +Fetchlim=max(50,dx*2);%dx*2;%600;%dx*2*10; +dlo=hwSea_lim; %minimum water depth to calculate wave. below this you don't calculate it + + + + +extra=1;%if (angle<90 | angle>270);extra=1;else;extra=1;end + + + + +%The standard way +if extra==0 +F=calculatefetch(MASK,ndir,dx,angle); +end +% +% %For Georgia +%extrafetch=0;%[m} +%F=calculatefetchWITHEXTRASonlyOCEAN(MASK,ndir,dx,angle,extrafetch); +% F1=calculatefetchWITHEXTRASonlyOCEAN(MASK,ndir,dx,angle,extrafetch); +% F2=calculatefetchWITHEXTRASonlyOCEAN(MASK,ndir,dx,mod(angle+90,360),extrafetch); +% F3=calculatefetchWITHEXTRASonlyOCEAN(MASK,ndir,dx,mod(angle+180,360),extrafetch); +% F4=calculatefetchWITHEXTRASonlyOCEAN(MASK,ndir,dx,mod(angle+270,360),extrafetch); +% F=(F1+F2+F3+F4)/4; + +%For the idealize basin +%extrafetch=10000;%[m} +if extra==1; +F=calculatefetchWITHEXTRAS(MASK,ndir,dx,angle,extrafetch,Lbasin,h-0.1-range/4); +end +% F1=calculatefetchWITHEXTRAS(MASK,ndir,dx,angle,extrafetch); +% F2=calculatefetchWITHEXTRAS(MASK,ndir,dx,mod(angle+90,360),extrafetch); +% F3=calculatefetchWITHEXTRAS(MASK,ndir,dx,mod(angle+180,360),extrafetch); +% F4=calculatefetchWITHEXTRAS(MASK,ndir,dx,mod(angle+270,360),extrafetch); +% F=(F1+F2+F3+F4)/4; + +Fo=F; + +% %For all the modified ways. Creates a buffer on the side +% %boundaries. Just used as a mask, the actual value is not +% %importnat, just need to be larger than fetchlim. +% %%%%%%%%%%%%%%%%%%%%%$%$#$&^$#^$&#^$#^$#&^$$*&^%&*%$*^%$%&*$*&%%$&*%$&%*$%&* +% %%%%%%%%%%%%%%%%%%%%%$%$#$&^$#^$&#^$#^$#&^$$*&^%&*%$*^%$%&*$*&%%$&*%$&%*$%&* +% %%%%%%%%%%%%%%%%%%%%%$%$#$&^$#^$&#^$#^$#&^$$*&^%&*%$*^%$%&*$*&%%$&*%$&%*$%&* +% [N,M]=size(h); +% Fo(2+floor(N*0.5):end-1,1:20)=9999; +% Fo(2+floor(N*0.5):end-1,end-20:end)=9999; +% %%%%%%%%%%%%%%%%%%% +% %%%%%%%%%%%%%%%%%%%%%$%$#$&^$#^$&#^$#^$#&^$$*&^%&*%$*^%$%&*$*&%%$&*%$&%*$%&* + + +F(Fo<=Fetchlim)=0; + +%usa questo per isolared la mudflat +%%%%%%%%%%%%%%%%%%%%%%%%%%% +if extra==1; +MASK(end-Lbasin:end,:)=1; +F(end-Lbasin:end,:)=extrafetch; +Fo(end-Lbasin:end,:)=extrafetch; +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%diffuse the fetch field +alphadiffusefetch=0.1; %messo 10 for the VCR wave validation 10;%0;%%%QUESTO ERA 1 FINO AD APRILE 23 2018!!!!! +F=diffusefetch(MASK,F,alphadiffusefetch,dx); +F(Fo<=Fetchlim | MASK==0)=0; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% figure +% imagesc(F) +% colormap('jet') +% caxis([0 max(F(:))]) +% pause +% + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +a=find(Fo>Fetchlim & h>dlo & F>0 & MASK==1);%h>dlo & %a=find(Fo>dx*2);%h>dlo & %a=find(h>dlo); +D=h(a);Ff=F(a); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +% %%%TRCUCCOZO TO AVOID depths too small +% hbedsheatresslim=0.5; +% h(h=tan(10/180*pi); +% value=value.*facNL; + + +%The standrd +%DO NOT CREEP AT THE POND EDGE +value(Spond(p(a))==1 & Spond(q(a))==0) = 0; +value(Spond(p(a))==0 & Spond(q(a))==1) = 0; + +% %trucco per evitare the i canali si chiudana +% %DO NOT CREEP AT THE POND EDGE +% %"it is a pond but it would be a channel, i.e., it is not a deep channel" %cit. Giulio M +% value((Spond(p(a))==1 & VEG(p(a))==1) & Spond(q(a))==0) = 0; +% value(Spond(p(a))==0 & (Spond(q(a))==1 & VEG(q(a))==1)) = 0; + + +S(p(a))=S(p(a))+value; %exit from that cell +i=[i;G(q(a))]; j=[j;G(p(a))]; s=[s;-value]; %gain from the neigborh cell +end + +%summary of the material that exits the cell +i=[i;G(p)]; j=[j;G(p)]; s=[s;1+S(p)]; +ds2 = sparse(i,j,s);%solve the matrix inversion +P=ds2\rhs;z(G>0)=full(P(G(G>0))); + + diff --git a/calculatefetch.m b/calculatefetch.m new file mode 100644 index 0000000..3235788 --- /dev/null +++ b/calculatefetch.m @@ -0,0 +1,78 @@ +function F=calculatefetch(A,ndir,dx,angle) +angle=angle-180; +[N,M]=size(A); +F=zeros(N,M); +A(isnan(A))=0; %any NaN is a boundary, that is, a 0 + +%Every boundary is a wall, need to do this so the fetch does not warp around! +A(1,1:end)=0;A(end,1:end)=0;A(1:end,1)=0;A(1:end,end)=0; + +di=1+mod(floor(angle/360*ndir+0.5),ndir); +alfa=(di-1)/ndir*2*pi; +m=max(abs(cos(alfa)),abs(sin(alfa))); +if (di<=(ndir*1/8) | (di>ndir*3/8 & di<=ndir*5/8) | di>ndir*7/8) + IND=[1+mod(round([1:N]'*cos(alfa)/m)*ones(1,M),N)]+[mod(ones(N,1)*[1:M]+round([1:N]'*sin(alfa)/m)*ones(1,M)-1,M)]*N; +else + IND=([1+mod(ones(M,1)*[1:N]+round([1:M]'*cos(alfa)/m)*ones(1,N)-1,N)]+[mod(+round([1:M]'*sin(alfa)/m)*ones(1,N),M)]*N); +end + +F(IND)=cumsumreset(A(IND))/m*dx; + +%at the boundary, impose the fetch of the nearby cell +F(1,1:end)=F(2,1:end);F(end,1:end)=F(end-1,1:end);F(1:end,1)=F(1:end,2);F(1:end,end)=F(1:end,end-1); + +% F1=[F(:,1) F(:,1:end-1)]; +% F2=[F(1,:); F(1:end-1,:)]; +% F3=[F(:,2:end) F(:,end)]; +% F4=[F(2:end,:); F(end,:)]; +% F=(F+F1+F2+F3+F4)/5; + + +%NM=max(M,N); +% P1=zeros(NM,ndir);P2=zeros(NM,ndir); +% for i=1:ndir; +% dir=(i-1)/ndir*2*pi; +% m=max(abs(cos(dir)),abs(sin(dir))); +% mm(i)=m; +% P1(:,i)=-round([1:NM]*cos(dir)/m); +% P2(:,i)=round([1:NM]*sin(dir)/m); +% end +% + +% di=1+mod(floor(angle/360*ndir+0.5),ndir); +% alfa=(di-1)/ndir*2*pi; +% m=max(abs(cos(alfa)),abs(sin(alfa))); +% if (di<=(ndir*1/8) | (di>ndir*3/8 & di<=ndir*5/8) | di>ndir*7/8) +% IND=[1+mod(-P1(1:N,di)*ones(1,M),N)]+[mod(ones(N,1)*[1:M]-1-P2(1:N,di)*ones(1,M),M)]*N; +% else +% IND=([1+mod(ones(M,1)*[1:N]-1-P1(1:M,di)*ones(1,N),N)]+[mod(-P2(1:M,di)*ones(1,N),M)]*N); +% end +% +% F(IND)=cumsumreset(A(IND),1)/mm(di)*dx; + + + + + +% di=1+mod(floor(angle/360*ndir+0.5),ndir); +% alfa=(i-1)/ndir*2*pi; +% m=max(abs(cos(alfa)),abs(sin(alfa))); +% +% if (di<=(ndir*1/8) | (di>ndir*3/8 & di<=ndir*5/8) | di>ndir*7/8) +% IND=[1+mod(round([1:N]'*cos(alfa)/m)*ones(1,M),N)]+[mod(ones(N,1)*[1:M]-1-round([1:N]'*sin(alfa)/m)*ones(1,M),M)]*N; +% else +% IND=([1+mod(ones(M,1)*[1:N]-1+round([1:M]'*cos(alfa)/m)*ones(1,N),N)]+[mod(-round([1:M]'*sin(alfa)/m)*ones(1,N),M)]*N); +% end +% +% F(IND)=cumsumreset(A(IND),1)/mm(di)*dx; + + +% P1=zeros(NM,ndir);P2=zeros(NM,ndir); +% for i=1:ndir; +% dir=(i-1)/ndir*2*pi; +% m=max(abs(cos(dir)),abs(sin(dir))); +% mm(i)=m; +% P1(:,i)=-round([1:NM]*cos(dir)/m); +% P2(:,i)=round([1:NM]*sin(dir)/m); +% end + diff --git a/calculatefetchWITHEXTRAS.m b/calculatefetchWITHEXTRAS.m new file mode 100644 index 0000000..6576e16 --- /dev/null +++ b/calculatefetchWITHEXTRAS.m @@ -0,0 +1,167 @@ +function F=calculatefetchWITHEXTRAS(A,ndir,dx,angle,extrafetch,Lbasin,h) +%angle=315-180-0. +angleo=angle; +%angle=mod(angle-180,360); +angle=angle-180; +[N,M]=size(A); +F=zeros(N,M); +A(isnan(A))=0; %any NaN is a boundary, that is, a 0c + +%Every boundary is a wall, need to do this so the fetch does not warp around! +A(1,1:end)=0;A(end,1:end)=0;A(1:end,1)=0;A(1:end,end)=0; + + +%di=1+mod(floor(angle/360*ndir+0.5),ndir); +di=1+mod(floor(angle/360*ndir),ndir);%modified Dec 2019 +alfa=(di-1)/ndir*2*pi; +m=max(abs(cos(alfa)),abs(sin(alfa))); +if (di<=(ndir*1/8) | (di>ndir*3/8 & di<=ndir*5/8) | di>ndir*7/8) + IND=[1+mod(round([1:N]'*cos(alfa)/m)*ones(1,M),N)]+[mod(ones(N,1)*[1:M]+round([1:N]'*sin(alfa)/m)*ones(1,M)-1,M)]*N; +else + IND=([1+mod(ones(M,1)*[1:N]+round([1:M]'*cos(alfa)/m)*ones(1,N)-1,N)]+[mod(+round([1:M]'*sin(alfa)/m)*ones(1,N),M)]*N); +end + +%F(IND)=cumsumreset(A(IND),1)/m*dx; + +%error +if angleo<0 | angleo>360;A='error' ;end + +% %if (angleo<44 | angleo>316) +% if (angleo<134 | angleo>226) +% F(IND)=cumsumresetEXTRA(A(IND),extrafetch/dx)/m*dx; +% else +% F(IND)=cumsumreset(A(IND))/m*dx; +% end + +padding=Lbasin*2;%must be larger than fetchlim!!! +if (angleo<45 | angleo>315) +F(IND)=cumsumresetEXTRA(A(IND),extrafetch/dx)/m*dx; + if angleo<44; + %ll=length(F(2:end-1-floor(N*0.5),end)); + %F(2+floor(N*0.5):end-1,end-padding:end)=2*extrafetch*[1:ll]'/ll*ones(padding+1,1)'; + %F(end,end-padding:end)=extrafetch; + a=find(h(:,end)<0);if a>0;Lside=a(end)-1;else;Lside=1;end + F(Lside:end,end-100:end)=max(extrafetch,F(Lside:end,end-100:end)); %TOGLI PER EVITRARE IL FETCH ALTO NELLE MUDFLAT SEAWARD + %F(Lside:end,end-10:end)=max(extrafetch,F(Lside:end,end-10:end)); + %F(:,end)=extrafetch; + else + %ll=length(F(2:end-1-floor(N*0.5),1)); + %F(2+floor(N*0.5):end-1,1:1+padding)=2*extrafetch*[1:ll]'/ll*ones(padding+1,1)'; + %F(1,end-padding:end)=extrafetch; + %F(1,:)=extrafetch; + a=find(h(:,1)<0);if a>0;Lside=a(end)-1;else;Lside=1;end + F(Lside:end,1:101)=max(extrafetch,F(Lside:end,1:101)); %TOGLI PER EVITRARE IL FETCH ALTO NELLE MUDFLAT SEAWARD + %F(Lside:end,1:11)=max(extrafetch,F(Lside:end,1:11)); + end + + +elseif (angleo>=45 & angleo<90)%134) +a=find(h(:,end)<0);if a>0;Lside=N-a(end)-1;else;Lside=N-1;end +F(IND)=cumsumresetEXTRAlateral1(A(IND),extrafetch/dx,Lside)/m*dx; +%(IND)=cumsumreset(A(IND))/m*dx; +%a=find(MASK( + F(end-Lside:end,:)=extrafetch; %offshore boudnary + %ll=length(F(2:end-1-N/2,end)); + %F(2+N/2:end-1,end-padding:end)=extrafetch*[1:ll]'/ll*ones(padding+1,1)'; + +elseif (angleo>275 & angleo<=315) +a=find(h(:,1)<0);if a>0;Lside=N-a(end)-1;else;Lside=N-1;end +F(IND)=cumsumresetEXTRAlateral1(A(IND),extrafetch/dx,Lside)/m*dx; +%F(IND)=cumsumreset(A(IND))/m*dx; + F(end-Lside:end,:)=extrafetch; %offshore boudnary + %ll=length(F(2:end-1-N/2,1)); + %F(2+N/2:end-1,1:1+padding)=extrafetch*[1:ll]'/ll*ones(padding+1,1)'; + + + +else +F(IND)=cumsumreset(A(IND))/m*dx; +end + +% padding=5; +% if (angleo<44 | angleo>316) +% F(IND)=cumsumresetEXTRA(A(IND),extrafetch/dx)/m*dx; +% if angleo<44; +% ll=length(F(2:end-1,end)); +% F(2:end-1,end-padding:end)=extrafetch*[1:ll]'/ll*ones(padding+1,1)'; +% else +% ll=length(F(2:end-1,1)); +% F(2:end-1,1:1+padding)=extrafetch*[1:ll]'/ll*ones(padding+1,1)'; +% end +% elseif (angleo>44 & angleo<90)%134) +% F(IND)=cumsumresetEXTRAlateral1(A(IND),extrafetch/dx)/m*dx; +% F(end-padding:end,:)=2*extrafetch; +% elseif (angleo>270 & angleo<316) +% F(IND)=cumsumresetEXTRAlateral1(A(IND),extrafetch/dx)/m*dx; +% F(end-padding:end,:)=extrafetch; +% else +% F(IND)=cumsumreset(A(IND))/m*dx; +% end + + + + +% +% figure +% imagesc(F) +% pause +% + +%at the boundary, impose the fetch of the nearby cell +F(1,1:end)=F(2,1:end);F(end,1:end)=F(end-1,1:end);F(1:end,1)=F(1:end,2);F(1:end,end)=F(1:end,end-1); + +% F1=[F(:,1) F(:,1:end-1)]; +% F2=[F(1,:); F(1:end-1,:)]; +% F3=[F(:,2:end) F(:,end)]; +% F4=[F(2:end,:); F(end,:)]; +% F=(F+F1+F2+F3+F4)/5; + + +%NM=max(M,N); +% P1=zeros(NM,ndir);P2=zeros(NM,ndir); +% for i=1:ndir; +% dir=(i-1)/ndir*2*pi; +% m=max(abs(cos(dir)),abs(sin(dir))); +% mm(i)=m; +% P1(:,i)=-round([1:NM]*cos(dir)/m); +% P2(:,i)=round([1:NM]*sin(dir)/m); +% end +% + +% di=1+mod(floor(angle/360*ndir+0.5),ndir); +% alfa=(di-1)/ndir*2*pi; +% m=max(abs(cos(alfa)),abs(sin(alfa))); +% if (di<=(ndir*1/8) | (di>ndir*3/8 & di<=ndir*5/8) | di>ndir*7/8) +% IND=[1+mod(-P1(1:N,di)*ones(1,M),N)]+[mod(ones(N,1)*[1:M]-1-P2(1:N,di)*ones(1,M),M)]*N; +% else +% IND=([1+mod(ones(M,1)*[1:N]-1-P1(1:M,di)*ones(1,N),N)]+[mod(-P2(1:M,di)*ones(1,N),M)]*N); +% end +% +% F(IND)=cumsumreset(A(IND),1)/mm(di)*dx; + + + + + +% di=1+mod(floor(angle/360*ndir+0.5),ndir); +% alfa=(i-1)/ndir*2*pi; +% m=max(abs(cos(alfa)),abs(sin(alfa))); +% +% if (di<=(ndir*1/8) | (di>ndir*3/8 & di<=ndir*5/8) | di>ndir*7/8) +% IND=[1+mod(round([1:N]'*cos(alfa)/m)*ones(1,M),N)]+[mod(ones(N,1)*[1:M]-1-round([1:N]'*sin(alfa)/m)*ones(1,M),M)]*N; +% else +% IND=([1+mod(ones(M,1)*[1:N]-1+round([1:M]'*cos(alfa)/m)*ones(1,N),N)]+[mod(-round([1:M]'*sin(alfa)/m)*ones(1,N),M)]*N); +% end +% +% F(IND)=cumsumreset(A(IND),1)/mm(di)*dx; + + +% P1=zeros(NM,ndir);P2=zeros(NM,ndir); +% for i=1:ndir; +% dir=(i-1)/ndir*2*pi; +% m=max(abs(cos(dir)),abs(sin(dir))); +% mm(i)=m; +% P1(:,i)=-round([1:NM]*cos(dir)/m); +% P2(:,i)=round([1:NM]*sin(dir)/m); +% end + diff --git a/cumsumreset.m b/cumsumreset.m new file mode 100644 index 0000000..e29528a --- /dev/null +++ b/cumsumreset.m @@ -0,0 +1,5 @@ +function F=cumsumreset(A) + +a=find(A==0); +A(a) = 1-diff([0; a]); +F=cumsum(A,1); diff --git a/cumsumresetEXTRA.m b/cumsumresetEXTRA.m new file mode 100644 index 0000000..bf774d2 --- /dev/null +++ b/cumsumresetEXTRA.m @@ -0,0 +1,27 @@ +function F=cumsumreset(A,extrafetch) + +% a=find(A==0); +% A(a) = 1-diff([0; a]); +% F=cumsum(A,1); + +S=A; +S(1,2:end-1)=extrafetch; +S(S==0)=NaN; +S(S==1)=0; +G=cumsum(S,1); +G(isnan(G))=0; + + +a=find(A==0); +A(a) = 1-diff([0; a]); +F=cumsum(A,1); + +F=F+G; +% % S=S(2:end,2:end-1); +% % [x,y]=find(S==0); +% % +% % figure;plot(x,y,'.') +% % figure;plot(x) +% figure;imagesc(S); +% figure;imagesc(G); +% pause \ No newline at end of file diff --git a/cumsumresetEXTRAlateral1.m b/cumsumresetEXTRAlateral1.m new file mode 100644 index 0000000..e4570d5 --- /dev/null +++ b/cumsumresetEXTRAlateral1.m @@ -0,0 +1,33 @@ +function F=cumsumreset(A,extrafetch,Lbasin) + +% a=find(A==0); +% A(a) = 1-diff([0; a]); +% F=cumsum(A,1); + +S=A; + +%ll=length(S(1,2:end-1)); +%S(1,2:end-1)=extrafetch*[1:ll]/ll; + +S(1,end-Lbasin:end-1)=extrafetch; +%S(1,:)=extrafetch; + +S(S==0)=NaN; +S(S==1)=0; +G=cumsum(S,1); +G(isnan(G))=0; + + +a=find(A==0); +A(a) = 1-diff([0; a]); +F=cumsum(A,1); + +F=F+G; +% % S=S(2:end,2:end-1); +% % [x,y]=find(S==0); +% % +% % figure;plot(x,y,'.') +% % figure;plot(x) +% figure;imagesc(S); +% figure;imagesc(G); +% pause \ No newline at end of file diff --git a/diffuseedgesediments.m b/diffuseedgesediments.m new file mode 100644 index 0000000..5960307 --- /dev/null +++ b/diffuseedgesediments.m @@ -0,0 +1,30 @@ +function F=diffusefetch(A,F,alpha,dx); + + +D=A.*(alpha*3600*24)/(dx^2); + +G=0*F; +p=find(A==1);%exclude the NOLAND CELLS (A==0) +NN=length(p);G(p)=[1:NN];rhs=F(p);[m,n]=size(G);i=[];j=[];s=[]; + +S=0*G; +[row col]=ind2sub(size(A),p); +for k = [m -1 1 -m] +%avoid to the the cells out of the domain (risk to make it periodic...) +if k==m;aa=find(col+1<=n);end;if k==-m;aa=find(col-1>0);end;if k==-1;aa=find(row-1>0);end;if k==1;aa=find(row+1<=m);end; + +q=p+k;%the translated cell +a=aa(A(q(aa))==1 );%only inclued the cells in whcih you can creep to + +value=min(D(p(a)),D(q(a)));%value=(D(p(a))+D(q(a)))/2.*facNL; + +S(p(a))=S(p(a))+value; %exit from that cell +i=[i;G(q(a))]; j=[j;G(p(a))]; s=[s;-value]; %gain from the neigborh cell +end + +%summary of the material that exits the cell +i=[i;G(p)]; j=[j;G(p)]; s=[s;1+S(p)]; +ds2 = sparse(i,j,s);%solve the matrix inversion +P=ds2\rhs;F(G>0)=full(P(G(G>0))); + + diff --git a/diffusefetch.m b/diffusefetch.m new file mode 100644 index 0000000..3d44509 --- /dev/null +++ b/diffusefetch.m @@ -0,0 +1,31 @@ +function F=diffusefetch(A,F,alpha,dx); + + +D=A*(alpha*3600*24)/(dx^2); + + +G=0*F; +p=find(A==1);%exclude the NOLAND CELLS (A==0) +NN=length(p);G(p)=[1:NN];rhs=F(p);[m,n]=size(G);i=[];j=[];s=[]; + +S=0*G; +[row col]=ind2sub(size(A),p); +for k = [m -1 1 -m] +%avoid to the the cells out of the domain (risk to make it periodic...) +if k==m;aa=find(col+1<=n);end;if k==-m;aa=find(col-1>0);end;if k==-1;aa=find(row-1>0);end;if k==1;aa=find(row+1<=m);end; + +q=p+k;%the translated cell +a=aa(A(q(aa))==1 );%only inclued the cells in whcih you can creep to + +value=min(D(p(a)),D(q(a)));%value=(D(p(a))+D(q(a)))/2.*facNL; + +S(p(a))=S(p(a))+value; %exit from that cell +i=[i;G(q(a))]; j=[j;G(p(a))]; s=[s;-value]; %gain from the neigborh cell +end + +%summary of the material that exits the cell +i=[i;G(p)]; j=[j;G(p)]; s=[s;1+S(p)]; +ds2 = sparse(i,j,s);%solve the matrix inversion +P=ds2\rhs;F(G>0)=full(P(G(G>0))); + + diff --git a/drainpond.m b/drainpond.m new file mode 100644 index 0000000..2d93fd8 --- /dev/null +++ b/drainpond.m @@ -0,0 +1,38 @@ +function [S,AC]=drainpond(z,A,S,N,M,dx,dt,zntw,distdr); + +AC=A*0;%the channel netwrk. 0 if not a channel netwrk 1 is a channel network + +%need to add to the driange systems those channel that are created in a transitory way%%%%% +CC = bwconncomp((z>-zntw),8); %find everything that migth be a channel: deep and connected to the main channel network +for i=1:CC.NumObjects + ind=CC.PixelIdxList{i}; + a=find(A(ind)==2);%if that blob includes a cell connect to a channel (A==2) + if a>0;AC(ind)=1;end %all the elements become a channel +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%G=AC; +%%%%%%%%%MAKE IT FATTER BY A CERTAIN LAG +fat=floor(distdr/dx); %how much thick to make the rim +for lag=1:fat + p = find(AC==1);[row col]=ind2sub(size(A),p); + for k = [N -1 1 -N] + %avoid to the the cells out of the domain (risk to make it periodic...) + if k==N;a=find(col+1<=M);end; + if k==-N;a=find(col-1>0);end; + if k==-1;a=find(row-1>0);end; + if k==1;a=find(row+1<=N);end; + q=p+k;%the translated cell + AC(q(a))=1; %make the translated cell a channel network + end +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +CC = bwconncomp((S==1),8); %find the isolated ponds +for i=1:CC.NumObjects + ind=CC.PixelIdxList{i}; + a=find(AC(ind)==1);%if that pond includes a cell within the channel network + if a>0;S(ind)=0;end %all the elements in the pond become drained (S==0) +end + diff --git a/excludeboundarycell.m b/excludeboundarycell.m old mode 100755 new mode 100644 diff --git a/findisolatedponds.m b/findisolatedponds.m new file mode 100644 index 0000000..62be5a6 --- /dev/null +++ b/findisolatedponds.m @@ -0,0 +1,137 @@ +function [S,AC,DIF]=findisolatedponds(Z,A,N,M,dx,Zntw,Zsill,distdr,minponddepth); + +% %Find everything that is isolated at least by the ponddepth threhsold +% Zold=Z; +% +% %This is to avoid pond driange on the sides! +% ZB=A*0;ZB(:,1)=1;ZB(:,end)=1;ZB(1,:)=1;ZB(end,:)=1; +% Z(ZB==1 & A==1)=10; +% Z(A==0)=10; +% +% Z=max(Zntw,Z);%below Znet is just a mudflat +% Z=min(Zsill,Z);%flooding sill height +% ZZ=Z;ZZ(isnan(ZZ))=0; +% %ZZ=imfill(ZZ,8); %THIS DOES ALL THE HEAVY LIFTING!!! ZZ=imfill(ZZ); %THIS DOES ALL THE HEAVY LIFTING!!! +% ZZ=imfill(ZZ,4); %THIS DOES ALL THE HEAVY LIFTING!!! ZZ=imfill(ZZ); %THIS DOES ALL THE HEAVY LIFTING!!! +% DIF = ZZ-Z;%the depth of pond filling +% S=DIF>minponddepth;%what constitutes a pond +% +% AC=(Zold0.001)=ZZ(DIF>0.001)-Zold(DIF>0.001); %tHIS IS THE ACTUAL water depth in the pond + S=DIF>minponddepth;%what constitutes a pond + + +AC=(Zoldminponddepth;%what constitutes a pond +% +% %AC=(-zold=minponddepth);%%IF<10^-8 +% + + + + + + + + + + + + + + + + + +% figure +% imagesc(S) +% pause +%AC=A*0;%the channel netwrk. 0 if not a channel netwrk 1 is a channel network +%Need to add to the driange systems those channel that are created as the model evolves%%%%% +% CC = bwconncomp((z>-zntw),8); %find everything that migth be a channel: deep and connected to the main channel network +% for i=1:CC.NumObjects +% ind=CC.PixelIdxList{i}; +% a=find(A(ind)==2);%if that blob includes a cell connect to a channel (A==2) +% if a>0;AC(ind)=1;end %all the elements become a channel +% end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%G=AC; +% %%%%%%%%%MAKE IT FATTER BY A CERTAIN LAG +% if distdr>0 +% fat=floor(distdr/dx); %how much thick to make the rim +% for lag=1:fat +% p = find(AC==1);[row col]=ind2sub(size(A),p); +% for k = [N -1 1 -N] +% %avoid to the the cells out of the domain (risk to make it periodic...) +% if k==N;a=find(col+1<=M);end; +% if k==-N;a=find(col-1>0);end; +% if k==-1;a=find(row-1>0);end; +% if k==1;a=find(row+1<=N);end; +% q=p+k;%the translated cell +% AC(q(a))=1; %make the translated cell a channel network +% end +% end +% end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +% CC = bwconncomp((S==1),8); %find the isolated ponds +% for i=1:CC.NumObjects +% ind=CC.PixelIdxList{i}; +% a=find(AC(ind)==1);%if that pond includes a cell within the channel network +% if a>0;S(ind)=0;end %all the elements in the pond become drained (S==0) +% end + diff --git a/getwaterdepth.m b/getwaterdepth.m old mode 100755 new mode 100644 index 42ef54b..a662071 --- a/getwaterdepth.m +++ b/getwaterdepth.m @@ -1,18 +1,19 @@ -function [h,ho,fTide,dtide,dHW,wl]=getwaterdepth(range,msl,z,kro); +function [h,ho,fTide,dtide,dsurge,dHW,wl,wlo]=getwaterdepth(range,msl,z,kro,hpRIV); DH=range; %total water displacement -dHW=max(0,z+msl+range/2);%water depth at MHW -dtide=max(0,z+msl+range/2);%water depth at MHW +dHW=max(0,-z+hpRIV+msl+range/2);%water depth at MHW +dtide=max(0,-z+hpRIV+msl+range/2);%water depth at MHW +dsurge=max(0,-z+hpRIV+msl);%water depth at MHW fTide=min(1,max(10^-3,dHW/DH));%hydroperiod -h=0.5*(max(0,dHW)+max(0,dHW-DH)); %for the flow -%ho=h; %the original water depth without the limiter -h(h(msl+range/2))=0; -wl=h-z; -% figure; -% imagesc(ho) -% pause \ No newline at end of file diff --git a/imposeboundarydepth.m b/imposeboundarydepth.m new file mode 100644 index 0000000..dfc99bf --- /dev/null +++ b/imposeboundarydepth.m @@ -0,0 +1,23 @@ +function d=imposeboundarydepth(A,d,opt,dimposed); + +if opt ==1 + +[m,n] = size(A); +%boundary cell impose depth +p = find(A==1 | A==21);[row col]=ind2sub(size(A),p); +for k = [m -1 1 -m] +if k==m;aa=find(col+1<=n);end;if k==-m;aa=find(col-1>0);end;if k==-1;aa=find(row-1>0);end;if k==1;aa=find(row+1<=m);end; +q=p+k;%the translated cell +a=aa(A(q(aa))==2);%only select the cells where you imposed co (A==2) +d(q(a))=d(p(a)); +end + + +elseif opt==2 +d=d;% +%or you can imposed the depth you want +%d(A==2 | A==22)=dimposed; + +end + + \ No newline at end of file diff --git a/initializegeometry.m b/initializegeometry.m new file mode 100644 index 0000000..dd9163c --- /dev/null +++ b/initializegeometry.m @@ -0,0 +1,43 @@ +function [N,M,dx,A,z,Active,x,y,msl]=initializegeometry(P); + +% %geometric parameters + +dx=5*2; +N=1000/2*1;%x +M=600/2*1;%y + +%%%%%%%%%%%cell types +A=ones(N,M); + +%sea b.c. +A(end,:)=2; + +%bathymetry +%initial profile. ELEVATION AT MSL +x=[0:N-1]*dx/1000;y=[0:M-1]*dx/1000; + + +slope=1/1*0.5/1000;z=[0:N-1]*slope*dx;z=z'*ones(1,M)-1+1;%sloping + +%channel in the middle +%z(end-100:end,M/2-5:M/2+5)=50; + + +%random +z=z+2*(rand(N,M)-0.5)*0.2;% *0.5 *0.2; + +z(end-1:end,:)=10; + +z=-z; + +z(A==0)=NaN; +%%%%%%%%%%%%%%%%%% + + + + +msl=0; +zbedo=-z-msl; +Active=zbedo0 +% z(q(er(aa(i))))=z(q(er(aa(i))))+dz; +% S(q(er(aa(i))))=1;%also update the pond ID +% pondlost=pondlost-dz; +% end +% end + +end + +deltaY=zoriginal-z; +deltaY2=deltaY; + + + + + +% function [A,d]=isolatedpondexpansion(A,d,m,n,dx,dt,MF); +% +% p = find(A==1);[row col]=ind2sub(size(A),p); +% +% for k = [m -1 1 -m] +% +% %avoid to the the cells out of the domain (risk to make it periodic...) +% if k==m;aa=find(col+1<=n);end; +% if k==-m;aa=find(col-1>0);end; +% if k==-1;aa=find(row-1>0);end; +% if k==1;aa=find(row+1<=m);end; +% +% q=p+k;%the translated cell +% er=aa(A(q(aa))==3); +% +% rng('shuffle');r=rand(length(er),1);% you only need rng at the beginnign of the loop +% a=find(r<0.1/365*dt*MF/dx); +% +% A(p(er(a)))=3; +% end \ No newline at end of file diff --git a/isolatedpondexpansion.m b/isolatedpondexpansion.m new file mode 100644 index 0000000..75611d4 --- /dev/null +++ b/isolatedpondexpansion.m @@ -0,0 +1,80 @@ +function [S,deltaY2,pondlost]=isolatedpondexpansion(z,S,A,N,M,dx,dt,zpondcr,maxdpond,aPEXP,pondlost); + +zoriginal=z; + +p=find(S==1);%find the existing ponds + +[row col]=ind2sub([N M],p); + +for k = [N -1 1 -N] + +%avoid to the the cells out of the domain (risk to make it periodic...) +if k==N;a=find(col+1<=M);end; +if k==-N;a=find(col-1>0);end; +if k==-1;a=find(row-1>0);end; +if k==1;a=find(row+1<=N);end; + +q=p+k;%the translated cell +er=a(A(q(a))==1 & S(q(a))==0); %find the neightr cell that is a non-pond cell + + +%rng('shuffle'); +r=rand(length(er),1);% you only need rng at the beginnign of the loop +aa=find(r0 + z(q(er(aa(i))))=z(q(er(aa(i))))-dz; + S(q(er(aa(i))))=1;%also update the pond ID + pondlost=pondlost+dz; + end +end + +end + +deltaY=zoriginal-z; +deltaY2=deltaY; + +% figure +% imagesc(deltaY2) +% pause + + + +% function [A,d]=isolatedpondexpansion(A,d,m,n,dx,dt,MF); +% +% p = find(A==1);[row col]=ind2sub(size(A),p); +% +% for k = [m -1 1 -m] +% +% %avoid to the the cells out of the domain (risk to make it periodic...) +% if k==m;aa=find(col+1<=n);end; +% if k==-m;aa=find(col-1>0);end; +% if k==-1;aa=find(row-1>0);end; +% if k==1;aa=find(row+1<=m);end; +% +% q=p+k;%the translated cell +% er=aa(A(q(aa))==3); +% +% rng('shuffle');r=rand(length(er),1);% you only need rng at the beginnign of the loop +% a=find(r<0.1/365*dt*MF/dx); +% +% A(p(er(a)))=3; +% end \ No newline at end of file diff --git a/mainevolutionstep.m b/mainevolutionstep.m old mode 100755 new mode 100644 index a63c6ea..19b569b Binary files a/mainevolutionstep.m and b/mainevolutionstep.m differ diff --git a/mycmap.cpt b/mycmap.cpt old mode 100755 new mode 100644 index ce38f5b..b6cf936 --- a/mycmap.cpt +++ b/mycmap.cpt @@ -1 +1,2 @@ --1 205 165 0 0 0 0 0 0 255 255 255 1 0 0 0 +-1 205 165 0 0 0 0 0 +0 255 255 255 1 0 0 0 diff --git a/mycmap1.cpt b/mycmap1.cpt old mode 100755 new mode 100644 diff --git a/mycmapCLR.cpt b/mycmapCLR.cpt old mode 100755 new mode 100644 index ce38f5b..b6cf936 --- a/mycmapCLR.cpt +++ b/mycmapCLR.cpt @@ -1 +1,2 @@ --1 205 165 0 0 0 0 0 0 255 255 255 1 0 0 0 +-1 205 165 0 0 0 0 0 +0 255 255 255 1 0 0 0 diff --git a/pcolorCENTER.m b/pcolorCENTER.m new file mode 100644 index 0000000..9d462ba --- /dev/null +++ b/pcolorCENTER.m @@ -0,0 +1,17 @@ +function pcolorCENTER(x,y,z,dx); +[N,M]=size(x); + +X=ones(N*2,M); +Y=ones(N*2,M); +Z=ones(N*2,M); + +i=[1:N]; +%for i=1:N + X(1+2*(i-1),:)=x(i,:); + X(2+2*(i-1),:)=x(i,:)+dx; + Y(1+2*(i-1),:)=y(i,:); + Y(2+2*(i-1),:)=y(i,:); + Z(1+2*(i-1),:)=z(i,:); + Z(2+2*(i-1),:)=z(i,:); +%end +pcolor(X,Y,Z) \ No newline at end of file diff --git a/periodicY.m b/periodicY.m new file mode 100644 index 0000000..26620af --- /dev/null +++ b/periodicY.m @@ -0,0 +1,82 @@ +function [a,q]=periodicY(k,N,M,p); +%q is the translated cell +%CAN BE OPTIMIZED + +[row col]=ind2sub([N M],[1:N*M]'); + +%row = rem([1:N*M]'-1,N)+1; +%row=[1:N]'*ones(1,M);row=reshape(row,1,N*M)'; +%col = ([1:N*M]'-row)/N + 1; + +if k==-1; + q=p+k; + a=find(row(p)>1); +end; + +if k==1; + q=p+k; + a=find(row(p)0); + G=[1:N*M]'; + Q=G+k; + s1=find(col==M); + s2=sub2ind([N M],row(s1),1+0*col(s1)); + %s2 = row(s1) + (1+0*col(s1)-1)*N; + Q(s1)=G(s2); + q=Q(p); +end; + +if k==-N;a=find(p>0); + G=[1:N*M]'; + Q=G+k; + s1=find(col==1); + s2=sub2ind([N M],row(s1),M+0*col(s1)); + %s2 = row(s1) + (M+0*col(s1)-1)*N; + Q(s1)=G(s2); + q=Q(p); +end; + + + + +%%%%%%%%%%%%%%%ORIGINAL +% function [a,q]=periodicY(k,N,M,p); +% [row col]=ind2sub([N M],p); +% +% q=p+k;%the translated cell +% +% +% if k==-1;a=find(row>1);end; +% if k==1;a=find(row0); +% s1=find(col==M); +% s2=find(col==1); +% q(s1)=p(s2); +% end; +% if k==-N;a=find(p>0);s1=find(col==1);s2=find(col==M);q(s1)=p(s2);end; +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + + + + + + + +% if k==N;a=find(p>0); +% s1=find(col==M); +% s2=find(col==1); +% %q(s1)=p(s1-M*N+M+N-1); +% s2-(s1-length(p)+N-1) +% q(s1)=p(s1-length(p)+N-1); +% end; +% pause +% if k==-N;a=find(p>0);s1=find(col==1);s2=find(col==M);q(s1)=p(s2);end; + +%if k==N;a=find(p>0);s=find(col==M);col(s)=1;q=sub2ind(N,M,row,col);end; +%if k==-N;a=find(p>0);s1=find(col==1);s2=find(col==M);q(s1)=q(s2);end; +%if k==-N;a=find(col-1>0);end; diff --git a/ponddynamics.m b/ponddynamics.m new file mode 100644 index 0000000..f6d41c7 --- /dev/null +++ b/ponddynamics.m @@ -0,0 +1,4 @@ +function [A,d,AC]=ponddynamics(A,d,dx,dt,MF,m,n,dntwr); +% A=pondformation(A,dx,dt,MF); +% [A,d]=isolatedpondexpansion(A,d,m,n,dx,dt,MF); +% [A,AC]=drainpond(d,A,m,n,dx,dt,MF,dntwr); \ No newline at end of file diff --git a/pondformation.m b/pondformation.m new file mode 100644 index 0000000..8136066 --- /dev/null +++ b/pondformation.m @@ -0,0 +1,26 @@ +function [deltaY2,pondlost]=pondformation(A,dx,dt,Epondform,z,zpondcr,maxdpond,zsill,pondlost); + +zoriginal=z; + +p=find(A==1); %how many points can become a pond, how much normal cells are present. + +dz=max(z(p)-zpondcr,0);%base level of the pond +dz=min(maxdpond,dz);%maximum scour depth + +%a=find(dz>0);%only create new ponds when the scour can be made! +a=find(dz>0 & z(p)0 + %S(p(aa(i)))=1; + z(p(aa(i)))=z(p(aa(i)))-dz(aa(i)); + pondlost=pondlost+dz(aa(i)); + end +end + + +deltaY=zoriginal-z; +deltaY2=deltaY; \ No newline at end of file diff --git a/run_MM2D_example1createarticlefigure3.m b/run_MM2D_example1createarticlefigure3.m new file mode 100644 index 0000000..06eedcd --- /dev/null +++ b/run_MM2D_example1createarticlefigure3.m @@ -0,0 +1,260 @@ +clear; close all;clc + + +%NOTES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Qs is the total lateral transport, E is the vertical erosion flux + +%lateral boundary options +% 0 is closed (no flux if nothign specified. no-gradient if AW equal to 1 or -1) +% 1 is periodic + +%A +%0: not in the domain +%1: a normal cell +%2: the open sea boundary +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + +%%%%%%%%%%%%%%PARAMETERS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +gR=[2.5 2.5 2.5 2.5 5 5 5 5 7.5 7.5 7.5 7.5 10 10 10 10]; +gC=[60 30 15 7.5 60 30 15 7.5 60 30 15 7.5 60 30 15 7.5]; + +figure('units','normalized','outerposition',[0 0 1 1]) %set(gca,'Color','k') +for contatore=1:16 + +rng(3)%to get always the same random numbers. you can pick any number in there + +P.g=9.81; %gravity [m/s2] +P.rho=1030; %water density [kg/m3] +P.rhos=2650; %sediment density (quartz-mica) [kg/m3] +P.ss=(P.rhos-P.rho)/P.rho; %relative density +P.kro=0.1;%[m] minimum water depth for hydrodynamics. NEEDS TO BE SMALLER THAN hwSea_lim +P.DiffSmud=1; %[-]coefficient for tidal dispersion [-]. 1 DO NOT CHANGE +P.DoMUD=1;%base diffusivity of suspedned mud [m2/s]. Process not related to tides (e.g. wind and waves, other ocean circulation) + +%Sea level rise +P.RSLR=gR(contatore)/1000/365; %from mm/yr to m/day (the time unit is the day!) + +%Tide +P.Ttide=12.5/24; %tidal period [day] +P.Trange=0.7; %Tidal Trange [m] +P.TrangeVEG=P.Trange;%tidal range for vegetation [m]. Generally same of tidal range + +%Wind for sea waves +P.wind=7;%reference wind speed [m/s] + +%Edge erosion +P.aw=0.3/365; %wave edge erodability m/yr/W/m2 +P.maxedgeheight=2; +P.fox=0;%fraction of edge eroded material that is oxidized. + +%Wind waves and swell numerics +P.hwSwelltransport_lim=1; +P.hwSea_lim=0.2;%0.5; %limiter water deth of sea waves %THIS IS USED TO FIND THE "EDGE"%NEEDS TO BE LARGER THAN KO!!!! + +%SSC at the sea boundary +P.co2=gC(contatore)/1000;%40/1000; %Sea boundary SSC for mud [g/l] +%P.co2=50/1000;%40/1000; %Sea boundary SSC for mud [g/l] + +%Manning coeffinent unvegeated (same for sand and mud) +P.Cb=0.02; + +%Mud +P.d50_2=0.02/1000/1000;%mud grain size [m] +P.ws2=0.2/1000;%0.2/1000;% m/s +P.por2=0.7;P.rbulk2=P.rhos*(1-P.por2);%P.por2=0.7 + +%Mud parameters +P.me=0.1*10^-4*24*3600; %per day!!! +P.taucr=0.2; +P.crMARSH=0.1/365;%creep coefficient vegetated +P.crMUD=3.65/365;%creep coeffcinet +P.alphaMUD=0.25; %coefficient for bedload downslope of mud. added April 2019. Similar to P.alphaSAND + +%Vegetation parameters +P.dBlo=0;%-0.2;%-0.2;%0;%-(0.237*P.TrangeVEG-0.092); +P.dBup=P.TrangeVEG/2;%-0.2; +P.Cv=0.1;%Manning for vegetated ara +P.wsB=1/1000;%Mud Settling velocity for vegetated ara +P.taucrVEG=0.5;%Critical sheak stress for vegetated areas + +%Organic accretion by vegetation +P.AccreteOrganic=1; +P.Korg=6/1000/365;% [mm/yr] + +%ON/OFF processes +P.VEGETATION=1;%vegeation resistance and settling (DOES NOT controll organic accretion) +P.computeSeaWaves=1; +P.computeEdgeErosionSea=1; +P.computetide=1; + +%Various boundary condtions +P.periodic=0; +P.imposeseaboundarydepthmorphoALL=1; %to use when a channel mouth is at a boundary + +%Pond dynamics +P.calculateponddynamics=1; + P.Epondform=0.4*10^-5;%4*10^-4/10;%probabiliy of new pond formation per year (area/area) + P.zpondcr=-0.2;%P.Trange/4;%base of new pond formation with respect to MSL + P.minponddepth=0.1;%1; %minimum depth to define a pond after you identified the "lakes" + P.maxdpond=max(0.2,max(P.minponddepth*1.1,0.15*P.Trange));%0.5;%0.5;%maximum depth scour of a new pond + P.zntwrk=(P.Trange/2)*0.5;%(P.Trange/2)*0.2;%P.Trange/2*0.9;%P.Trange/2-0.3;%0.3; %depth above msl that defines the channel network. the smaller the harder to drain! + P.aPEXP=0.015*10;%isolated pond expansion rate m/yr + P.ponddeeprate=0.003;%m/yr + +%Global numerical limiters +limitdeltaz=2;%[m] +limitmaxup=1;%[m] + +%Time parameters +tmax=200;% number of time steps +tINT=1;%how many time steps you want to do the plot (if 1 you plot every time step). Does not affect the computation + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%time series +numberserie=10000; +dtOserie=ones(numberserie,1)*365; + +time=cumsum(dtOserie)/365; %converted to years, just to plot. Does not affect the computation +time=[time(2:end);time(end)]; + +%SeaWave direction +angleWINDserie=rand(numberserie,1)*360; %every time step a random direction +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + + +%%%%%%%%%%%%%%%Geometry Initilization%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%[N,M,dx,A,z,Active,x,y,msl]=initializegeometry(P);%use this to create a new marsh +load REFERENCEMARSHr07; %use this to load an existign marsh + +%savegeometry(N,M,dx,A,z,Active,x,y,msl,'name_SAVEYOURCURRENTCONFIGURATION');%execute this once to save you current configuration + +%Store value for mass balance check%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +sumzIN=sum(z(A==1)); +FLXz=zeros(4,1); +KBTOT=0; +zOX=0; +pondloss=0; + +IO.z=z; +IO.msl=msl; +IO.Active=Active; +fIO.FLXz=FLXz; +fIO.pondloss=pondloss; +fIO.KBTOT=KBTOT; +fIO.zOX=zOX; + +zbedo=z-msl; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + +makevideo=0;%put one if you want to create a video +%v=VideoWriter('nameofvideo','Motion JPEG AVI'); + +%%%%%%%%%%%%%%%%%%%%%%MAIN LOOP%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if makevideo==1;open(v);end +s=0;step=0;tic; +for t=1:tmax; %iteration over the tmax time stpes + +%SeaWave direction +angleWIND=angleWINDserie(t);%rand(1)*360; %every time step a random direction +%Length of time step +dtO=dtOserie(t); + +if t==1;dto=0.00001;else;dto=dtO;end + +dti=0;dt=dto; +while dtilimitdeltaz | maxup>limitmaxup + if firstattemp==1;else;dt=dt/2*min(limitdeltaz/maxdeltaz,limitmaxup/maxup);end;firstattemp=0; + if t<=2;dt=min(0.5*365,dt);end + [IOtemp,fIOtemp,maxdeltaz,maxup,PLT]=mainevolutionstep(A,P,dx,dt,IO,fIO,angleWIND,t); + step=step+1; %this is how many time you called the function mainevolution step + end + + %the partial updating step was succefull! Keep going + IO=IOtemp; + fIO=fIOtemp; + dti=dti+dt;%how much you moved forward + dt=min(dt*2,max(0,dto-dti));%the remaining time in the time step +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + + + +%%%PLOT%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if mod(t,tINT)==0;s=s+1; +%read the variables +names = fieldnames(IO); +for i=1:length(names);eval([names{i} '=IO.' names{i} ';' ]);end + +%read the fluxes +names = fieldnames(fIO); +for i=1:length(names);eval([names{i} '=fIO.' names{i} ';' ]);end + +%read the plot +names = fieldnames(PLT); +for i=1:length(names);eval([names{i} '=PLT.' names{i} ';' ]);end + +zbed=z-msl; + +ax1 = subplot(4,4,contatore); +IM=imagesc(x,y,zbed');axis equal;set(IM,'alphadata',~(A'==0));set(gca,'YDir','normal'); +cmp=demcmap([-3 P.Trange/2],256); %-3 1 P.Trange/2 +colormap(ax1,cmp) +caxis([-3 P.Trange/2]); +title(strcat(num2str(time(t)),' years ',num2str(step))) + + + + +%THIS IS JUST TO MAKE THE SEDIMENT BUDGET AND CHECK THAT THE SEDIMENT +%CONSERVE SEDIMENT!!! NOT NECESSARY FOR THE COMPUTATION +QseaTide=FLXz(2); +QmouthTide=FLXz(4); +sumFLUX2=-QseaTide*dx-QmouthTide*dx; +sumz=sum(z(A==1)); +%NOTE: Thsi is the equivalent volumetric flux, not the mass flux +checksum= [(sumzIN-sumz)+sumFLUX2/dx^2]-pondloss+KBTOT-zOX ; +if abs(checksum(1))>0.1 ;checksum,pause;end% if this is not zero, the code has an error somewhere + + +%this is to make the video +if makevideo==1;V=getframe(figure(1));writeVideo(v,V);end + + +pause(0.1) +end%end of plot +end; %end of panel plotting +end;%of of repetitie cycle + + + +%this is to make the video +if makevideo==1;close(v);end diff --git a/run_MM2D_example2startfromslopingdomain.m b/run_MM2D_example2startfromslopingdomain.m new file mode 100644 index 0000000..b7de21c --- /dev/null +++ b/run_MM2D_example2startfromslopingdomain.m @@ -0,0 +1,254 @@ +clear; close all;clc + + +%NOTES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Qs is the total lateral transport, E is the vertical erosion flux + +%lateral boundary options +% 0 is closed (no flux if nothign specified. no-gradient if AW equal to 1 or -1) +% 1 is periodic + +%A +%0: not in the domain +%1: a normal cell +%2: the open sea boundary +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + +%%%%%%%%%%%%%%PARAMETERS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +figure('units','normalized','outerposition',[0 0 1 1]) %set(gca,'Color','k') + + +rng(3)%to get always the same random numbers. you can pick any number in there + +P.g=9.81; %gravity [m/s2] +P.rho=1030; %water density [kg/m3] +P.rhos=2650; %sediment density (quartz-mica) [kg/m3] +P.ss=(P.rhos-P.rho)/P.rho; %relative density +P.kro=0.1;%[m] minimum water depth for hydrodynamics. NEEDS TO BE SMALLER THAN hwSea_lim +P.DiffSmud=1; %[-]coefficient for tidal dispersion [-]. 1 DO NOT CHANGE +P.DoMUD=1;%base diffusivity of suspedned mud [m2/s]. Process not related to tides (e.g. wind and waves, other ocean circulation) + +%Sea level rise +P.RSLR=2.5/1000/365; %from mm/yr to m/day (the time unit is the day!) + +%Tide +P.Ttide=12.5/24; %tidal period [day] +P.Trange=0.7; %Tidal Trange [m] +P.TrangeVEG=P.Trange;%tidal range for vegetation [m]. Generally same of tidal range + +%Wind for sea waves +P.wind=7;%reference wind speed [m/s] + +%Edge erosion +P.aw=0.3/365; %wave edge erodability m/yr/W/m2 +P.maxedgeheight=2; +P.fox=0;%fraction of edge eroded material that is oxidized. + +%Wind waves and swell numerics +P.hwSwelltransport_lim=1; +P.hwSea_lim=0.2;%0.5; %limiter water deth of sea waves %THIS IS USED TO FIND THE "EDGE"%NEEDS TO BE LARGER THAN KO!!!! + +%SSC at the sea boundary +P.co2=60/1000;%40/1000; %Sea boundary SSC for mud [g/l] + +%Manning coeffinent unvegeated (same for sand and mud) +P.Cb=0.02; + +%Mud +P.d50_2=0.02/1000/1000;%mud grain size [m] +P.ws2=0.2/1000;%0.2/1000;% m/s +P.por2=0.7;P.rbulk2=P.rhos*(1-P.por2);%P.por2=0.7 + +%Mud parameters +P.me=0.1*10^-4*24*3600; %per day!!! +P.taucr=0.2; +P.crMARSH=0.1/365;%creep coefficient vegetated +P.crMUD=3.65/365;%creep coeffcinet +P.alphaMUD=0.25; %coefficient for bedload downslope of mud. added April 2019. Similar to P.alphaSAND + +%Vegetation parameters +P.dBlo=0;%-0.2;%-0.2;%0;%-(0.237*P.TrangeVEG-0.092); +P.dBup=P.TrangeVEG/2;%-0.2; +P.Cv=0.1;%Manning for vegetated ara +P.wsB=1/1000;%Mud Settling velocity for vegetated ara +P.taucrVEG=0.5;%Critical sheak stress for vegetated areas + +%Organic accretion by vegetation +P.AccreteOrganic=1; +P.Korg=6/1000/365;% [mm/yr] + +%ON/OFF processes +P.VEGETATION=1;%vegeation resistance and settling (DOES NOT controll organic accretion) +P.computeSeaWaves=1; +P.computeEdgeErosionSea=1; +P.computetide=1; + +%Various boundary condtions +P.periodic=0; +P.imposeseaboundarydepthmorphoALL=1; %to use when a channel mouth is at a boundary + +%Pond dynamics +P.calculateponddynamics=1; + P.Epondform=0.4*10^-5;%4*10^-4/10;%probabiliy of new pond formation per year (area/area) + P.zpondcr=-0.2;%P.Trange/4;%base of new pond formation with respect to MSL + P.minponddepth=0.1;%1; %minimum depth to define a pond after you identified the "lakes" + P.maxdpond=max(0.2,max(P.minponddepth*1.1,0.15*P.Trange));%0.5;%0.5;%maximum depth scour of a new pond + P.zntwrk=(P.Trange/2)*0.5;%(P.Trange/2)*0.2;%P.Trange/2*0.9;%P.Trange/2-0.3;%0.3; %depth above msl that defines the channel network. the smaller the harder to drain! + P.aPEXP=0.015*10;%isolated pond expansion rate m/yr + P.ponddeeprate=0.003;%m/yr + +%Global numerical limiters +limitdeltaz=2;%[m] +limitmaxup=1;%[m] + +%Time parameters +tmax1000;% number of time steps +tINT=1;%how many time steps you want to do the plot (if 1 you plot every time step). Does not affect the computation + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%time series +numberserie=10000; +dtOserie=ones(numberserie,1)*365; + +time=cumsum(dtOserie)/365; %converted to years, just to plot. Does not affect the computation +time=[time(2:end);time(end)]; + +%SeaWave direction +angleWINDserie=rand(numberserie,1)*360; %every time step a random direction +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + + +%%%%%%%%%%%%%%%Geometry Initilization%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +[N,M,dx,A,z,Active,x,y,msl]=initializegeometry(P);%use this to create a new marsh +%load REFERENCEMARSHr07; %use this to load an existign marsh + +%savegeometry(N,M,dx,A,z,Active,x,y,msl,'name_SAVEYOURCURRENTCONFIGURATION');%execute this once to save you current configuration + +%Store value for mass balance check%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +sumzIN=sum(z(A==1)); +FLXz=zeros(4,1); +KBTOT=0; +zOX=0; +pondloss=0; + +IO.z=z; +IO.msl=msl; +IO.Active=Active; +fIO.FLXz=FLXz; +fIO.pondloss=pondloss; +fIO.KBTOT=KBTOT; +fIO.zOX=zOX; + +zbedo=z-msl; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + +makevideo=0;%put one if you want to create a video +%v=VideoWriter('nameofvideo','Motion JPEG AVI'); + +%%%%%%%%%%%%%%%%%%%%%%MAIN LOOP%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if makevideo==1;open(v);end +s=0;step=0;tic; +for t=1:tmax; %iteration over the tmax time stpes + +%SeaWave direction +angleWIND=angleWINDserie(t);%rand(1)*360; %every time step a random direction +%Length of time step +dtO=dtOserie(t); + +if t==1;dto=0.00001;else;dto=dtO;end + +dti=0;dt=dto; +while dtilimitdeltaz | maxup>limitmaxup + if firstattemp==1;else;dt=dt/2*min(limitdeltaz/maxdeltaz,limitmaxup/maxup);end;firstattemp=0; + if t<=2;dt=min(0.5*365,dt);end + [IOtemp,fIOtemp,maxdeltaz,maxup,PLT]=mainevolutionstep(A,P,dx,dt,IO,fIO,angleWIND,t); + step=step+1; %this is how many time you called the function mainevolution step + end + + %the partial updating step was succefull! Keep going + IO=IOtemp; + fIO=fIOtemp; + dti=dti+dt;%how much you moved forward + dt=min(dt*2,max(0,dto-dti));%the remaining time in the time step +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + + + +%%%PLOT%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if mod(t,tINT)==0;s=s+1; +%read the variables +names = fieldnames(IO); +for i=1:length(names);eval([names{i} '=IO.' names{i} ';' ]);end + +%read the fluxes +names = fieldnames(fIO); +for i=1:length(names);eval([names{i} '=fIO.' names{i} ';' ]);end + +%read the plot +names = fieldnames(PLT); +for i=1:length(names);eval([names{i} '=PLT.' names{i} ';' ]);end + +zbed=z-msl; + +IM=imagesc(x,y,zbed');axis equal;set(IM,'alphadata',~(A'==0));set(gca,'YDir','normal'); +cmp=demcmap([-3 P.Trange/2],256); %-3 1 P.Trange/2 +colormap(cmp) +caxis([-3 P.Trange/2]); +title(strcat(num2str(time(t)),' years ',num2str(step))) + + + + +%THIS IS JUST TO MAKE THE SEDIMENT BUDGET AND CHECK THAT THE SEDIMENT +%CONSERVE SEDIMENT!!! NOT NECESSARY FOR THE COMPUTATION +QseaTide=FLXz(2); +QmouthTide=FLXz(4); +sumFLUX2=-QseaTide*dx-QmouthTide*dx; +sumz=sum(z(A==1)); +%NOTE: Thsi is the equivalent volumetric flux, not the mass flux +checksum= [(sumzIN-sumz)+sumFLUX2/dx^2]-pondloss+KBTOT-zOX ; +if abs(checksum(1))>0.1 ;checksum,pause;end% if this is not zero, the code has an error somewhere + + +%this is to make the video +if makevideo==1;V=getframe(figure(1));writeVideo(v,V);end + + +pause(0.1) +end%end of plot +end; %end of panel plotting + + + +%this is to make the video +if makevideo==1;close(v);end diff --git a/savegeometry.m b/savegeometry.m new file mode 100644 index 0000000..4072bb9 --- /dev/null +++ b/savegeometry.m @@ -0,0 +1,3 @@ +function savegeometry(N,M,dx,A,z,Active,x,y,msl,name); + +save(name) diff --git a/savegeometry_3sediments.m b/savegeometry_3sediments.m new file mode 100644 index 0000000..70ec79b --- /dev/null +++ b/savegeometry_3sediments.m @@ -0,0 +1,3 @@ +function savegeometry(N,M,dx,A,AW,Yb,plyr,zb,z,Y1,Y2,Y3,flyr1,flyr2,flyr3,flyrb1,flyrb2,flyrb3,Active,x,y,msl,SPCLcell,name); + +save(name) diff --git a/seaboundaryNeumanbedelevationALLBOUNDARY.m b/seaboundaryNeumanbedelevationALLBOUNDARY.m old mode 100755 new mode 100644 index c3b1f86..d540f88 --- a/seaboundaryNeumanbedelevationALLBOUNDARY.m +++ b/seaboundaryNeumanbedelevationALLBOUNDARY.m @@ -1,4 +1,4 @@ -function [z]=seaboundaryNeumanbedelevationALLBOUNDARY(A,z) +function [Y2]=seaboundaryNeumanbedelevationALLBOUNDARY(A,Y2) [N,M]=size(A); p=find(A==2);%exclude the NOLAND CELLS (A==0) @@ -9,7 +9,9 @@ [a,q]=excludeboundarycell(k,N,M,p); a=a(A(q(a))==1);%only inclued the cells in whcih you can creep to -z(p(a))=z(q(a)); + +Y2(p(a))=Y2(q(a)); + end diff --git a/sedtran.m b/sedtran.m old mode 100755 new mode 100644 index 8adf8de..d832b13 --- a/sedtran.m +++ b/sedtran.m @@ -1,6 +1,17 @@ -function [EmD,SSM,FLX]=sedtran(d,A,DiffS,h,ho,E,WS,dx,dt,rbulk,co,Ux,Uy,FLX,fTide,Ttide,kro); +function [EmD,SSM,FLX]=sedtran(numeric,d,A,DoMUD,DiffS,h,ho,E,WS,dx,dt,rbulk,co,Ux,Uy,FLX,fTide,Ttide,URx,URy,UresX,UresY,periodic,computeriver,computetide,residualcurrents,kro,MUD,coMUD,Qsmouthi); -A(ho<=0)=0; %eliminate the cells in which the water depth is too small + +QmouthRiver=FLX(1); +QseaTide=FLX(2); +QseaRiver=FLX(3); +QmouthTide=FLX(4); +%consider the pond cells as A==0 +%A(A==3)=1; %but do not update it! +%A(A==22)=2; + +A(ho<=0)=0; %ATTENTION: eliminate the cells in which the water depth is too small. They cannot evolve at all!!!! They will not erode nor accrete!!!!!!!!!!!!!!!!!!!!!! + +%h=max(0.1,h); %rivermouthfront p = find(A>0);%exclude the NOLAND CELLS (A==0) @@ -12,23 +23,35 @@ %boundary conditions imposed SSC a=find(A==2);rhs(G(a))=co.*h(a).*fTide(a); -Dxx=(DiffS*Ttide/2*(abs(Ux.*Ux))*(24*3600).^2)/(dx^2).*h;%.*(ho>kro);%.*(hgross>0.01);%% the h is not the coefficient for diffusion, it the h in the mass balance eq. -Dyy=(DiffS*Ttide/2*(abs(Uy.*Uy))*(24*3600).^2)/(dx^2).*h;%.*(ho>kro);%.*(hgross>0.01); +%iver and mud (if mud, you impose the SSC at the inlet) +if MUD==1; +a=find(A==10);rhs(G(a))=coMUD.*h(a).*fTide(a); +end + -%the factor 24*3600 is used to convert the Ux and Uy from m/s to m/day +Dxx=(DoMUD*24*3600 +DiffS*Ttide/2.*(abs(Ux.*Ux))*(24*3600).^2)/(dx^2).*h;%.*h;%.*(ho>kro);%.*(hgross>0.01);%% the h is not the coefficient for diffusion, is the h in the mass balance eq. +Dyy=(DoMUD*24*3600 +DiffS*Ttide/2.*(abs(Uy.*Uy))*(24*3600).^2)/(dx^2).*h;%.*h;%.*(ho>kro);%.*(hgross>0.01); [row col]=ind2sub(size(A),p); for k = [N -1 1 -N] %go over the 4 directions for the gradients %avoid to the the cells out of the domain (risk to make it periodic...) +if periodic==0 [a,q]=excludeboundarycell(k,N,M,p); +elseif periodic==1; +[a,q]=periodicY(k,N,M,p); %for the long-shore +end a=a(A(q(a))==1 | A(q(a))==2 | A(q(a))==10);%exlcude the translated cell that are NOLAND cells %go over the extradiagonal direction for each gradient direction if (k==N | k==-N);D=Dyy;else;D=Dxx;end; -DD=(D(p(a))+D(q(a)))/2; +if numeric==1 +DD=(D(p(a))+D(q(a)))/2;%THE ORIGINAL FOR MARSH fa piu' accumulo sui lati dei canali, cioe canali piou stretti +elseif numeric==0 +DD=min(D(p(a)),D(q(a)));%fa meno levees on the sides. sopratutto con la sand +end Fin=0*DD;Fin(A(p(a))==1)=1; % (A(q(a))==1 | A(q(a))==10)=1; to conserve the mass a the river mouth = no input Fout=0*DD;Fout(A(q(a))==1)=1; %needed not to affect the b.c. -> Do not choose 2 and 1p @@ -36,6 +59,27 @@ %tidal dispersion component value=DD./h(p(a))./fTide(p(a)); +%river flow component +if computeriver==1 +if (k==N);UR=URy(p(a));up=find(UR>0);F=UR(up);end %East-west +if (k==-N);UR=URy(p(a));up=find(UR<0);F=-UR(up);end +if (k==1);UR=URx(p(a));up=find(UR>0);F=UR(up);end %North-south +if (k==-1);UR=URx(p(a));up=find(UR<0);F=-UR(up);end +value(up)=value(up)+F*3600*24/dx; +end + +if residualcurrents==1;%tidal residual currents ansd transport. + % (I imposed no residual currents are the opne boundary to avoid + % calcuitating the fluxes to get the mass balance at 100%) +if (k==N);UR=UresY(p(a));up=find(UR>0);F=UR(up);end +if (k==-N);UR=UresY(p(a));up=find(UR<0);F=-UR(up);end +if (k==1);UR=UresX(p(a));up=find(UR>0);F=UR(up);end +if (k==-1);UR=UresX(p(a));up=find(UR<0);F=-UR(up);end +value(up)=value(up)+F*3600*24/dx; +end + + + S(p(a))=S(p(a))+value.*Fin; %exit from that cell i=[i;G(q(a))]; j=[j;G(p(a))]; s=[s;-value.*Fout]; %gain from the neigborh cell @@ -49,6 +93,12 @@ settling(a)=0;%do not settle in the b.c. (to not change the SSC) S(p(a))=1;%to impose the b.c. +%river boundary +if MUD==1; %if mud, handlew this as an imposed SSC +a=find(A(p)==10);%find the co b.c. +settling(a)=0;%do not settle in the b.c. (to not change the SSC) +S(p(a))=1;%to impose the b.c. +end i=[i;G(p)]; j=[j;G(p)]; s=[s;S(p)+settling]; ds2 = sparse(i,j,s);%solve the matrix inversion @@ -67,28 +117,87 @@ + + + + + + + + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%OUTPUT FOR SSM BALANCE. DOES NOT AFFECT COMPUTATION!!! p = find(A==1);[row col]=ind2sub(size(A),p); -%sea boundary -Q=0; +%Sea boundary%%%%%%%%%%%%%%% +Q=0;QR=0; for k = [N -1 1 -N] +if periodic==0 [a,q]=excludeboundarycell(k,N,M,p); +elseif periodic==1; +[a,q]=periodicY(k,N,M,p); %for the long-shore +end %only get the cell that discharge into a b.c. A==2 a=a(A(q(a))==2); %for each gradient direction if (k==N | k==-N);D=Dyy;else;D=Dxx;end; -DD=(D(p(a))+D(q(a)))/2; +if numeric==1 +DD=(D(p(a))+D(q(a)))/2;%THE ORIGINAL FOR MARSH fa piu' accumulo sui lati dei canali, cioe canali piou stretti +elseif numeric==0 +DD=min(D(p(a)),D(q(a)));%fa meno levees on the sides. sopratutto con la sand +end +%Tide Q=Q+sum(DD.*(SSM(p(a))./h(p(a))./fTide(p(a))-SSM(q(a))./h(q(a))./fTide(q(a)))); %exit from that cell +end + +QseaTide=QseaTide+dx*dt*Q/rbulk; %(note the the divided dx is for the gradient, not for the cell width!) +QseaRiver=QseaRiver+dt*3600*24*QR/rbulk; +%%%%%%%%%%%%%%%%%%% + + +%River mouth%%%%%%%%%%%%%%%% +Q=0; +for k = [N -1 1 -N] +if periodic==0 +[a,q]=excludeboundarycell(k,N,M,p); +elseif periodic==1; +[a,q]=periodicY(k,N,M,p); %for the long-shore +end +%exlcude the translated cell that are the b.c. +a=a(A(q(a))==10); +if (k==N | k==-N);D=Dyy;else;D=Dxx;end +if numeric==1%mud +DD=(D(p(a))+D(q(a)))/2;%THE ORIGINAL FOR MARSH fa piu' accumulo sui lati dei canali, cioe canali piou stretti +elseif numeric==0%sand +DD=min(D(p(a)),D(q(a)));%fa meno levees on the sides. sopratutto con la sand end -FLX=FLX+dx*dt*Q/rbulk; %(note the the divided dx is for the gradient, not for the cell width!) +%Tide +Q=Q+sum(DD.*(SSM(p(a))./h(p(a))./fTide(p(a))-SSM(q(a))./h(q(a))./fTide(q(a)))); %exit from that cell +%River +if (k==N | k==-N);QR=QR+sum((SSM(p(a)).*URy(p(a))))*sign(k);end +if (k==1 | k==-1);QR=QR+sum((SSM(p(a)).*URx(p(a))))*sign(k);end +end + +QmouthTide=QmouthTide+dx*dt*Q/rbulk; +QmouthRiver=QmouthRiver+dt*3600*24*QR/rbulk; %%%%%%%%%%%%%%%%%%% +%SUM OF THE RIVER INPUT (IMPOSED FROM THE OUTISE USING Qsmouthi, which is cl +QmouthRiver=FLX(1)+Qsmouthi*dt*24*3600/rbulk; +FLX=[QmouthRiver;QseaTide;QseaRiver;QmouthTide]; + + + + diff --git a/sumSedcolum.m b/sumSedcolum.m new file mode 100644 index 0000000..7d28a48 --- /dev/null +++ b/sumSedcolum.m @@ -0,0 +1,3 @@ +function sumY=sumSedcolum(Yb,flyrb,flyr,dlyr,Yi) + +sumY=[Yb+Yi]; \ No newline at end of file diff --git a/tidalFLOW.m b/tidalFLOW.m old mode 100755 new mode 100644 index 373ad6c..d3a3b99 --- a/tidalFLOW.m +++ b/tidalFLOW.m @@ -1,8 +1,13 @@ -function [U,Ux,Uy,q,P]=flowBasin(A,MANN,h,d,dx,DH,T,kro); +function [U,Ux,Uy,q,P]=flowBasin(A,MANN,h,ho,d,dx,DH,T,periodic,kro); Uo=1; +A(A==22)=1; %this behaves as normal flow %but do not update A! +%consider the pond cells as A==0 +A(A==3)=1; %the isoalted pond behaves as normal cell (btu different depth...) %but do not update A! + MANN(isnan(MANN))=0.1;% csi=h.^(1/3)./MANN.^2./Uo*24*3600; + D=csi.*h.^2/(dx^2); G=0*d;a=find(A~=0);NN=length(a);G(a)=[1:NN]; @@ -11,7 +16,7 @@ [N,M] = size(G);i=[];j=[];s=[]; %boundary conditions imposed water level -a=find(A==2); +a=find(A==2 | A==21); i=[i;G(a)]; j=[j;G(a)]; s=[s;ones(size(a))];rhs(G(a))=0;%water level zero S=0*G; @@ -20,12 +25,17 @@ for k = [N -1 1 -N] %the translated cells +if periodic==0 [a,q]=excludeboundarycell(k,N,M,p); +elseif periodic==1; +[a,q]=periodicY(k,N,M,p); %for the long-shore +end a=a(A(q(a))>0);%exlcude the translated cell that are NOLAND cells DD=(D(p(a))+D(q(a)))/2;%.*(fM(p(a))+fM(q(a)))/2; %THA BEST!!!! BESTA! WITH THIS MORE STABLE + S(p(a))=S(p(a))+DD; %exit from that cell i=[i;G(q(a))]; j=[j;G(p(a))]; s=[s;-DD]; %gain from the neigborh cell end @@ -39,28 +49,37 @@ P(A==2)=0; %need when swtinching q and p + D=D./h*dx; Ux=0*A;Uy=0*A; U1=0*A;Um1=0*A;UN=0*A;UmN=0*A; p = find(A==1 | A==10 | A==2);[row col]=ind2sub(size(A),p); for k = [N -1 1 -N] %the translated cell +if periodic==0 [a,q]=excludeboundarycell(k,N,M,p); - +elseif periodic==1; +[a,q]=periodicY(k,N,M,p); %for the long-shore +end a=a(A(q(a))>0);%exlcude the translated cell that are NOLAND cells +%DD=(D(p(a))+D(q(a)))/2; DD=min(D(p(a)),D(q(a)));%.*(fM(p(a))+fM(q(a)))/2; MEGLIO -if (k==1); U1(p(a))=U1(p(a))+sign(k)*(P(p(a))-P(q(a))).*DD; -elseif (k==-1); Um1(p(a))=Um1(p(a))+sign(k)*(P(p(a))-P(q(a))).*DD; -elseif (k==N); UN(p(a))=UN(p(a))+sign(k)*(P(p(a))-P(q(a))).*DD; -elseif (k==-N); UmN(p(a))=UmN(p(a))+sign(k)*(P(p(a))-P(q(a))).*DD; + +if (k==1); U1(p(a))=U1(p(a))+sign(k)*(P(p(a))-P(q(a))).*DD;%./h(p(a))*dx; +elseif (k==-1); Um1(p(a))=Um1(p(a))+sign(k)*(P(p(a))-P(q(a))).*DD;%./h(p(a))*dx; +elseif (k==N); UN(p(a))=UN(p(a))+sign(k)*(P(p(a))-P(q(a))).*DD;%./h(p(a))*dx; +elseif (k==-N); UmN(p(a))=UmN(p(a))+sign(k)*(P(p(a))-P(q(a))).*DD;%./h(p(a))*dx; end end -Uy=max(abs(U1),abs(Um1)).*sign(U1+Um1); -Ux=max(abs(UN),abs(UmN)).*sign(UN+UmN); + +Uy=max(abs(U1),abs(Um1)).*sign(U1+Um1);%MEGLIO +Ux=max(abs(UN),abs(UmN)).*sign(UN+UmN);%MEGLIO + + U=sqrt(Ux.^2+Uy.^2); q=U.*h; @@ -72,3 +91,6 @@ + + + diff --git a/totalsedimenterosionMUDsine.m b/totalsedimenterosionMUDsine.m old mode 100755 new mode 100644 index 8df14c8..bc891e0 --- a/totalsedimenterosionMUDsine.m +++ b/totalsedimenterosionMUDsine.m @@ -1,15 +1,9 @@ -function [E]=totalsedimenterosionMUDsine(U,MANN,VEG,fTide,taucr,tcrgradeint,taucrVEG,me,h,lev,TrangeVEG); +function [E,Etide]=totalsedimenterosionMUDsine(U,MANN,VEG,fTide,UR,Uwave_sea,Uwave,Tp_sea,Tp_swell,fMFswell,fMFsea,fMFriver,taucr,taucrVEG,me,h,lev,TrangeVEG,computeSeaWaves,computeSwellwave,computeRiver); fUpeak=pi/2; taucro=U*0+taucr; taucro(VEG==1)=taucrVEG; -%increase tcr with depth (asusming an existing vertical distribution. -%USE with cauton, only ok for simulation of small marsh domain -xi=-lev-TrangeVEG/2;xi(xi<0)=0; -taucro=taucro+xi*tcrgradeint; -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %%%%%%%%%%%%%%%%%%%%%tidal current erosion ncyc=10; @@ -19,4 +13,25 @@ tauC=1030*9.81*MANN.^2.*h.^(-1/3).*Utide.^2; E=E+1/(ncyc+1)*me.*(sqrt(1+(tauC./taucro).^2)-1); end +Etide=E; + + +%Sea waves erosion +if computeSeaWaves==1 +Uwave_sea=max(0,Uwave_sea); + +ko=0.1/1000*3; +aw=Tp_sea.*Uwave_sea/(2*pi); +fw=0.00251*exp(5.21*(aw/ko).^-0.19);fw(aw/ko SIMPLE, s +% VO = TSMOVAVG(VI, 'e', TIMEPER, DIM) => EXPONENTIAL, e +% VO = TSMOVAVG(VI, 't', NUMPER, DIM) => TRIANGULAR, t +% VO = TSMOVAVG(VI, 'w', WEIGHTS, DIM) => WEIGHTED, w +% VO = TSMOVAVG(VI, 'm', NUMPER, DIM) => MODIFIED, m +% +% NOTE: The Simple and Weighted moving averages have a new calling +% syntax. This is different from previous versions of this function. +% The simple moving average does not require a LEAD input and the +% weighted moving average does not require a PIVOT point. [END NOTE] +% +% VI is assumed to be a row vector or row-oriented matrix. DIM is an +% optional input, and if it is not included as an input, the default +% value 2 is assumed. DIM = 2 means that each row is a variable and +% each column is an observation (row-oriented matrix). If a column- +% oriented matrix is supplied for VI (i.e. each column is a variable +% and each row is an observation) then a DIM must be included and set +% to the value DIM = 1. VO is identical in form to VI. +% +% Simple Moving Average: +% LAG is the parameter indicating the number of previous data points to +% be used in conjunction with the current data point when calculating +% the moving average. LAG can also be thought of as the window size or +% number of periods of the moving average. +% +% Exponential Moving Average: +% Exponential moving average is a weighted moving average, where TIMEPER is +% the time period of the exponential moving average. Exponential moving +% averages reduce the lag by applying more weight to recent prices. For +% example, a 10 period exponential moving average weights the most recent +% price by 18.18%. +% +% Exponential Percentage = 2/(TIMEPER + 1) or 2/(WINDOW_SIZE + 1) +% +% Triangular Moving Average: +% Triangular moving average is a double-smoothing of the data. The first +% simple moving average is calculated with a window width of +% CEIL((NUMPER+1)/2). Then a second simple moving average is calculated +% on the first moving average with the same window size. +% +% Weighted Moving Average: +% A weighted moving average is calculated with a weight vector, WEIGHTS. +% The length of the weight vector determines the size of the window. If +% larger weight factors are used for more recent prices and smaller factors +% for previous prices, the trend will be more responsive to recent changes. +% +% Modified Moving Average: +% The modified moving average calculation is similar to the simple moving +% average calculation. NUMPER can be thought of as the LAG of the simple +% moving average. The first modified moving average value, VO(NUMPER), is +% calculated the same way the first simple moving average value is +% calculated. All subsequent values are calculated by adding the new +% price and then subtracting the last average from the resulting sum. +% +% See also MEAN, PERAVG. + +% Reference: Achelis, Steven B., Technical Analysis From A To Z, +% Second Printing, McGraw-Hill, 1995, pg. 184-192 + +% Author: P. Wang +% Copyright 1999-2007 The MathWorks, Inc. +% $Revision: 1.1.4.2.2.1 $ $Date: 2007/03/28 18:03:40 $ + +% Transpose input if necessary. If nargin == 3, assume +% the user supplied a row oriented matrix. +if nargin == 4 + dim = varargin{2}; + if dim == 1 + % If input is column-oriented, transpose. + vin = vin'; + + elseif dim == 2 + % Default, do nothing. + + else + error('ftseries:tsmovavg:InvalidDIM', ... + 'Invalid DIM indicator. DIM = 1 (column vector); Dim = 2 (row vectors).'); + end + +elseif (nargin < 3) || (nargin > 4) + error('ftseries:tsmovavg:InvalidNumOfInputs', ... + 'Invalid number of inputs.'); +end + +% vinVars = # variables; observ = # obervations +[vinVars, observ] = size(vin); + +% Mov Avg modes and calculations +switch lower(mode(1)) + case 's' % Simple + % Lag window size + lag = varargin{1}; + + % Calculate simple mov avg + vout = simplema(vin, lag, vinVars, observ); + + case 'e' % Exponential + % Get the time periods + if nargin < 3 + timePer = 5; + + elseif nargin <= 4 + timePer = varargin{1}; + + else + error('ftseries:tsmovavg:InvalidNumOfInputsExponential', ... + 'Invalid number of inputs.'); + end + + if numel(timePer) ~= 1 || ((timePer <= 0) || (timePer >= observ)) + error('ftseries:tsmovavg:InvalidLag', ... + 'Lag must be scalar greater than 0 or less than the number of observations.'); + end + + % Calculate the exponential percentage + k = 2 / (timePer + 1); + + % + % EMA = (k * (CurrentPrice-PreviousPeriodEMA)) / PreviousPeriodEMA + % EMA = (k * (vin-vout)) / vout + % = K*vin + ((1-k) * vout) + % + + % Calculate the simple moving average for the first 'exp mov avg' + % value. + vout = nan(vinVars, observ); + vout(:, timePer) = sum(vin(:, 1:timePer), 2)/timePer; + + % K*vin; 1-k + kvin = vin(:, timePer:observ) * k; + oneK = 1-k; + + % First period calculation + vout(:, timePer) = kvin(:, 1) + (vout(:, timePer) * oneK); + + % Remaining periods calculation + for idx = timePer+1:observ + vout(:, idx) = kvin(:, idx-timePer+1) + (vout(:, idx-1) * oneK); + end + + case 't' % Triangular + % Get the time periods + if nargin < 3 + numPer = 5; + + elseif nargin <= 4 + numPer = varargin{1}; + + else + error('ftseries:tsmovavg:InvalidNumOfInputsTriangular', ... + 'Invalid number of inputs.'); + end + + % Window size + maPer = ceil((numPer + 1) / 2); + + % First moving average + movAvg1 = simplema(vin, maPer, vinVars, observ); + + % Second moving average + vout = simplema(movAvg1, maPer, vinVars, observ); + + case 'w' % Weighted + % Get the weights + if (nargin >= 3) && (nargin <= 4) + weights = varargin{1}; + + % Determine the length of the weights (window size) + [wi, lenWghts] = size(weights); + if wi ~= 1 + error('ftseries:tsmovavg:WeightMustBeRowVector', ... + 'The weights must a row oriented vector.'); + end + + else + error('ftseries:tsmovavg:InvalidNumOfInputsWeighted', ... + 'Invalid number of inputs.'); + end + + % Sum of weights + sumWghts = sum(weights); + + % Preallocate a vector of nans + vout = ones(size(vin))*nan; + + % repmat the weights vector to have the same number of rows as the input + weights = repmat(weights, size(vin, 1), 1); + + % Calculate the weighted mov avg + for idx = lenWghts:observ + vout(:, idx) = sum(weights.*vin(:, idx-lenWghts+1:idx), 2) / sumWghts; + end + + case 'm' % Modified + % Get the number of periods + if nargin < 3 + numPer = 5; + + elseif nargin <= 4 + numPer = varargin{1}; + + else + error('ftseries:tsmovavg:InvalidNumOfInputsExponential', ... + 'Invalid number of inputs.'); + end + + % Calculate simple mov avg. The first point of the modified moving + % average is calculated the same way the first point of the simple + % moving average is calculated. However, all subsequent points are + % calculated using the modified mov avg formula. + vout = simplema(vin, numPer, vinVars, observ); + + % Calculate modified mov avg + for idx = numPer+1:observ + % Remaining periods calculation + vout(:, idx) = vout(:, idx-1) + (vin(:, idx)-vout(:, idx-1))/numPer; + end + + otherwise + error('ftseries:tsmovavg:InvalidMode', ... + 'Invalid mode for moving average calculation.'); +end + +% Transpose output (if appropriate) so it is the same +% size as the input. +if exist('dim', 'var') && dim == 1 + vout = vout'; +end + +% ------------------------------------------------------------------------- +function simOut = simplema(vin,lag,vinVars,observ) +% Simple moving average + +if numel(lag) ~= 1 || ((lag <= 0) || (lag >= observ)) + error('ftseries:simplema:InvalidLag', ... + 'Lag must be scalar greater than 0 or less than the number of observations.'); +end + +% Preallocate a vector of nans +simOut = nan(size(vin)); + +for idx = 1:vinVars + % Simple moving average + ma = filter(ones(1,lag)/lag, 1, vin(idx,:)); + + % Fill in the NaN's vector + simOut(idx, lag:end) = ma(lag:end); +end + + +% [EOF] \ No newline at end of file diff --git a/wavediffusionstep.m b/wavediffusionstep.m new file mode 100644 index 0000000..dc0380e --- /dev/null +++ b/wavediffusionstep.m @@ -0,0 +1,137 @@ +function [E]=wavediffusionstep(E,i,c,cg,d,ho,sigma,dx,angle,Cbr,Cbed,alpha,angle0,N,M,A,AW,Awup,p,ii,jj,periodic,Holateral) + + +%if there is at least an active cell..., then calculate the diffusion +%if sum(Active==1)>0 +if sum(ho>0)>0 + + +%Use the fomrula of Mase +%use chain rule +%get this: + 1/2*CCg*Eyy (CCg)y*Ey +%solved the first one as diffusion +%solve second one as edvection + +%when you stretch the coordinates, you also stretch the dx!!! +stretchdx=1/cos(angle/180*pi); + + DD=(cg.*c)'; + DD(d<0.1 | A(i,:)==0)=0; + DDp=([DD(2:end); DD(end)]-[DD(1); DD(1:end-1)])/2;DDp(1)=0;DDp(end)=0; + s=[];S=0*p; + %%%%%%%%%%%%%%%%%% + for k=[-1 1]; + +% if periodic==0; q=p+k;if k==1;a=[1:M-1];elseif k==-1;a=[2:M];end +% elseif periodic==1; q=1+mod(p+k-1,M);a=[1:M]; +% end + q=1+mod(p+k-1,M);a=[1:M];%%KINDA force the periodic options in chosiing the cell (Dleaware wave diffusion) + + %Lateral diffraction (diffusion) + sigma=min(2*pi/1,max(2*pi/20,sigma)); + cg=min(100,max(1,cg)); + if length(sigma)>1 + value=alpha/dx^2 /2./sigma(p(a))' *0.5.*(DD(p(a))+DD(q(a)))/2 *dx*stretchdx./cg(p(a))'*(cos(angle/180*pi))^2; %*stretchdx + else + value=alpha/dx^2 /2./sigma *0.5.*(DD(p(a))+DD(q(a)))/2 *dx*stretchdx./cg(p(a))'*(cos(angle/180*pi))^2; %*stretchdx + end + + up=[];F=0; + if (k==1);UR=DDp;up=find(UR>0);F=UR(up);end + if (k==-1);UR=DDp;up=find(UR<0);F=-UR(up);end + if length(sigma)>1 + value(up)=value(up)+alpha/dx^2 /2./sigma(up)' .*F *dx*stretchdx./cg(up)'*(cos(angle/180*pi))^2; + else + value(up)=value(up)+alpha/dx^2 /2./sigma .*F *dx*stretchdx./cg(up)'*(cos(angle/180*pi))^2; + end + %note that you gor a diveded by dx^2. One dx is ffrom DDp (dx not included there) + %the other dx is from the derivate of E + + %NOTE: + %in both case you have *dx/cg, which is the time elapsed + %you need to strech the dx, becuase with the same cg it needs to + %travle farther if it is at anlgle + %then you get the cos^2 form the Mase formulae + + %REMOVE THE LATERAL DIFFUSION + if periodic==0;%DO IT ONLY IF IT IS NOT PERIODIC!!!!!! + value=value.*(AW(p(a))==0)';%zero diffusion at the lateral boundary if not periodic + end + + S(p(a))=S(p(a))+value; %exit from that cell + s=[s;-value]; %gain from the neigborh cell + end + + + + %Shoaling, bed friction, and breaking + s=[s;S+1]; + %s=[s;S+( (cg(i,p)+dx*stretchdx*(bedfr(i,p))+brk+whitecap)./cg(i+1,p))']; + + A = sparse(ii(abs(s)>0),jj(abs(s)>0),s(abs(s)>0));%solve the matrix inversion + %A = sparse(ii,jj,s);%solve the matrix inversion + + E=(A\E')'; + + + +% %Set the lateral bouundaries if not periodics. Two options: +% %no-gradient (+-1) or zero wave heigth (+-2) +% if periodic==0 %DO IT ONLY IF IT IS NOT PERIODIC!!!!!! +% if angle>0%angle<0 +% a=find(AW(:)==-1); %find the incomign boundary +% if a>0; +% b=find(AW(1+mod(-1+a-1,M))==0); +% E(a)=E(a(b)); +% end +% a=find(AW(:)==-2); %find the incomign boundary +% if a>0; +% b=find(AW(1+mod(-1+a-1,M))==0); %the first free cel; +% E(a)=E(a(b))*0+Holateral^2/16; +% end +% elseif angle<0%angle>0 +% a=find(AW(:)==1); %find the incomign boundary +% if a>0; +% b=find(AW(1+mod(-1+a+1,M))==0); +% E(a)=E(a(b)); +% end +% a=find(AW(:)==2); %find the incomign boundary +% if a>0; +% b=find(AW(1+mod(-1+a+1,M))==0); +% E(a)=E(a(b))*0+Holateral^2/16; +% end +% end +% end +% + + + + + + + + + + + + + + + + + + + %value=alpha/dx^2 /(2*(2*pi/T)) *((DD(p(a))+DD(q(a)))/2-1/2*DD(p(a))) *dx*stretchdx./cg(i+1,p(a))'; + %value=alpha/dx^2 /(2*(2*pi/T)) *0.5*(DD(p(a))+DD(q(a)))/2 *dx.*STRDX'./cg(i+1,p(a))'; %*stretchdx + %value=alpha/dx^2 /(2*(2*pi/T)).*min(DD(p(a)),DD(q(a))) *dx*stretchdx./cg(i+1,p(a))'; %GOOD + %value=alpha/dx^2 /(2*(2*pi/T)) *min(DD(p(a)),DD(q(a))) *dx*stretchdx; + %value=alpha/dx^2 /(2*(2*pi/T)) *min(DD(p(a)),DD(q(a))) *dx*stretchdx./cg(i+1,p(a))'; + %value=alpha/dx^2 ./max(kwave(i,p(a)),kwave(i,q(a)))' *dx*stretchdx; + %value=alpha/dx^2 ./kwave(i+1,p(a))' *dx*stretchdx; + + + + + +end + \ No newline at end of file diff --git a/wavek.m b/wavek.m new file mode 100644 index 0000000..3926f1f --- /dev/null +++ b/wavek.m @@ -0,0 +1,39 @@ +function K=wavek(F,H); +% function K=wavek(F,H); +% where K is wavenumber (rad/m) +% F is frequency (Hz) +% H is depth (m) +% +% Copyright (C) 2001, Lee Gordon, NortekUSA LLC +% KK=0*H; +% +% a=find(H>0.1); +% H=H(a); + +g=9.80171; + +% This routine use an approximate equation, then sharpens the result with +% one interpolation. The result is good to around 1 part in 10^-6. + +% The equation came out of a textbook, but I have long since forgotton +% which one. If you know, please tell me! lgordon@nortekusa.com + +e1=4*pi^2*F.^2.*H/g; %f4 = omega^2 * h1/g +e2 = 1+0.6666666*e1 + 0.355555555*e1.^2 + 0.1608465608*e1.^3 + ... + 0.0632098765*e1.^4 + 0.0217540484*e1.^5 + 0.0065407983*e1.^6; +e3 = +e1.^2 + e1./e2; +K1=sqrt(e3)./H; + +%compute error as basis for interpolation + +o1=sqrt(g*K1.*tanh(K1.*H)); +e1=o1.^2.*H/g; +e2 = 1+0.6666666*e1 + 0.355555555*e1.^2 + 0.1608465608*e1.^3 + ... + 0.0632098765*e1.^4 + 0.0217540484*e1.^5 + 0.0065407983*e1.^6; +e3 = +e1.^2 + e1./e2; +K2=sqrt(e3)./H; + +%interpolate +K=2*K1-K2; + +% KK(a)=K; \ No newline at end of file