Add files for new release

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(r<aw*Pedge(edg)*dt/dx);
+%edgG=edg;save edgG edgG
+
+%these cells will erode (the "high" cells)
+edger=edg(a);
+
+deltaY2=MASK*0;
+EdgeERY2=MASK*0;
+
+
+
+ %m  max heigth that is eroded, to avoid strange erosion of channels
+[N,M]=size(z);
+for i=1:length(edger)
+
+    %these are the adjacent cells
+    [I,J] = ind2sub(size(z),edger(i));
+    p=[sub2ind(size(z),mod(I+1-1,N)+1,J) sub2ind(size(z),mod(I-1-1,N)+1,J) sub2ind(size(z),I,mod(J+1-1,M)+1) sub2ind(size(z),I,mod(J-1-1,M)+1)];
+
+    %standard way
+    %a=find(MASK(p)==1 & A(p)==1); %only choses the adjacent cells if they are mudlfat and if they are standard cells
+
+    %to arode also at the opne boundary
+    a=find(MASK(p)==1); 
+
+
+ %pause(1)
+    if a>0 %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
+
+
+

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<hbedsheatresslim)=hbedsheatresslim;
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+%[Hs,Tp]=YeV_correction(Ff,wind,D);%[Hs,Tp]=YeV(Ff,wind,min(3,D));  %TRUCCO PER EVITARE LARGE WAVES IN CHANELS
+[Hs,Tp]=YeV(Ff,wind,D);%[Hs,Tp]=YeV(Ff,wind,min(3,D));  %TRUCCO PER EVITARE LARGE WAVES IN CHANELS
+HS(a)=Hs;
+TP(a)=Tp;TP(TP==0)=1;
+
+
+%do not diffuse in cells outside the MASK
+%HS=diffusefetch(MASK,HS,alpha,dx);
+%TP=diffusefetch(MASK,TP,alpha,dx);
+
+
+%hlimbedshearstress=0.5;
+%h=max(hlimbedshearstress,h);% to reduce the bed shear stress for very small water depth
+
+
+kwave=0*h;
+kk=wavek(1./TP(a),h(a));%kk=wavekFAST(1./Tp,D);
+kwave(a)=kk;
+kwave(kwave==0)=1;
+
+Um=pi*HS./(TP.*sinh(kwave.*h));
+
+cg=(2*pi./kwave./TP)*0.5.*(1+2*kwave.*h./(sinh(2*kwave.*h)));
+PW=cg*1030*9.8.*HS.^2/16;
+
+Um(MASK==0)=0;
+PW(MASK==0)=0;
+
+
+
+% h=[0.1:0.1:3.5];
+% [Hs,Tp]=YeV(3000,7,h);%[Hs,Tp]=YeV(Ff,wind,min(3,D));  %TRUCCO PER EVITARE LARGE WAVES IN CHANELS
+% kk=wavek(1./Tp,h);%kk=wavekFAST(1./Tp,D);
+% Um=pi*Hs./(Tp.*sinh(kk.*h));
+% ko=0.1/1000*3;
+% aw=Tp.*Um/(2*pi);
+% fw=0.00251*exp(5.21*(aw/ko).^-0.19);fw(aw/ko<pi/2)=0.3;
+% figure;plot(h,0.5*1000*0.015*Um.^2,h,0.5*1000*fw.*Um.^2)

YeV.m

@@ -0,0 +1,8 @@
+function [Hs,Tp]=YeV(fetch,wind,h)
+g=9.8;
+delta=h*g./wind.^2;
+chi=fetch*g./wind.^2;
+epsilon=3.64*10^-3*(tanh(0.493*delta.^0.75).*tanh(3.13*10^-3*chi.^0.57./tanh(0.493*delta.^0.75))).^1.74;
+ni=0.133*(tanh(0.331*delta.^1.01).*tanh(5.215*10^-4*chi.^0.73./tanh(0.331*delta.^1.01))).^-0.37;
+Hs=4*sqrt(wind.^4.*epsilon/g^2);
+Tp=wind./ni/g;

bedcreepponds.m

@@ -0,0 +1,70 @@
+function z=bedcreepponds(z,A,Active,Yreduction,crMUD,crMARSH,dx,dt,VEG,S,Qs,rbulk2,alphaMUD);
+
+A(Active==0)=0;
+
+%%%%%%%!!!!!!!$$$$$$
+%%%%%%%!!!!!!!$$$$$$
+%A(A==2)=1;%%trucco per fare creep also at the boundary
+
+
+%%downslope dependnets on sediment transport
+Qs=Qs/rbulk2;
+
+creep=A*0;
+creep(VEG==0)=crMUD+(alphaMUD*3600*24*Qs(VEG==0));          
+creep(VEG==1)=crMARSH;    
+
+% %%%%%%%!!!!!!!$$$$$$
+% %trucco per tenere il seaward basso
+% creep(end-100:end,:)=creep(end-100:end,:)+100;%
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+D=(creep)/(dx^2)*dt;%.*Yreduction;
+
+%consider the pond cells as A==0
+%A(S==1)=0; %DO NOT CREEP INTO PONDS!!!! NOT USED ANYMOREEEE!!!
+
+G=0*z;
+p=find(A==1);%exclude the NOLAND CELLS(A==0)!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+NN=length(p);G(p)=[1:NN];rhs=z(p);[N,M]=size(G);i=[];j=[];s=[];
+
+Spond=S;%%%%ATTENZIONE
+S=0*G; %This S ia a different emanign that ponds. it is just to store values
+
+[row col]=ind2sub(size(A),p);
+for k = [N -1 1 -N] 
+[a,q]=excludeboundarycell(k,N,M,p);
+a=a(A(q(a))==1);%only inclued the cells in whcih you can creep!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+value=(D(p(a))+D(q(a)))/2; %The orginal!!!!
+%value=min(D(p(a)),D(q(a))); %used ro make the wall of the marsh more vertical
+value=value.*min(Yreduction(p(a)),Yreduction(q(a)));
+
+
+% gradF=abs(z(p(a))-z(q(a)))/dx;
+% facNL=gradF>=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)));
+
+

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
+