function [s,par]=transus(s,par,varargin);
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University !
% Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski,            !
% Jaap van Thiel de Vries, Robert McCall                                  !
%                                                                         !
% d.roelvink@unesco-ihe.org                                               !
% UNESCO-IHE Institute for Water Education                                !
% P.O. Box 3015                                                           !
% 2601 DA Delft                                                           !
% The Netherlands                                                         !
%                                                                         !
% This library is free software; you can redistribute it and/or           !
% modify it under the terms of the GNU Lesser General Public              !
% License as published by the Free Software Foundation; either            !
% version 2.1 of the License, or (at your option) any later version.      !
%                                                                         !
% This library is distributed in the hope that it will be useful,         !
% but WITHOUT ANY WARRANTY; without even the implied warranty of          !
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU        !
% Lesser General Public License for more details.                         !
%                                                                         !
% You should have received a copy of the GNU Lesser General Public        !
% License along with this library; if not, write to the Free Software     !
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307     !
% USA                                                                     !
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

persistent aref cc ccb ccv ccvt chain costhm cu cub cumchain cuq3d cv cvb cvq3d dcdz delta delta_x deltas dsig dsigt dzbds eswbed eswmax exp_ero fac ftheta i isig j jg pbbedu pbbedv psi_x sbmtot sbtot sdif shields sigs sinthm ssmtot sub suq3d sus svb svq3d svs uau uav ueu_sed uev_sed um veu_sed vev_sed vm vmag2 ; 

 params;
 spaceparams;
 xmpi_module;
 interp;
 paramsconst;
% use vsmumod



 if isempty(i), i=0; end;
 if isempty(isig), isig=0; end;
 if isempty(j), j=0; end;
 if isempty(jg), jg=0; end;
 if isempty(exp_ero), exp_ero=0; end;

 if isempty(chain), chain=[]; end;
 if isempty(cumchain), cumchain=[]; end;
 if isempty(vmag2), vmag2=[]; end;
 if isempty(uau), uau=[]; end;
 if isempty(uav), uav=[]; end;
 if isempty(um), um=[]; end;
 if isempty(vm), vm=[]; end;
 if isempty(ueu_sed), ueu_sed=[]; end;
 if isempty(uev_sed), uev_sed=[]; end;
 if isempty(veu_sed), veu_sed=[]; end;
 if isempty(vev_sed), vev_sed=[]; end;
 if isempty(ccvt), ccvt=[]; end;
 if isempty(dcdz), dcdz=[]; end;
 if isempty(dsigt), dsigt=[]; end;
 if isempty(aref), aref=[]; end;
 if isempty(cc), cc=[]; end;
 if isempty(ccb), ccb=[]; end;
 if isempty(cu), cu=[]; end;
 if isempty(cv), cv=[]; end;
 if isempty(sus), sus=[]; end;
 if isempty(svs), svs=[]; end;
 if isempty(cub), cub=[]; end;
 if isempty(cvb), cvb=[]; end;
 if isempty(sub), sub=[]; end;
 if isempty(svb), svb=[]; end;
 if isempty(pbbedu), pbbedu=[]; end;
 if isempty(pbbedv), pbbedv=[]; end;
 if isempty(suq3d), suq3d=[]; end;
 if isempty(svq3d), svq3d=[]; end;
 if isempty(eswmax), eswmax=[]; end;
 if isempty(eswbed), eswbed=[]; end;
 if isempty(sigs), sigs=[]; end;
 if isempty(deltas), deltas=[]; end;
 if isempty(dsig), dsig=[]; end;
 if isempty(ccv), ccv=[]; end;
 if isempty(sdif), sdif=[]; end;
 if isempty(cuq3d), cuq3d=[]; end;
 if isempty(cvq3d), cvq3d=[]; end;
 if isempty(fac), fac=[]; end;

 if isempty(sinthm), sinthm=[]; end;
 if isempty(costhm), costhm=[]; end;

% Lodewijk: for direction of sediment transport (bed slope effect)
 if isempty(delta), delta=0; end;
 if isempty(delta_x), delta_x=0; end;
 if isempty(shields), shields=0; end;
 if isempty(ftheta), ftheta=0; end;
 if isempty(psi_x), psi_x=0; end;
 if isempty(sbtot), sbtot=0; end;
% Lodewijk: for magnitude of sediment transport (bed slope effect)
 if isempty(ssmtot), ssmtot=0; end;
 if isempty(dzbds), dzbds=0; end;
 if isempty(sbmtot), sbmtot=0; end;

%include 's.ind'
%include 's.inp'

if(~ ~isempty(vmag2));
vmag2=zeros(s.nx+1,s.ny+1);
uau=zeros(s.nx+1,s.ny+1);
uav=zeros(s.nx+1,s.ny+1);
ueu_sed=zeros(s.nx+1,s.ny+1);
uev_sed=zeros(s.nx+1,s.ny+1);
veu_sed=zeros(s.nx+1,s.ny+1);
vev_sed=zeros(s.nx+1,s.ny+1);
cu=zeros(s.nx+1,s.ny+1);
cv=zeros(s.nx+1,s.ny+1);
cc=zeros(s.nx+1,s.ny+1);
ccb=zeros(s.nx+1,s.ny+1);
fac=zeros(s.nx+1,s.ny+1,par.ngd);
sus=zeros(s.nx+1,s.ny+1);
svs=zeros(s.nx+1,s.ny+1);
cub=zeros(s.nx+1,s.ny+1);
cvb=zeros(s.nx+1,s.ny+1);
sub=zeros(s.nx+1,s.ny+1);
svb=zeros(s.nx+1,s.ny+1);
pbbedu=zeros(s.nx+1,s.ny+1);
pbbedv=zeros(s.nx+1,s.ny+1);
ccvt=zeros(s.nx+1,s.ny+1);
dcdz=zeros(s.nx+1,s.ny+1);
dsigt=zeros(s.nx+1,s.ny+1);
dsig=zeros(s.nx+1,s.ny+1,par.kmax);
ccv=zeros(s.nx+1,s.ny+1,par.kmax);
sdif=zeros(s.nx+1,s.ny+1,par.kmax);
um=zeros(s.nx+1,s.ny+1);
vm=zeros(s.nx+1,s.ny+1);
deltas=zeros(s.nx+1,s.ny+1);
sigs=zeros(s.nx+1,s.ny+1);
eswmax=zeros(s.nx+1,s.ny+1);
eswbed=zeros(s.nx+1,s.ny+1);
suq3d=zeros(s.nx+1,s.ny+1);
svq3d=zeros(s.nx+1,s.ny+1);
cuq3d=zeros(s.nx+1,s.ny+1,par.kmax);
cvq3d=zeros(s.nx+1,s.ny+1,par.kmax);
aref=zeros(s.nx+1,s.ny+1);
chain=zeros(1,par.kmax);
cumchain=zeros(1,par.kmax);
sinthm=zeros(s.nx+1,s.ny+1);
costhm=zeros(s.nx+1,s.ny+1);
% Lodewijk
delta_x   = 0.0d0;
% Lodewijk
shields   = 0.0d0;
% Lodewijk
ftheta    = 0.0d0;
% Lodewijk
psi_x     = 0.0d0;
% Lodewijk
sbtot     = 0.0d0;
% Lodewijk
delta     = 0.0d0;
uau(:)       = 0.0d0;
uav(:)       = 0.0d0;
um(:)        = 0.0d0;
vm(:)        = 0.0d0;
chain(:)     = 0.0d0;
cumchain(:)  = 0.0d0;
fac(:)       = 1.0d0;
exp_ero   = 0.0d0;
% generate sigma grid shape...
for isig=2:par.kmax;
chain(isig) = par.sigfac.^(isig-1);
cumchain(isig) = cumchain(isig-1)+chain(isig);
end; isig=fix(par.kmax+1);
end;

% use eulerian velocities
vmag2(:)     = s.ue.^2+s.ve.^2;
cu(:)        = 0.0d0;
cv(:)        = 0.0d0;
cub(:)       = 0.0d0;
cvb(:)       = 0.0d0;
s.dcsdx     = 0.0d0;
s.dcsdy     = 0.0d0;

sinthm(:) = sin(s.thetamean-s.alfaz);
costhm(:) = cos(s.thetamean-s.alfaz);

% short wave runup
if(par.swrunup==1 && par.struct==1);
[s,par]=hybrid(s,par);
end;

% compute turbulence due to wave breaking
if(par.lwt==1 || par.turb == turb_bore_averaged || par.turb == turb_wave_averaged);
[s,par]=waveturb(s,par);
end;

if(par.swave==1);
% include wave skewness and assymetry in sediment advection velocity
if(par.waveform==waveform_ruessink_vanrijn);
[s,par]=rvr(s,par);
elseif(par.waveform==waveform_vanthiel) ;
[s,par]=vt(s,par);
end;

end;

% calculate equilibrium concentration/sediment source
% Soulsby van Rijn
if((par.form==form_soulsby_vanrijn) ||(par.form==form_vanthiel_vanrijn)||(par.form==form_vanrijn1993));
% Soulsby van Rijn and Van Thiel de Vries & Reniers 2008 formulations
[s,par]=sedtransform(s,par);
end;

% compute reduction factor for sediment sources due to presence of hard layers
for jg = 1:par.ngd;
for j=1:s.ny+1;
for i=1:s.nx+1;
exp_ero = par.morfac.*par.dt./(1.0d0-par.por).*s.hh(i,j).*(s.ceqsg(i,j,jg).*s.pbbed(i,j,1,jg)./s.tsg(i,j,jg)+ s.ceqbg(i,j,jg).*s.pbbed(i,j,1,jg)./par.dt);
% limit erosion to available sediment on top of hard layer wwvv changed tiny into epsilon
fac(i,j,jg) = min(1.0d0,s.structdepth(i,j).*s.pbbed(i,j,1,jg)./max(eps,exp_ero) );
%if (fac(i,j,jg)*exp_ero > dzbed(i,j,1)*pbbed(i,j,1,jg)) then
%   limit erosion to available sand in top layer
%   fac(i,j,jg) = min(fac(i,j,jg),dzbed(i,j,1)*pbbed(i,j,1,jg)/max(tiny(0.0d0),exp_ero) )
%   write(*,*)'WARNING: expected erosion from top layer is larger than available sand in top layer'
%endif
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1+1);
end; jg =fix(par.ngd+1);

% compute diffusion coefficient
s.dc = par.facdc.*(par.nuh+par.nuhfac.*s.hh.*(s.dr./par.rho).^(1.0d0./3.0d0));

for jg = 1:par.ngd;
cc = s.ccg(:,:,jg);
if(s.d50(jg)>0.002d0);
% RJ: set ceqsg to zero for gravel.
% Dano: try without this fix cc = 0.0d0 ! Can be used to test total transport mode
end;
%
% X-direction
%
% Get ua in u points and split out in u and v direction
uau([1:s.nx],:) = 0.5.*(s.ua(1:s.nx,:).*costhm([1:s.nx],:)+s.ua(2:s.nx+1,:).*costhm([1:s.nx],:));
uav([1:s.nx],:) = 0.5.*(s.ua(1:s.nx,:).*sinthm([1:s.nx],:)+s.ua(2:s.nx+1,:).*sinthm([1:s.nx],:));
if(par.nz>1);
ueu_sed([1:s.nx],:) = 0.5.*(s.ue_sed(1:s.nx,:)+s.ue_sed(2:s.nx+1,:));
else;
ueu_sed(:)=s.ueu;
end;
veu_sed([1:s.nx],:) = 0.5.*(s.ve_sed(1:s.nx,:)+s.ve_sed(2:s.nx+1,:));
if(xmpi_isbot);
veu_sed(s.nx+1,:) = veu_sed(s.nx,:);
end;

% Compute vmagu including ua
%         s%vmagu = sqrt((s%uu+uau)**2+(s%vu+uav)**2)
s.vmagu = sqrt((ueu_sed+uau).^2+(veu_sed+uav).^2);
% sediment advection velocity for suspended load and bed load respectively
% REMARK: when vreps does not equal vv; no mass conservation
%         s%ureps = s%ueu+uau
%         s%urepb = s%ueu+uau  ! RJ maybe reduce this velocity?
s.ureps = ueu_sed+uau;
s.urepb = ueu_sed+uau;
%
for j=1:s.ny+1;
for i=1:s.nx;
if(s.ureps(i,j)>0.0d0);
% test cu(i,j)=cc(i,j)
cu(i,j)=par.thetanum.*cc(i,j)+(1.0d0-par.thetanum).*cc(min(i+1,s.nx),j);
cub(i,j)=par.thetanum.*s.pbbed(i,j,1,jg).*s.ceqbg(i,j,jg)+(1.0d0-par.thetanum).*s.pbbed(min(i+1,s.nx),j,1,jg).*s.ceqbg(min(i+1,s.nx),j,jg);
%cub(i,j)=par%thetanum*ccb(i,j)+(1.0d0-par%thetanum)*ccb(min(i+1,nx),j)
elseif(s.ureps(i,j)<0.0d0) ;
cu(i,j)=par.thetanum.*cc(i+1,j)+(1.0d0-par.thetanum).*cc(max(i,2),j);
cub(i,j)=par.thetanum.*s.pbbed(i+1,j,1,jg).*s.ceqbg(i+1,j,jg)+(1.0d0-par.thetanum).*s.pbbed(max(i,2),j,1,jg).*s.ceqbg(max(i,2),j,jg);
%cub(i,j)=par%thetanum*ccb(i+1,j)+(1.0d0-par%thetanum)*ccb(max(i,2),j)
else;
cu(i,j)=0.5d0.*(cc(i,j)+cc(i+1,j));
cub(i,j)=0.5d0.*(s.pbbed(i,j,1,jg).*s.ceqbg(i,j,jg)+s.pbbed(i+1,j,1,jg).*s.ceqbg(i+1,j,jg));
%cub(i,j)=0.5d0*(ccb(i,j)+ccb(i+1,j))
end;
s.dcsdx(i,j)=(cc(i+1,j)-cc(i,j))./s.dsu(i,j);
end; i=fix(s.nx+1);
end; j=fix(s.ny+1+1);
% wwvv dcdx(nx:1,:) is still untouched, correct this ofr the parallel case
%Robert
cu(s.nx+1,:) = cc(s.nx+1,:);
% wwvv fix this in parallel case
% wwvv in parallel version, there will be a discrepancy between the values
% of dzbdx(nx+1,:).
%wwvv so fix that
%
sus(:) = 0.0d0;
sub(:) = 0.0d0;
%
% suspended load, Lodewijk: no bed slope effect (yet)
sus=par.sus.*(cu.*s.ureps.*s.hu-s.dc.*s.hu.*s.dcsdx).*s.wetu;
% bed load, Lodewijk: no bed slope effect (yet)
sub=par.bed.*(cub.*s.urepb.*s.hu).*s.wetu;
%
% Originally bed slope effect of XBeach : par%bdslpeffmag = 1
% Original one, but only on bed load : par%bdslpeffmag = 2
if(par.bdslpeffmag == bdslpeffmag_roelv_total);
if(par.bermslope>0);
fmask=find((s.h./s.hu>1.0d0|(s.hu<1.0d0& par.instat==instat_stat_table)));NOTfmask=find(~(s.h./s.hu>1.0d0|(s.hu<1.0d0& par.instat==instat_stat_table)));
sus(fmask) = sus(fmask)-par.sus(fmask).*(10.0d0.*par.facsl.*cu(fmask).*s.vmagu.*s.hu.*(s.dzbdx-par.bermslope)).*s.wetu;
;
sus(NOTfmask) = sus(NOTfmask)-par.sus(NOTfmask).*(par.facsl.*cu(NOTfmask).*s.vmagu.*s.hu.*s.dzbdx).*s.wetu;
%end %where;
else;
sus = sus-par.sus.*(par.facsl.*cu.*s.vmagu.*s.hu.*s.dzbdx).*s.wetu;
end;
end;

if(par.bdslpeffmag == bdslpeffmag_roelv_total || par.bdslpeffmag == bdslpeffmag_roelv_bed);
if(par.bermslope>0);
fmask=find((s.h./s.hu>1.0d0 |(s.hu<1.0d0& par.instat==instat_stat_table)));NOTfmask=find(~(s.h./s.hu>1.0d0 |(s.hu<1.0d0& par.instat==instat_stat_table)));
sub(fmask) = sub(fmask)-par.bed.*(10.0d0.*par.facsl.*cub(fmask).*s.vmagu.*s.hu.*(s.dzbdx-par.bermslope)).*s.wetu;
;
sub(NOTfmask) = sub(NOTfmask)-par.bed.*(par.facsl.*cub(NOTfmask).*s.vmagu.*s.hu.*s.dzbdx).*s.wetu;
%end %where;
else;
sub = sub-par.bed.*(par.facsl.*cub.*s.vmagu.*s.hu.*s.dzbdx).*s.wetu;
end;
end;
%
%
% Y-direction
%
% Jaap: get ua in v points and split out in u and v direction
if(s.ny>0);
uau(:,[1:s.ny]) = 0.5.*(s.ua(:,1:s.ny).*costhm(:,[1:s.ny])+s.ua(:,2:s.ny+1).*costhm(:,[1:s.ny]));
uav(:,[1:s.ny]) = 0.5.*(s.ua(:,1:s.ny).*sinthm(:,[1:s.ny])+s.ua(:,2:s.ny+1).*sinthm(:,[1:s.ny]));
% Jaap
uau(:,s.ny+1) = uau(:,s.ny);
% Jaap
uav(:,s.ny+1) = uav(:,s.ny);
if(par.nz>1);
vev_sed(:,[1:s.ny]) = 0.5.*(s.ve_sed(:,1:s.ny)+s.ve_sed(:,2:s.ny+1));
else;
vev_sed(:)=s.vev;
end;
uev_sed(:,[1:s.ny]) = 0.5.*(s.ue_sed(:,1:s.ny)+s.ue_sed(:,2:s.ny+1));
if(xmpi_isright);
uev_sed(:,[1:s.ny+1]) = uev_sed(:,[1:s.ny]);
end;
else;
uau=s.ua.*costhm;
uav=s.ua.*sinthm;
uev_sed(:)=s.ue_sed;
vev_sed(:)=s.ve_sed;
end;
% Jaap: compute vmagv including ua
%         s%vmagv = sqrt((s%uv+uau)**2+(s%vv+uav)**2)
s.vmagv = sqrt((uev_sed+uau).^2+(vev_sed+uav).^2);
% sediment advection velocity for suspended load and bed load respectively
% REMARK: when vreps does not equal vv; no mass conservation
%         s%vreps = s%vev+uav
%         s%vrepb = s%vev+uav   ! RJ maybe reduce this velocity? Should be s%vv instead of s%vev?
s.vreps = vev_sed+uav;
% RJ maybe reduce this velocity? Should be s%vv instead of s%vev?
s.vrepb = vev_sed+uav;
%
if(s.ny>0);
for j=1:s.ny;
for i=1:s.nx+1;
if(s.vreps(i,j)>0);
cv(i,j)=par.thetanum.*cc(i,j)+(1.0d0-par.thetanum).*cc(i,min(j+1,s.ny));
cvb(i,j)=par.thetanum.*s.pbbed(i,j,1,jg).*s.ceqbg(i,j,jg)+(1.0d0-par.thetanum).*s.pbbed(i,min(j+1,s.ny),1,jg).*s.ceqbg(i,min(j+1,s.ny),jg);
%cvb(i,j)=par%thetanum*ccb(i,j)+(1.0d0-par%thetanum)*ccb(i,min(j+1,ny))
elseif(s.vreps(i,j)<0) ;
cv(i,j)=par.thetanum.*cc(i,j+1)+(1.0d0-par.thetanum).*cc(i,max(j,2));
cvb(i,j)=par.thetanum.*s.pbbed(i,j+1,1,jg).*s.ceqbg(i,j+1,jg)+(1.0d0-par.thetanum).*s.pbbed(i,max(j,2),1,jg).*s.ceqbg(i,max(j,2),jg);
%cvb(i,j)=par%thetanum*ccb(i,j+1)+(1.0d0-par%thetanum)*ccb(i,max(j,2))
else;
%Jaap: cc instead of cv
cv(i,j)=0.5d0.*(cc(i,j)+cc(i,j+1));
cvb(i,j)=0.5d0.*(s.pbbed(i,j,1,jg).*s.ceqbg(i,j,jg)+s.pbbed(i,j+1,1,jg).*s.ceqbg(i,j+1,jg));
%cvb(i,j)=0.5d0*(ccb(i,j)+ccb(i,j+1))
end;
%Jaap
s.dcsdy(i,j)=(cc(i,j+1)-cc(i,j))./s.dnv(i,j);

end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1);
% wwvv dcdy(:,ny+1) is not filled in, so in parallel case:
%Robert
cv(:,s.ny+1) = cc(:,s.ny+1);
% wwvv in parallel version, there will be a discrepancy between the values
% of dzbdy(:,ny+1).
% wwvv so fix that
else;
cv = cc;
cvb = s.ceqbg(:,:,jg);
% s%ny>0
end;
%
% Compute sedimnent transport in v-direction
%
svs(:) = 0.0d0;
svb(:) = 0.0d0;
% Suspended load
svs=par.sus.*(cv.*s.vreps.*s.hv-s.dc.*s.hv.*s.dcsdy).*s.wetv;
% Bed load
svb=par.bed.*(cvb.*s.vrepb.*s.hv).*s.wetv;
%
% Originally bed slope magnitude effect of XBeach : par%bdslpeffmag = 1
% Original one, but only on bed load : par%bdslpeffmag = 2
if(par.bdslpeffmag == bdslpeffmag_roelv_total);
svs = svs-par.sus.*(par.facsl.*cv.*s.vmagv.*s.hv.*s.dzbdy).*s.wetv;
end;
if(par.bdslpeffmag == bdslpeffmag_roelv_total || par.bdslpeffmag == bdslpeffmag_roelv_bed  );
svb = svb-par.bed.*(par.facsl.*cvb.*s.vmagv.*s.hv.*s.dzbdy).*s.wetv;
end;
%
%
% Bed slope magnitude effect (as Souslby intended) and change direction transport (see Van Rijn 1993 (section 7.2.6))
%
%
if(par.bdslpeffmag == bdslpeffmag_souls_total || par.bdslpeffmag == bdslpeffmag_souls_bed);
for j=1:s.ny+1;
for i=1:s.nx+1;
if((abs(sub(i,j)) > 0.000001d0) ||(abs(svb(i,j)) > 0.000001d0));
sbmtot = sqrt(  sub(i,j).^2.0d0  +  svb(i,j).^2.0d0   );
dzbds = s.dzbdx(i,j).*sub(i,j)./sbmtot + s.dzbdy(i,j).*svb(i,j)./sbmtot;
% dzbdn = s%dzbdx*Svb/Sbtot + s%dzbdy*Sub/Sbtot
sub(i,j) = sub(i,j).*(1.0d0 - par.facsl.*dzbds);
svb(i,j) = svb(i,j).*(1.0d0 - par.facsl.*dzbds);
%
end;
if(((abs(sus(i,j)) > 0.000001d0) ||(abs(svs(i,j)) > 0.000001d0))&& par.bdslpeffmag == bdslpeffmag_souls_total);
ssmtot = sqrt(  sus(i,j).^2.0d0  +  svs(i,j).^2.0d0  );
dzbds = s.dzbdx(i,j).*sus(i,j)./ssmtot + s.dzbdy(i,j).*svs(i,j)./ssmtot;
% dzbdn = s%dzbdx*Svb/Sbtot + s%dzbdy*Sub/Sbtot
sus(i,j) = sus(i,j).*(1.0d0 - par.facsl.*dzbds);
svs(i,j) = svs(i,j).*(1.0d0 - par.facsl.*dzbds);
end;
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1+1);
end;
%
% Lodewijk: modify the direction of the bed load transport based on the bed slope, see Van Rijn 1993 (section 7.2.6)
if(par.bdslpeffdir == bdslpeffdir_talmon);
for j=1:s.ny+1;
for i=1:s.nx+1;
if(((abs(s.urepb(i,j)) > 0.0001d0) ||(abs(s.vrepb(i,j)) > 0.0001d0))&&((abs(s.taubx(i,j)) > 0.0001d0) ||(abs(s.tauby(i,j)) > 0.0001d0)));
if(s.urepb(i,j) < 0.0d0);
% Angle between fluid velocity vector and the s%x-axis
delta_x = atan(s.vrepb(i,j)./s.urepb(i,j))+par.px;
else;
% Angle between fluid velocity vector and the s%x-axis
delta_x = atan(s.vrepb(i,j)./s.urepb(i,j));
end;
delta =(par.rhos-par.rho)./par.rho;
shields = sqrt(s.taubx(i,j).^2 + s.tauby(i,j).^2)./(delta.*par.rho.*par.g.*s.d50(jg));
% shields = (urepb(i,j)**2.0d0+vrepb(i,j)**2.0d0)*s%cf(i,j)/(par%g*D50(jg)*delta)
% Talmon
ftheta = 1.0d0./(9.0d0.*(s.d50(jg)./s.hh(i,j)).^0.3d0.*sqrt(shields));
psi_x = atan((sin(delta_x)-ftheta.*s.dzbdy(i,j))  ./(cos(delta_x)-ftheta.*s.dzbdx(i,j))  );
psi_x = par.bdslpeffdirfac.*(psi_x - delta_x)+delta_x;
% Magnitude of sediment transport without direction modifcation
sbtot = sqrt(  sub(i,j).^2.0d0  +  svb(i,j).^2.0d0  );
% Decompose the sediment transport again, know with the knowledge of the direction of the sediment transport vector
sub(i,j) = sbtot .* cos(psi_x);
svb(i,j) = sbtot .* sin(psi_x);
else;
sub(i,j) = 0.0d0;
svb(i,j) = 0.0d0;
end;
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1+1);
end;
%
%
%
for j=1:s.ny+1;
for i=1:s.nx;
if(sub(i,j)>0.0d0);
pbbedu(i,j) = s.pbbed(i,j,1,jg);
elseif(sub(i,j)<0.0d0) ;
pbbedu(i,j)= s.pbbed(i+1,j,1,jg);
else;
pbbedu(i,j)=0.5d0.*(s.pbbed(i,j,1,jg)+s.pbbed(i+1,j,1,jg));
end;
end; i=fix(s.nx+1);
end; j=fix(s.ny+1+1);
%
sub = pbbedu.*sub;
%
for j=1:s.ny;
for i=1:s.nx+1;
if(svb(i,j)>0);
pbbedv(i,j)=s.pbbed(i,j,1,jg);
elseif(svb(i,j)<0) ;
pbbedv(i,j)=s.pbbed(i,j+1,1,jg);
else;
pbbedv(i,j)=0.5d0.*(s.pbbed(i,j,1,jg)+s.pbbed(i,j+1,1,jg));
end;
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1);
%
svb = pbbedv.*svb;
%
% BRJ: implicit concentration update (compute sources first, sink must be computed after updating actual sed.conc.)
%
if(s.ny>0);
for j=2:s.ny;
for i=2:s.nx;
% Changed to hh from hold by RJ (13072009) !**2/max(hh(i,j),par%hmin)
s.ero(i,j,jg) = fac(i,j,jg).*s.hh(i,j).*s.ceqsg(i,j,jg).*s.pbbed(i,j,1,jg)./s.tsg(i,j,jg);
% depo_ex(i,j,jg) = max(hold(i,j),0.01d0)*cc(i,j)/Tsg(i,j,jg)
% BRJ: the volume in the water column is updated and not the volume concentration.
cc(i,j) =(par.dt.*s.tsg(i,j,jg))./(par.dt+s.tsg(i,j,jg)).*(s.hold(i,j).*cc(i,j)./par.dt -((sus(i,j).*s.dnu(i,j)-sus(i-1,j).*s.dnu(i-1,j)+svs(i,j).*s.dsv(i,j)-svs(i,j-1).*s.dsv(i,j-1)).*s.dsdnzi(i,j)-s.ero(i,j,jg)));

% Jaap: negative cc's are possible...
cc(i,j)=max(cc(i,j),0.0d0);
cc(i,j)=min(cc(i,j),par.cmax.*s.hh(i,j));
s.depo_ex(i,j,jg) = cc(i,j)./s.tsg(i,j,jg);
end; i=fix(s.nx+1);
end; j=fix(s.ny+1);

else;
j=1;
for i=2:s.nx;
% Changed to hh from hold by RJ (13072009) !**2/max(hh(i,j),par%hmin)
s.ero(i,j,jg) = fac(i,j,jg).*s.hh(i,j).*s.ceqsg(i,j,jg).*s.pbbed(i,j,1,jg)./s.tsg(i,j,jg);
% depo_ex(i,j,jg) = max(hold(i,j),0.01d0)*cc(i,j)/Tsg(i,j,jg)
% BRJ: the volume in the water column is updated and not the volume concentration.
cc(i,j) =(par.dt.*s.tsg(i,j,jg))./(par.dt+s.tsg(i,j,jg)).*(s.hold(i,j).*cc(i,j)./par.dt -((sus(i,j).*s.dnu(i,j)-sus(i-1,j).*s.dnu(i-1,j)).*s.dsdnzi(i,j)-s.ero(i,j,jg)));

% Jaap: negative cc's are possible...
cc(i,j)=max(cc(i,j),0.0d0);
cc(i,j)=min(cc(i,j),par.cmax.*s.hh(i,j));

s.depo_ex(i,j,jg) = cc(i,j)./s.tsg(i,j,jg);
end; i=fix(s.nx+1);
end;


cc = cc./s.hh;

% do lateral boundaries...
if(xmpi_istop);
cc(1,:)=cc(2,:);
s.ero(1,:,jg)=s.ero(2,:,jg);
s.depo_ex(1,:,jg)=s.depo_ex(2,:,jg);
end;
if(xmpi_isleft && s.ny>0);
cc(:,1)=cc(:,2);
s.ero(:,1,jg)=s.ero(:,2,jg);
s.depo_ex(:,1,jg)=s.depo_ex(:,2,jg);
end;
if(xmpi_istop);
cc(s.nx+1,:)=cc(s.nx+1-1,:);
s.ero(s.nx+1,:,jg)=s.ero(s.nx,:,jg);
s.depo_ex(s.nx+1,:,jg)=s.depo_ex(s.nx,:,jg);
end;
if(xmpi_isright && s.ny>0);
cc(:,s.ny+1)=cc(:,s.ny+1-1);
s.ero(:,s.ny+1,jg)=s.ero(:,s.ny,jg);
s.depo_ex(:,s.ny+1,jg)=s.depo_ex(:,s.ny,jg);
end;

% wwvv fix the first and last rows and columns of cc in parallel case
#ifdef usempi;
[cc]=xmpi_shift_ee(cc);
#end; jg =fix(par.ngd+1);
%
% wwvv border columns and rows of ccg Svg and Sug have to be communicated
s.ccg(:,:,jg) = cc;
s.svsg(:,:,jg) = svs;
s.susg(:,:,jg) = sus;
s.svbg(:,:,jg) = svb;
s.subg(:,:,jg) = sub;

% number of sediment fractions
end

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function [s,par]=bed_update(s,par,varargin);
persistent dafac dzavt dzb dzg dzleft dzmax dzt edg edg1 edg2 fac hav hinterland i id ie ii indsub indsus indsvb indsvs j j1 jd jdz je jg n_aval ndz savailable sdz sout ; 
mpi_sum=[];
 params;
 interp;
 spaceparams;
 xmpi_module;



 if isempty(i), i=0; end;
 if isempty(j), j=0; end;
 if isempty(j1), j1=0; end;
 if isempty(jg), jg=0; end;
 if isempty(ii), ii=0; end;
 if isempty(ie), ie=0; end;
 if isempty(id), id=0; end;
 if isempty(je), je=0; end;
 if isempty(jd), jd=0; end;
 if isempty(jdz), jdz=0; end;
 if isempty(ndz), ndz=0; end;
 if isempty(hinterland), hinterland=0; end;
 if isempty(indsus), indsus=[]; end;
 if isempty(indsub), indsub=[]; end;
 if isempty(indsvs), indsvs=[]; end;
 if isempty(indsvb), indsvb=[]; end;
 if isempty(dzb), dzb=0; end;
 if isempty(dzmax), dzmax=0; end;
 if isempty(dzt), dzt=0; end;
 if isempty(dzleft), dzleft=0; end;
 if isempty(sdz), sdz=0; end;
 if isempty(dzavt), dzavt=0; end;
 if isempty(fac), fac=0; end;
 if isempty(savailable), savailable=0; end;
 if isempty(dafac), dafac=0; end;
 if isempty(sout), sout=[]; end;
 if isempty(hav), hav=[]; end;
 if isempty(edg), edg=zeros(1,par.ngd); end;
 if isempty(edg1), edg1=zeros(1,par.ngd); end;
 if isempty(edg2), edg2=zeros(1,par.ngd); end;
 if isempty(dzg), dzg=zeros(1,par.ngd); end;
 dz=[];
 pb=[];
 if isempty(n_aval), n_aval=0; end;

%include 's.ind'
%include 's.inp'

if(~ ~isempty(sout));
sout=zeros(s.nx+1,s.ny+1);
hav=zeros(s.nx+1,s.ny+1);
indsus=(zeros(s.nx+1,s.ny+1,par.ngd));
indsub=(zeros(s.nx+1,s.ny+1,par.ngd));
indsvs=(zeros(s.nx+1,s.ny+1,par.ngd));
indsvb=(zeros(s.nx+1,s.ny+1,par.ngd));
end;

% Super fast 1D
if(s.ny==0);
j1 = 1;
else;
j1 = 2;
end;
s.dzbnow  = 0.0d0;
dzb    = 0.0d0;
if(par.t>=par.morstart && par.t < par.morstop && par.morfac > .999d0);
%
% bed_predict
%
% reduce sediment transports when hard layer comes to surface
% this step is mainly necessary at the transition from hard layers to sand
if(par.struct == 1);

for jg = 1:par.ngd;
indsus(:) = 0;
indsub(:) = 0;
indsvs(:) = 0;
indsvb(:) = 0;
sout(:)   = 0.0d0;
for j=j1:s.ny+1;
for i=2:s.nx+1;
% fluxes at i,j
% bed load s%u-direction
if(s.subg(i,j,jg) > 0.0d0);
indsub(i,j,jg) = 1;
sout(i,j) = sout(i,j) + s.subg(i,j,jg).*s.dnu(i,j);
end;
% fluxes at i-1,j
% bed load s%u-direction
if(s.subg(i-1,j,jg) < 0.0d0 );
sout(i,j) = sout(i,j) - s.subg(i-1,j,jg).*s.dnu(i-1,j);
end;
if(par.sourcesink==0);
% fluxes at i,j
% suspended load s%u-direction
if(s.susg(i,j,jg) > 0.0d0 );
indsus(i,j,jg) = 1;
sout(i,j) = sout(i,j) + s.susg(i,j,jg).*s.dnu(i,j);
end;
% fluxes at i-1,j
% suspended load s%u-direction
if(s.susg(i-1,j,jg) < 0.0d0 );
sout(i,j) = sout(i,j) - s.susg(i-1,j,jg).*s.dnu(i-1,j);
end;
end;
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1+1);
if(s.ny>0);
for j=j1:s.ny+1;
for i=2:s.nx+1;
% bed load s%v-direction
if(s.svbg(i,j,jg) > 0.0d0 );
indsvb(i,j,jg) = 1;
sout(i,j) = sout(i,j) + s.svbg(i,j,jg).*s.dsv(i,j);
end;
% fluxes at i,j-1
% bed load s%v-direction
if(s.svbg(i,j-1,jg) < 0.0d0 );
sout(i,j) = sout(i,j) - s.svbg(i,j-1,jg).*s.dsv(i,j-1);
end;
if(par.sourcesink==0);
% suspended load s%v-direction
if(s.svsg(i,j,jg) > 0.0d0 );
indsvs(i,j,jg) = 1;
sout(i,j) = sout(i,j) + s.svsg(i,j,jg).*s.dsv(i,j);
end;
% fluxes at i,j-1
% suspended load s%v-direction
if(s.svsg(i,j-1,jg) < 0.0d0 );
sout(i,j) = sout(i,j) - s.svsg(i,j-1,jg).*s.dsv(i,j-1);
end;
% sourcesink = 0
end;
%s%nx+1
end; i=fix(s.nx+1+1);
%s%ny+1
end; j=fix(s.ny+1+1);
%s%ny>0
end;
%
for j=j1:s.ny+1;
for i=2:s.nx+1;
savailable = s.structdepth(i,j).*s.pbbed(i,j,1,jg)./par.morfac./par.dt.*(1.0d0-par.por)./s.dsdnzi(i,j);
% reduction factor for cell outgoing sediment transports wwvv changed tiny into epsilon
fac  = min(1.0d0,savailable./max(sout(i,j),eps) );
% fix sediment transports for the presence of a hard layer; remind indSus etc are 1 in cases of cell outgoing transports
% updated S         oell outgoing transports                  cell incoming transports
if(fac < 1.0d0);
s.subg(i,j,jg)   = fac.*indsub(i,j,jg).*s.subg(i,j,jg)         +(1-indsub(i,j,jg)).*s.subg(i,j,jg);
s.subg(i-1,j,jg) = fac.*(1-indsub(i-1,j,jg)).*s.subg(i-1,j,jg) + indsub(i-1,j,jg).*s.subg(i-1,j,jg);
if(s.ny>0);
s.svbg(i,j,jg)   = fac.*indsvb(i,j,jg).*s.svbg(i,j,jg)         +(1-indsvb(i,j,jg)).*s.svbg(i,j,jg);
s.svbg(i,j-1,jg) = fac.*(1-indsvb(i,j-1,jg)).*s.svbg(i,j-1,jg) + indsvb(i,j-1,jg).*s.svbg(i,j-1,jg);
end;
if(par.sourcesink==0);
s.susg(i,j,jg)   = fac.*indsus(i,j,jg).*s.susg(i,j,jg)         +(1-indsus(i,j,jg)).*s.susg(i,j,jg);
s.susg(i-1,j,jg) = fac.*(1-indsus(i-1,j,jg)).*s.susg(i-1,j,jg) + indsus(i-1,j,jg).*s.susg(i-1,j,jg);
if(s.ny>0);
s.svsg(i,j,jg)   = fac.*indsvs(i,j,jg).*s.svsg(i,j,jg)         +(1-indsvs(i,j,jg)).*s.svsg(i,j,jg);
s.svsg(i,j-1,jg) = fac.*(1-indsvs(i,j-1,jg)).*s.svsg(i,j-1,jg) + indsvs(i,j-1,jg).*s.svsg(i,j-1,jg);
%s%ny = 0
end;
% sourcesink = 0
end;
%fac<1.0d0
end;
% s%nx+1
end; i=fix(s.nx+1+1);
%s%ny + 1
end; j=fix(s.ny+1+1);
%par%ngd
end; jg =fix(par.ngd+1);
%struct == 1
end;

if(s.ny>0);
for j=2:s.ny;
for i=2:s.nx;

% bed level changes per fraction in this morphological time step in meters sand including pores
% positive in case of erosion
if(par.sourcesink==0);
% dz from bed load transport gradients
% Dano, dz from sus transport gradients
dzg=par.morfac.*par.dt./(1.0d0-par.por).*(( s.susg(i,j,:).*s.dnu(i,j)-s.susg(i-1,j,:).*s.dnu(i-1,j) +s.svsg(i,j,:).*s.dsv(i,j)-s.svsg(i,j-1,:).*s.dsv(i,j-1) +s.subg(i,j,:).*s.dnu(i,j)-s.subg(i-1,j,:).*s.dnu(i-1,j)+s.svbg(i,j,:).*s.dsv(i,j)-s.svbg(i,j-1,:).*s.dsv(i,j-1) ).*s.dsdnzi(i,j)    );
elseif(par.sourcesink==1) ;
dzg=par.morfac.*par.dt./(1.0d0-par.por).*(s.ero(i,j,:)-s.depo_ex(i,j,:)   +( s.subg(i,j,:).*s.dnu(i,j)-s.subg(i-1,j,:).*s.dnu(i-1,j)+s.svbg(i,j,:).*s.dsv(i,j)-s.svbg(i,j-1,:).*s.dsv(i,j-1) ).*s.dsdnzi(i,j)    );
end;

% Simple bed update in case one fraction
if(par.ngd==1);

s.zb(i,j) = s.zb(i,j)-sum(sum(dzg));
% naamgeveing?
s.dzbnow(i,j) = s.dzbnow(i,j)-sum(sum(dzg));
s.dzbdt(i,j) = s.dzbnow(i,j)./par.dt;
s.sedero(i,j) = s.sedero(i,j)-sum(sum(dzg));
s.structdepth(i,j) = max(0.0d0,s.structdepth(i,j)-sum(sum(dzg)));

else;
% erosion/deposition rate of sand mass (m/s)
% positive in case of erosion
edg = dzg.*(1.0d0-par.por)./par.dt;


dz=s.dzbed(i,j,:);% =>
pb=s.pbbed(i,j,:,:);% =>

[par,s,i,j,dz,pb,edg]=update_fractions(par,s,i,j,dz,pb,edg,sum(sum(dzg)));

end;

% s%nx+1
end; i=fix(s.nx+1);
% s%ny+1
end; j=fix(s.ny+1);
else;
j=1;
for i=2:s.nx;
% bed level changes per fraction in this morphological time step in meters sand including pores
% positive in case of erosion
if(par.sourcesink==0);
% Dano, dz from sus transport gradients
dzg=par.morfac.*par.dt./(1.0d0-par.por).*(( s.susg(i,j,:).*s.dnu(i,j)-s.susg(i-1,j,:).*s.dnu(i-1,j) +s.subg(i,j,:).*s.dnu(i,j)-s.subg(i-1,j,:).*s.dnu(i-1,j) ).*s.dsdnzi(i,j) +(s.svsg(i,j,:)+s.svbg(i,j,:)).*par.lsgrad);
elseif(par.sourcesink==1) ;
dzg=par.morfac.*par.dt./(1.0d0-par.por).*(s.ero(i,j,:)-s.depo_ex(i,j,:)       +( s.subg(i,j,:).*s.dnu(i,j)-s.subg(i-1,j,:).*s.dnu(i-1,j) ).*s.dsdnzi(i,j) +(s.svsg(i,j,:)+s.svbg(i,j,:)).*par.lsgrad);
end;

% Simple bed update in case one fraction
if(par.ngd==1);

s.zb(i,j) = s.zb(i,j)-sum(sum(dzg));
s.dzbnow(i,j) = s.dzbnow(i,j)-sum(sum(dzg));
s.dzbdt(i,j) = s.dzbnow(i,j)./par.dt;
s.sedero(i,j) = s.sedero(i,j)-sum(sum(dzg));
s.structdepth(i,j) = max(0.0d0,s.structdepth(i,j)-sum(sum(dzg)));

% multiple fractions...
else;
% erosion/deposition rate of sand mass (m/s)
% positive in case of erosion
edg = dzg.*(1.0d0-par.por)./par.dt;

dz=s.dzbed(i,j,:);% =>
pb=s.pbbed(i,j,:,:);% =>

[par,s,i,j,dz,pb,edg]=update_fractions(par,s,i,j,dz,pb,edg,sum(sum(dzg)));

end;

% s%nx+1
end; i=fix(s.nx+1);
%s%ny>0
end;
#ifdef usempi;
[s.zb]=xmpi_shift_ee(s.zb);
#end;
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%
% Avalanching
%


if(par.avalanching==1);

for ii=1:round(par.morfac);
%aval=false
%Dano fix for MPI
n_aval=0;
s.dzbdx=0.0d0;
s.dzbdy=0.0d0;
for j=1:s.ny+1;
for i=1:s.nx;
s.dzbdx(i,j)=(s.zb(i+1,j)-s.zb(i,j))./s.dsu(i,j);
end; i=fix(s.nx+1);
end; j=fix(s.ny+1+1);
%
hav(:) = s.hh;
% Fix Hav for short wave runup:
if(par.swrunup == 1);
for j = 1:s.ny+1;
hinterland = 0;
for i = 1: s.nx+1;
if(hinterland == 0 && s.runup(j)+s.zs(round(s.iwl(j)),j) < s.zb(i,j) + par.eps);
hinterland = 1;
end;
if(s.wetz(i,j) == 1);
hav(i,j) = max(par.eps,(s.hh(i,j) + s.runup(j)));
elseif(hinterland == 0) ;
hav(i,j) =  max(par.eps,s.runup(j)+s.zs(round(s.iwl(j)),j)-s.zb(i,j) );
else;
hav(i,j) = par.eps;
end;
end; i = fix(s.nx+1+1);
end; j =fix(s.ny+1+1);
end;
%
%-1 Jaap -1 gives issues for bed updating at mpi boundaries
for i=2:s.nx;
for j=1:s.ny+1;
%if (max( max(hh(i,j),par%delta*H(i,j)), max(hh(i+1,j),par%delta*H(i+1,j)) )>par%hswitch+par%eps) then
% Jaap instead of s%hh
if(max(hav(i,j),hav(i+1,j))>par.hswitch+par.eps);
dzmax=par.wetslp;
% tricks: seaward of istruct (transition from sand to structure) wetslope is set to 0.03;
if(i>round(s.istruct(j)));
%dzmax = 0.03d0
dzmax = max(0.03d0,abs(s.dzbdx(i,j)).*0.99d0);
end;
else;
dzmax=par.dryslp;
end;

if(abs(s.dzbdx(i,j))>dzmax+1.0d-10 && s.structdepth(i+round(max(0.0d0,(abs(1.0d0).*sign(s.dzbdx(i,j)+eps)))),j)>par.eps);
n_aval=fix(n_aval+1);
dzb=(abs(1.0d0).*sign(s.dzbdx(i,j)+eps)).*(abs(s.dzbdx(i,j))-dzmax).*s.dsu(i,j);

if(dzb >= 0.0d0);
% index erosion point
ie = fix(i+1);
% index deposition point
id = fix(i);
% take into account varying grid sizes
dafac = s.dsdnzi(i,j)./s.dsdnzi(i+1,j);
% make sure dzb is not in conflict with par%dzmax
dzb=min(dzb,par.dzmax.*par.dt./s.dsu(i,j));
% make sure dzb is not larger than sediment layer thickness
dzb=min(dzb,s.structdepth(i+1,j));
else;
% index erosion point
ie = fix(i);
% index deposition point
id = fix(i+1);
% take into account varying grid sizes
dafac = s.dsdnzi(i+1,j)./s.dsdnzi(i,j);
dzb=max(dzb,-par.dzmax.*par.dt./s.dsu(i,j));
dzb=max(dzb,-s.structdepth(i,j));
end;


% Simple bed update in case one fraction
if(par.ngd == 1);

dzleft = abs(dzb);

s.zb(id,j) = s.zb(id,j)+dzleft.*dafac;
s.zb(ie,j) = s.zb(ie,j)-dzleft;
% naamgeveing?
s.dzbnow(id,j) = s.dzbnow(id,j)+dzleft.*dafac;
s.dzbnow(ie,j) = s.dzbnow(ie,j)-dzleft;
s.sedero(id,j) = s.sedero(id,j)+dzleft.*dafac;
s.sedero(ie,j) = s.sedero(ie,j)-dzleft;
s.structdepth(id,j) = max(0.0d0,s.structdepth(id,j)+dzleft.*dafac);
s.structdepth(ie,j) = max(0.0d0,s.structdepth(ie,j)-dzleft);

s.zs(id,j)  = s.zs(id,j)+dzleft.*dafac;
s.zs(ie,j)  = s.zs(ie,j)-dzleft;
s.dzav(id,j)= s.dzav(id,j)+dzleft.*dafac;
s.dzav(ie,j)= s.dzav(ie,j)-dzleft;

% multiple fractions...
else;

% now fix fractions....
dz = s.dzbed(ie,j,:);% =>
pb = s.pbbed(ie,j,:,:);% =>

% figure out how many depth layers (ndz) are eroded in point iii
sdz = 0;
ndz = 0;
while(sdz<abs(dzb));
ndz = fix(ndz+1);
sdz = sdz+dz(ndz);
end;

% now update bed and fractions by stepping through each layer seperately
dzleft = abs(dzb);
dzavt  = 0.0d0;

for jdz=1:ndz;

dzt = min(dz(jdz),dzleft);
dzleft = dzleft-dzt;

% erosion deposition per fraction (upwind or downwind); edg is positive in case of erosion
for jg=1:par.ngd;
% erosion    (dzt always > 0 )
edg2(jg) =  s.sedcal(jg).*dzt.*pb(jdz,jg).*(1.0d0-par.por)./par.dt;
% deposition (dzt always < 0 )
edg1(jg) = -s.sedcal(jg).*edg2(jg).*dafac;
end; jg=fix(par.ngd+1);

dzavt = dzavt + sum(sum(edg2)).*par.dt./(1.0d0-par.por);

% update bed in eroding point
[par,s,ie,j,s.dzbed(ie,j,:),s.pbbed(ie,j,:,:),edg2,dzavt]=update_fractions(par,s,ie,j,s.dzbed(ie,j,:),s.pbbed(ie,j,:,:),edg2,dzavt);

% update bed in deposition point
[par,s,id,j,s.dzbed(id,j,:),s.pbbed(id,j,:,:),edg1]=update_fractions(par,s,id,j,s.dzbed(id,j,:),s.pbbed(id,j,:,:),edg1,-dzavt.*dafac);

end; jdz=fix(ndz+1);

% update water levels and dzav
s.zs(ie,j)  = s.zs(ie,j)-dzavt;
s.dzav(ie,j)= s.dzav(ie,j)-dzavt;

s.zs(id,j)  = s.zs(id,j)+dzavt.*dafac;
s.dzav(id,j)= s.dzav(id,j)+dzavt.*dafac;

% yes/no multiple fractions
end;
end;
end; j=fix(s.ny+1+1);
end; i=fix(s.nx+1);
%JJ: update y slopes after avalanching in X-direction seems more appropriate
for j=1:s.ny;
for i=1:s.nx+1;
s.dzbdy(i,j)=(s.zb(i,j+1)-s.zb(i,j))./s.dnv(i,j);
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1);

%-1 Jaap -1 gives issues for bed updating at mpi boundaries
for j=2:s.ny;
for i=1:s.nx+1;
if(max(s.hh(i,j),s.hh(i,j+1))>par.hswitch+par.eps);
dzmax=par.wetslp;
else;
dzmax=par.dryslp;
end;
% Jaap
if(abs(s.dzbdy(i,j))>dzmax && s.structdepth(i,j+round(max(0.0d0,(abs(1.0d0).*sign(s.dzbdy(i,j)+eps)))))>par.eps);
n_aval=fix(n_aval+1);
dzb=(abs(1.0d0).*sign(s.dzbdy(i,j)+eps)).*(abs(s.dzbdy(i,j))-dzmax).*s.dnv(i,j);
%
if(dzb >= 0.0d0);
% index erosion point
je = fix(j+1);
% index deposition point
jd = fix(j);
% take into account varying grid sizes
dafac = s.dsdnzi(i,j)./s.dsdnzi(i,j+1);
dzb=min(dzb,par.dzmax.*par.dt./s.dnv(i,j));
dzb=min(dzb,s.structdepth(i,j+1));
else;
% index erosion point
je = fix(j);
% index deposition point
jd = fix(j+1);
% take into account varying grid sizes
dafac = s.dsdnzi(i,j+1)./s.dsdnzi(i,j);
dzb=max(dzb,-par.dzmax.*par.dt./s.dnv(i,j));
dzb=max(dzb,-s.structdepth(i,j));
end;

% Simple bed update in case one fraction
if(par.ngd == 1);

dzleft = abs(dzb);

s.zb(i,jd) = s.zb(i,jd)+dzleft.*dafac;
s.zb(i,je) = s.zb(i,je)-dzleft;
% naamgeveing?
s.dzbnow(i,jd) = s.dzbnow(i,jd)+dzleft.*dafac;
s.dzbnow(i,je) = s.dzbnow(i,je)-dzleft;
s.sedero(i,jd) = s.sedero(i,jd)+dzleft.*dafac;
s.sedero(i,je) = s.sedero(i,je)-dzleft;
s.structdepth(i,jd) = max(0.0d0,s.structdepth(i,jd)+dzleft.*dafac);
s.structdepth(i,je) = max(0.0d0,s.structdepth(i,je)-dzleft);

s.zs(i,jd)  = s.zs(i,jd)+dzleft.*dafac;
s.zs(i,je)  = s.zs(i,je)-dzleft;
s.dzav(i,jd)= s.dzav(i,jd)+dzleft.*dafac;
s.dzav(i,je)= s.dzav(i,je)-dzleft;

% multiple fractions...
else;

dz = s.dzbed(i,je,:);% =>
pb = s.pbbed(i,je,:,:);% =>

% figure out how many depth layers (ndz) are affected
sdz = 0;
ndz = 0;
while(sdz<abs(dzb));
ndz = fix(ndz+1);
sdz = sdz+dz(ndz);
end;

% now update bed and fractions by stepping through each layer seperately
dzleft = abs(dzb);
dzavt  = 0.0d0;

for jdz=1:ndz;
dzt = min(dz(jdz),dzleft);
dzleft = dzleft-dzt;

% erosion deposition per fraction (upwind or downwind); edg is positive in case of erosion
for jg=1:par.ngd;
% erosion    (dzt always > 0 )
edg2(jg) = s.sedcal(jg).*dzt.*pb(jdz,jg).*(1.0d0-par.por)./par.dt;
% deposition (dzt always < 0 )
edg1(jg) = -s.sedcal(jg).*edg2(jg).*dafac;
end; jg=fix(par.ngd+1);

dzavt = dzavt + sum(sum(edg2)).*par.dt./(1.0d0-par.por);

% upwind point
[par,s,i,je,s.dzbed(i,je,:),s.pbbed(i,je,:,:),edg2,dzavt]=update_fractions(par,s,i,je,s.dzbed(i,je,:),s.pbbed(i,je,:,:),edg2,dzavt);

% downwind point
[par,s,i,jd,s.dzbed(i,jd,:),s.pbbed(i,jd,:,:),edg1]=update_fractions(par,s,i,jd,s.dzbed(i,jd,:),s.pbbed(i,jd,:,:),edg1,-dzavt.*dafac);

end; jdz=fix(ndz+1);

% update water levels and dzav
s.zs(i,je)  = s.zs(i,je)-dzavt;
s.dzav(i,je)= s.dzav(i,je)-dzavt;

s.zs(i,jd)  = s.zs(i,jd)+dzavt.*dafac;
s.dzav(i,jd)= s.dzav(i,jd)+dzavt.*dafac;

%yes/no multiple fractions
end;
end;
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1);
% write(*,*)'t ',par%t,'ii ',ii,' n_aval ',n_aval,' rank ',xmpi_rank
% Dano: in parallel version bed update must take place AFTER EACH ITERATION
#ifdef usempi;
[n_aval  ,mpi_sum]=xmpi_allreduce(n_aval  ,mpi_sum);
#end; ii=fix(round(par.morfac)+1);
if(n_aval>0);
#ifdef usempi;
[s.zb]=xmpi_shift_ee(s.zb);
#end;
else;
 tempBreak=1;break;
end;

end

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function [par,s,i,j,dz,pb,edg,dzb]=update_fractions(par,s,i,j,dz,pb,edg,dzb,varargin);

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%     Copyright (C) 2009 Technische Universiteit Delft
%        Bram van Prooijen
%        b.c.vanprooijen@tudelft.nl
%      +31(0)15 2784070
%        Faculty of Civil Engineering and Geosciences
%        department of Hydraulic Engineering
%      PO Box 5048
%        2600 GA Delft
%        The Netherlands
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

persistent a ap b dzb_loc dzbt ed fac jd jg sm zbold ; 

 params;
 spaceparams;
 xmpi_module;




 if isempty(jg), jg=0; end;
 if isempty(jd), jd=0; end;
 if isempty(ed), ed=0; end;
 if isempty(zbold), zbold=0; end;
 if isempty(dzbt), dzbt=0; end;
 if isempty(fac), fac=0; end;
 if isempty(dzb_loc), dzb_loc=0; end;

 if isempty(ap), ap=[]; end;
 if isempty(b), b=[]; end;
 if isempty(sm), sm=[]; end;
 if isempty(a), a=[]; end;

if(~ ~isempty(ap));
ap=zeros(1,s.nd(1,1));
b=zeros(1,s.nd(1,1));
sm=zeros(s.nd(1,1),par.ngd);
a=zeros(s.nd(1,1),3);
end;
%TODO, dzb_loc is not initialized can be Nan, leading to infinite loop
dzb_loc = dzb;

%!!initialize Sm

ed = sum(sum(edg));

% do t_sub=1,nt_sub  !loop over subtimesteps
while(abs(dzb_loc) > 0.0d0);
% make sure erosion (dzg is positive) is limited to thickness of variable layer
dzbt     = min(dzb_loc,dz(par.nd_var));
% make sure deposition (dzg is negative) is limited to thickness of first layer below variable layer
dzbt     = max(dzbt,-par.frac_dz.*dz(par.nd_var+1));

% factor over mass change in cells to limit erosion and deposition
fac = dzbt./dzb_loc;

% update dzb
dzb_loc = dzb_loc-dzbt;

for jg=1:par.ngd;

a(:)=0.;
ap(:)=0.;
b(:)=0.;

sm(:,jg)=pb(:,jg).*dz.*(1.0d0-par.por);

%!!build matrix A
switch(par.nd_var);
case{1};
%in this case: A=0
case{2};
a(1,[1:3])=[0.0d0           ,  min(ed,0.0d0) ,  max(ed,0.0d0) ];
a(2,[1:3])=[-min(ed,0.0d0)  , -max(ed,0.0d0) ,  0.0d0         ];
case{3:1000};
a(1,[1:3])=[0.0d0           ,  min(ed,0.0d0) ,  max(ed,0.0d0) ];

a([2:par.nd_var-1],1)=-min(ed,0.0d0);
a([2:par.nd_var-1],2)=-abs(ed);
a([2:par.nd_var-1],3)=max(ed,0.0d0);

a(par.nd_var,[1:3])=[-min(ed,0.0d0) , -max(ed,0.0d0) ,   0.0d0    ];
end %select;

%!!determine RHS

% Ap makes sure that layer nd_var varies in thickness in time
% Ap = 0 with single fraction (in case(1))
ap(1) = sum(sum(a(1,[2:3]).*pb([1:2],jg)));
for jd = 2:par.nd_var;
ap(jd) = sum(sum(a(jd,:).*pb([jd-1:jd+1],jg)));
end; jd =fix(par.nd_var+1);
ap(par.nd_var+1) = sum(sum(a(par.nd_var+1,[1:2]).*pb([par.nd_var:par.nd_var+1],jg)));

% b represents the actual erosion and deposition in the top layer.
% However, the thickness of the top layer remains constant and instead layer nd_var will breath in thickness
b(1) = -edg(jg);

%!!update Sm

% Sm is the sediment mass per fraction per layer
% Sm(1:par%nd_var+1,jg) = Sm(1:par%nd_var+1,jg) + par%dt/nt_sub*(Ap(1:par%nd_var+1)+b(1:par%nd_var+1))
sm([1:par.nd_var+1],jg) = sm([1:par.nd_var+1],jg) + par.dt.*fac.*(ap([1:par.nd_var+1])+b([1:par.nd_var+1]));

%fractions
end; jg=fix(par.ngd+1);

% From Sm we can compute the new fraction ratios per layer and the layer thickness...
for jd=1:par.nd_var+1;
if(sum(sum(sm(jd,:)))>0);
pb(jd,:) = sm(jd,:)./sum(sum(sm(jd,:)));
else;
pb(jd,:) = pb(jd,:);
end; jd=fix(par.nd_var+1+1);
dz(jd) = sum(sum(sm(jd,:)))./(1.0d0-par.por);
end;

%!! modify grid

%merge two upper layers in case of erosion
if(dz(par.nd_var) < par.merge.*dz(par.nd_var+1));
forall(jg=1:par.ngd);
pb(par.nd_var,jg) =(dz(par.nd_var).*pb(par.nd_var,jg) + dz(par.nd_var+1).*pb(par.nd_var+1,jg))./(dz(par.nd_var)+dz(par.nd_var+1));
pb([par.nd_var+1:s.nd(i,j)-1],jg) = pb([par.nd_var+2:s.nd(i,j)],jg);
pb(s.nd(i,j),jg) = pb(s.nd(i,j),jg);
endforall;
s.z0bed(i,j) = s.z0bed(i,j)-dz(par.nd_var+1);
dz(par.nd_var) = dz(par.nd_var+1)+dz(par.nd_var);
end;
%split upper layer in case of sedimentation
if(dz(par.nd_var)>par.split.*dz(par.nd_var+1));
pb([par.nd_var+1:s.nd(i,j)],:) = pb([par.nd_var:s.nd(i,j)-1],:);
s.z0bed(i,j) = s.z0bed(i,j)+dz(par.nd_var+1);
dz(par.nd_var) = dz(par.nd_var)-dz(par.nd_var+1);
end;
% nt_sub
end;

pb = max(0.0d0,min(pb,1.0d0));

zbold = s.zb(i,j);
s.zb(i,j) = s.z0bed(i,j)+sum(sum(dz));
s.dzbnow(i,j) = s.dzbnow(i,j)+(s.zb(i,j)-zbold);
s.sedero(i,j) = s.sedero(i,j)+(s.zb(i,j)-zbold);
s.structdepth(i,j) = max(0.0d0,s.structdepth(i,j)+(s.zb(i,j)-zbold));

end %subroutine update_fractions

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function [s,par]=sedtransform(s,par,varargin);
persistent a alpha1 alpha2 alpha_ml asb ass b b2 beta ca cc1 cc2 cd ceqb ceqs dcf dcfin delta dhdx dhdy dster fallvelredfac fc fslope fw hfac hloc i j jg kl kscr ksw kvis m me ml onethird phi psi refa shieldscrit srfrhee srftotal sster taubw taucr tautotal te term1 ts twothird uandv ucr ucrb ucrc ucrs ucrw ue uorb urms2 used uw vero vmg w ws0 wster z0 ; 

 params;
 spaceparams;
 readkey_module;
 xmpi_module;



 if isempty(i), i=0; end;
 if isempty(j), j=0; end;
 if isempty(jg), jg=0; end;
 if isempty(z0), z0=0; end;
 if isempty(dcf), dcf=0; end;
 if isempty(dcfin), dcfin=0; end;
 if isempty(ml), ml=0; end;
 if isempty(te), te=0; end;
 if isempty(sster), sster=0; end;
 if isempty(cc1), cc1=0; end;
 if isempty(cc2), cc2=0; end;
 if isempty(wster), wster=0; end;
 if isempty(ass), ass=0; end;
 if isempty(kl), kl=0; end;
 if isempty(alpha_ml), alpha_ml=0; end;
 if isempty(alpha1), alpha1=0; end;
 if isempty(alpha2), alpha2=0; end;
 if isempty(beta), beta=0; end;
 if isempty(psi), psi=0; end;
 if isempty(delta), delta=0; end;
 if isempty(kvis), kvis=0; end;
 if isempty(onethird), onethird=0; end;
 if isempty(twothird), twothird=0; end;
 if isempty(phi), phi=0; end;
 if isempty(dster), dster=[]; end;
 if isempty(ws0), ws0=[]; end;
 if isempty(shieldscrit), shieldscrit=[]; end;
 if isempty(vmg), vmg=[]; end;
 if isempty(cd), cd=[]; end;
 if isempty(asb), asb=[]; end;
 if isempty(dhdx), dhdx=[]; end;
 if isempty(dhdy), dhdy=[]; end;
 if isempty(ts), ts=[]; end;
 if isempty(hfac), hfac=[]; end;
 if isempty(urms2), urms2=[]; end;
 if isempty(ucr), ucr=[]; end;
 if isempty(ucrc), ucrc=[]; end;
 if isempty(ucrw), ucrw=[]; end;
 if isempty(term1), term1=[]; end;
 if isempty(b2), b2=[]; end;
 if isempty(srftotal), srftotal=[]; end;
 if isempty(srfrhee), srfrhee=[]; end;
 if isempty(vero), vero=[]; end;
 if isempty(ucrb), ucrb=[]; end;
 if isempty(ucrs), ucrs=[]; end;
 if isempty(uandv), uandv=[]; end;
 if isempty(b), b=[]; end;
 if isempty(fslope), fslope=[]; end;
 if isempty(hloc), hloc=[]; end;
 if isempty(ceqs), ceqs=[]; end;
 if isempty(ceqb), ceqb=[]; end;
 if isempty(fallvelredfac), fallvelredfac=[]; end;
 if isempty(w), w=[]; end;
 if isempty(uorb), uorb=[]; end;
 if isempty(a), a=[]; end;
 if isempty(ksw), ksw=[]; end;
 if isempty(fw), fw=[]; end;
 if isempty(uw), uw=[]; end;
 if isempty(taubw), taubw=[]; end;
 if isempty(fc), fc=[]; end;
 if isempty(tautotal), tautotal=[]; end;
 if isempty(taucr), taucr=[]; end;
 if isempty(used), used=[]; end;
 if isempty(ue), ue=[]; end;
 if isempty(ca), ca=0; end;
 if isempty(refa), refa=0; end;
 if isempty(me), me=0; end;
 if isempty(m), m=0; end;
 if isempty(kscr), kscr=0; end;
%include 's.ind'
%include 's.inp'

if(~ ~isempty(vmg));
vmg=zeros(s.nx+1,s.ny+1);
term1=zeros(s.nx+1,s.ny+1);
b2=zeros(s.nx+1,s.ny+1);
cd=zeros(s.nx+1,s.ny+1);
asb=zeros(s.nx+1,s.ny+1);
ucr=zeros(s.nx+1,s.ny+1);
ucrc=zeros(s.nx+1,s.ny+1);
ucrw=zeros(s.nx+1,s.ny+1);
urms2=zeros(s.nx+1,s.ny+1);
hloc=zeros(s.nx+1,s.ny+1);
ts=zeros(s.nx+1,s.ny+1);
ceqs=zeros(s.nx+1,s.ny+1);
ceqb=zeros(s.nx+1,s.ny+1);
% Lodewijk
srftotal=zeros(s.nx+1,s.ny+1);
% Lodewijk
ucrb=zeros(s.nx+1,s.ny+1);
% Lodewijk
ucrs=zeros(s.nx+1,s.ny+1);
% Lodewijk
srfrhee=zeros(s.nx+1,s.ny+1);
% Lodewijk
vero=zeros(s.nx+1,s.ny+1);
% Lodewijk
fallvelredfac=zeros(s.nx+1,s.ny+1);
% Lodewijk
w=zeros(s.nx+1,s.ny+1,par.ngd);
dster=zeros(1,par.ngd);
ws0=zeros(1,par.ngd);
shieldscrit=zeros(1,par.ngd);
dhdx=zeros(s.nx+1,s.ny+1);
dhdy=zeros(s.nx+1,s.ny+1);
uandv=zeros(s.nx+1,s.ny+1);
b=zeros(s.nx+1,s.ny+1);
fslope=zeros(s.nx+1,s.ny+1);
hfac=zeros(s.nx+1,s.ny+1);

used=zeros(s.nx+1,s.ny+1);
ue=zeros(s.nx+1,s.ny+1);
uorb=zeros(s.nx+1,s.ny+1);
a=zeros(s.nx+1,s.ny+1);
ksw=zeros(s.nx+1,s.ny+1);
fw=zeros(s.nx+1,s.ny+1);
uw=zeros(s.nx+1,s.ny+1);
taubw=zeros(s.nx+1,s.ny+1);
fc=zeros(s.nx+1,s.ny+1);
tautotal=zeros(s.nx+1,s.ny+1);
taucr=zeros(s.nx+1,s.ny+1);

vmg(:) = 0.0d0;
onethird=1.0d0./3.0d0;
twothird=2.0d0./3.0d0;
% Angle of internal friction
phi = par.reposeangle./180.*par.px;
%do  only once, not necessary every time step % Robert:
for jg=1:par.ngd;
% cjaap: compute fall velocity with simple expression from Ahrens (2000)
te    = 20.0d0;
% Van rijn, 1993
kvis  = 4.0d0./(20.0d0+te).*1d-5;
sster = s.d50(jg)./(4.*kvis).*sqrt((par.rhos./par.rho-1).*par.g.*s.d50(jg));
cc1   = 1.06d0.*tanh(0.064d0.*sster.*exp(-7.5d0./sster.^2));
cc2    = 0.22d0.*tanh(2.34d0.*sster.^(-1.18d0).*exp(-0.0064d0.*sster.^2));
wster = cc1+cc2.*sster;
ws0(jg) = wster.*sqrt((par.rhos./par.rho-1.0d0).*par.g.*s.d50(jg));
% RJ: for modeling gravel
delta =(par.rhos-par.rho)./par.rho;
dster(jg)=(delta.*par.g./1.0d-12).^onethird.*s.d50(jg);
% Soulsby, 1997  (eq 76)
shieldscrit(jg) = 0.3d0./(1+1.2d0.*dster(jg))+0.055d0.*(1-exp(-0.02d0.*dster(jg)));
if(par.fallvelred==0);
w(:,:,jg) = ws0(jg);
end;
end; jg=fix(par.ngd+1);
end;


% Lodewijk: calculate fall velocity reduction coefficient based on concentration of previous time step
% Rowe (1987) made an estimation of the exponent alpha_ml by fitting a logarithmic function on a dataset of Richardson and Zaki (1954).
if(par.fallvelred==1);
for jg=1:par.ngd;
alpha_ml = 2.35d0.*(2.0d0+0.175d0.*(ws0(jg).*s.d50(jg)./kvis).^0.75d0)./(1.0d0+0.175d0.*(ws0(jg).*s.d50(jg)./kvis).^0.75d0);
fallvelredfac =(1.0d0-s.ccg(:,:,jg)).^(alpha_ml);
w(:,:,jg) = ws0(jg).*fallvelredfac;
end; jg=fix(par.ngd+1);
end;
%
% hloc   = max(hh,0.01d0) ! Jaap
hloc(:) = max(s.hh,0.01);
% Compute mean fall velocity
par.ws=0.0d0;
for jg=1:par.ngd;
par.ws=par.ws+w(1,1,jg);
end; jg=fix(par.ngd+1);
par.ws=par.ws./par.ngd;

%
% compute near bed turbulence
%
% due to short waves

if(par.swave==1);
% wave breaking induced turbulence due to short waves
for j=1:s.ny+1;
for i=1:s.nx+1;
% compute mixing length
% ML = 2*R(i,j)*par%Trep/(par%rho*c(i,j)*max(H(i,j),par%eps))
ml = sqrt(2.*s.r(i,j).*par.trep./(par.rho.*max(s.c(i,j),1d-10)));
% ML = 0.9d0*H(i,j)
ml = min(ml,hloc(i,j));
% exponential decay turbulence over depth
dcfin = exp(min(100.0d0,hloc(i,j)./max(ml,0.01d0)));
dcf = min(1.0d0,1.0d0./(dcfin-1.0d0));
% Jaap: new approach: compute kb based on waveturb result
s.kb(i,j) = s.kturb(i,j).*dcf;
%
if(par.turb == turb_bore_averaged);
s.kb(i,j) = s.kb(i,j).*par.trep./s.tbore(i,j);
end;

end; i=fix(s.nx+1+1);
% Jaap: rundown jet creating additional turbulence
if(par.swrunup==1);
s.kb(round(s.istruct(j)),j) = s.kb(round(s.istruct(j)),j) + par.jetfac.*(s.e(round(s.istruct(j)),j).*s.strucslope(j).*sqrt(par.g./s.hh(round(s.istruct(j)),j))).^twothird;
end;
end; j=fix(s.ny+1+1);
elseif(par.nonh==1) ;
for j=1:s.ny+1;
for i=1:s.nx+1;
%ML=max(s%rolthick(i,j),0.07d0*max(s%hh(i,j),par%hmin))
%s%kb(i,j) = s%kturb(i,j)*ML/max(s%hh(i,j),par%hmin) ! Simpler expression
% even simpler expression :-)
s.kb(i,j) = s.kturb(i,j);
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1+1);
%par%swave == 1
end;

% switch to include long wave stirring
if(par.lws==1);
vmg  = sqrt(s.ue.^2+s.ve.^2);
elseif(par.lws==0) ;
% vmg lags on actual mean flow; but long wave contribution to mean flow is included...
vmg =(1.0d0-1.0d0./par.cats./par.trep.*par.dt).*vmg +(1.0d0./par.cats./par.trep.*par.dt).*sqrt(s.ue.^2+s.ve.^2);
end;

urms2(:) = s.urms.^2.0d0+1.45d0.*s.kb;

for jg = 1:par.ngd;

ts       = par.tsfac.*hloc./w(:,:,jg);
s.tsg(:,:,jg) = max(ts,par.tsmin);
%
% calculate treshold velocity Ucr
%
% Soulsby van Rijn
if(par.form==form_soulsby_vanrijn);
if(s.d50(jg)<=0.0005d0);
ucr=0.19d0.*s.d50(jg).^0.1d0.*log10(4.0d0.*hloc./s.d90(jg));
%Dano see what happens with coarse material
else;
ucr=8.5d0.*s.d50(jg).^0.6d0.*log10(4.0d0.*hloc./s.d90(jg));
end;
% Van Thiel de Vries & Reniers 2008
elseif((par.form==form_vanthiel_vanrijn) ||(par.form==form_vanrijn1993)) ;
if(s.d50(jg)<=0.0005);
%Shields
ucrc=0.19d0.*s.d50(jg).^0.1d0.*log10(4.0d0.*hloc./s.d90(jg));
%Komar and Miller (1975)
ucrw(:)=0.24d0.*(delta.*par.g).^0.66d0.*s.d50(jg).^0.33d0.*par.trep.^0.33d0;
elseif(s.d50(jg)<=0.002) ;
%Shields
ucrc=8.5d0.*s.d50(jg).^0.6d0.*log10(4.0d0.*hloc./s.d90(jg));
%Komar and Miller (1975)
ucrw(:)=0.95d0.*(delta.*par.g).^0.57d0.*s.d50(jg).^0.43.*par.trep.^0.14;
elseif(s.d50(jg)>0.002) ;
%Maynord (1978) --> also Neill (1968) where 1.3d0 = 1.4d0
ucrc=1.3d0.*sqrt(delta.*par.g.*s.d50(jg)).*(hloc./s.d50(jg)).^(0.5d0.*onethird);
%Komar and Miller (1975)
ucrw(:)=0.95d0.*(delta.*par.g).^0.57d0.*s.d50(jg).^0.43.*par.trep.^0.14;
end;
b2 = vmg./max(vmg+sqrt(urms2),par.eps);
%Van Rijn 2007 (Bed load transport paper)
ucr = b2.*ucrc +(1-b2).*ucrw;
end;

% Lodewijk: implementation Rhee (2010), reduction sediment transport due dilatancy
srfrhee(:,:)  = 0.0d0;
srftotal(:,:) = 1.0d0;
if(par.dilatancy == 1);
% Erosion velocity, for now asume it equal to -s%dzbdt of the previous time step
vero(:,:) = max(0.0d0,-1.0d0.*s.dzbdt);
% Permeability, Adel 1987, which is based on the Ergun equation
kl = par.g./(160.0d0.*kvis).*(s.d15(jg).^2.0d0).*((par.por.^3.0d0)./(1.0d0-par.por).^2.0d0);
% bed porosity, estimated here as maximum porosity, to be added to user input
% A=3/4 for single particles and A=1/(1-n0) for a continuum
% Reduction factor on the critical Shields parameter by dilatancy (Van Rhee, 2010)
srfrhee(:,:) = vero(:,:)./kl.*(par.pormax-par.por)./(1.0d0-par.pormax).*par.rheea./delta;
end;
% Lodewijk: implementation bed slope reduction on critical flow velocity
if(par.bdslpeffini == bdslpeffini_none);
srftotal(:,:) = 1.0d0 + srfrhee(:,:);
elseif(par.bdslpeffini == bdslpeffini_total || par.bdslpeffini == bdslpeffini_bed) ;
for j=1:s.ny+1;
for i=1:s.nx+1;
% Prevent NaN values if too small values
if((abs(s.ue(i,j)) > 0.000001d0 || abs(s.ve(i,j)) > 0.000001d0) &&(abs(s.dzbdx(i,j))>0.000001d0 || abs(s.dzbdy(i,j))>0.000001d0));
% Angle between the x-axis and the flow velocity
% REMARK: also include waves in the velocity?
if(s.ue(i,j) < 0.0d0);
alpha1 = atan(s.ve(i,j)./s.ue(i,j)) + par.px;
else;
alpha1 = atan(s.ve(i,j)./s.ue(i,j));
end;
% Angle between the x-axis and the bed slope vector directed in down-slope direction, derived in thesis Lodewijk
if(s.dzbdy(i,j) >= 0.0d0);
alpha2 = -atan(s.dzbdx(i,j)./s.dzbdy(i,j))+1.5d0.*par.px;
else;
alpha2 = -atan(s.dzbdx(i,j)./s.dzbdy(i,j))+0.5d0.*par.px;
end;
% Angle between the flow direction and the up-slope directed vector
psi = alpha1-(alpha2-par.px);
% A smaller slope could result in a NaN for beta
if(abs(s.dzbdx(i,j))<0.000001d0);
% Beta purely based on dzbdy
beta = atan(abs(s.dzbdy(i,j)));
else;
% Maximum absolute bed slope angle, derived in thesis Lodewijk
beta = atan(abs(s.dzbdx(i,j)./sin(atan(s.dzbdx(i,j)./s.dzbdy(i,j)))));
end;
% Take min to exclude NaN's
beta = min(beta,phi);
if(par.dilatancy == 1);
% Soulsby (1997), modified by Lodewijk (see Thesis)
srftotal(i,j) =(cos(psi).*sin(beta)+sqrt(((srfrhee(i,j)).^2+2.*srfrhee(i,j).*cos(beta)+cos(beta).^2).*tan(phi).^2-sin(psi).^2.*sin(beta).^2))./tan(phi);
else;
% Soulsby (1997)
srftotal(i,j) =(cos(psi).*sin(beta)+sqrt(cos(beta).^2.*tan(phi).^2-sin(psi).^2.*sin(beta).^2))./tan(phi);
end;
end;
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1+1);
end;
% Calculate the new critical velocity based on the modification factors on the Shields parameter, bed slope only on bed load
% Lodewijk
ucrb(:,:) = ucr(:,:).*sqrt(srftotal);
if(par.bdslpeffini == bdslpeffini_total);
ucrs = ucrb;
else;
% Lodewijk, no effect on suspended load by bed slope
ucrs = ucr.*(1.0d0+sqrt(srfrhee));
end;


% Soulsby van Rijn
if(par.form==form_soulsby_vanrijn);
%
% drag coefficient
z0 = par.z0;
cd=(0.40d0./(log(max(hloc,10.0d0.*z0)./z0)-1.0d0)).^2;
%
% transport parameters
% bed load coefficent
asb=0.005d0.*hloc.*(s.d50(jg)./hloc./(delta.*par.g.*s.d50(jg))).^1.2d0;
% suspended load coeffient
ass=0.012d0.*s.d50(jg).*dster(jg).^(-0.6d0)./(delta.*par.g.*s.d50(jg)).^1.2d0;
%
% Make 0.018/Cd is always smaller than the flow friction coefficient
term1=(vmg.^2+0.018./cd.*par.sws.*urms2);
%
% reduce sediment suspensions for (inundation) overwash conditions with critical flow velocitties
% vmag2=min(vmag2,par%smax*par%C**2*D50(jg)*delta)
% vmag2=min(vmag2,par%smax*par%g/par%cf*D50(jg)*delta)            ! In terms of cf
% term1=sqrt(vmag2+0.018d0/Cd*urms2)     ! nearbed-velocity
% the two above lines are comment out and replaced by a limit on total velocity u2+urms2, robert 1/9 and ap 28/11
%
term1=min(term1,par.smax.*par.g./s.cf.*s.d50(jg).*delta);
term1=sqrt(term1);
%
%initialize ceqb
ceqb = 0.0d0.*term1;
%initialize ceqs
ceqs = 0.0d0.*term1;
for j=1:s.ny+1;
for i=1:s.nx;
% Lodewijk: sepperate bed load from suspended load since Ucr not anymore the same
if(term1(i,j)>ucrb(i,j) && hloc(i,j)>par.eps);
ceqb(i,j)=asb(i,j).*(term1(i,j)-ucrb(i,j)).^2.4d0;
end;
if(term1(i,j)>ucrs(i,j) && hloc(i,j)>par.eps);
ceqs(i,j)=ass.*(term1(i,j)-ucrs(i,j)).^2.4d0;
end;
end; i=fix(s.nx+1);
end; j=fix(s.ny+1+1);
%
% Van Thiel de Vries & Reniers 2008
elseif(par.form==form_vanthiel_vanrijn) ;
%
% transport parameters
%bed load coefficent
asb=0.015d0.*hloc.*(s.d50(jg)./hloc).^1.2d0./(delta.*par.g.*s.d50(jg)).^0.75d0;
%suspended load coeffient
ass=0.012d0.*s.d50(jg).*dster(jg).^(-0.6d0)./(delta.*par.g.*s.d50(jg)).^1.2d0;
%
% Jaap: par%sws to set short wave stirring to zero
% Jaap: Van Rijn use Peak orbital flow velocity --> 0.64 corresponds to 0.4 coefficient regular waves Van Rijn (2007)
term1=vmg.^2+0.64d0.*par.sws.*urms2;
% reduce sediment suspensions for (inundation) overwash conditions with critical flow velocitties
term1=min(term1,par.smax.*par.g./max(s.cf,1d-10).*s.d50(jg).*delta);
term1=sqrt(term1);
%
%initialize ceqb
ceqb = 0.0d0.*term1;
%initialize ceqs
ceqs = 0.0d0.*term1;
for j=1:s.ny+1;
for i=1:s.nx;
% Lodewijk: sepperate bed load from suspended load since Ucr not anymore the same
if(term1(i,j)>ucrb(i,j) && hloc(i,j)>par.eps);
ceqb(i,j)=asb(i,j).*(term1(i,j)-ucrb(i,j)).^1.5;
end;
if(term1(i,j)>ucrs(i,j) && hloc(i,j)>par.eps);
ceqs(i,j)=ass.*(term1(i,j)-ucrs(i,j)).^2.4;
end;
end; i=fix(s.nx+1);
end; j=fix(s.ny+1+1);
% Van Rijn (1993) a la Delft3D
elseif(par.form==form_vanrijn1993) ;
% General parameters
% similar to Van Thiel - van Rijn -> effective velocity
ue          = sqrt(vmg.^2+0.64d0.*par.sws.*urms2);
% velocity used for sediment transport
used        = vmg + s.ua;
%
% Bed load transport (Delft3D-FLOW manual; eq 11.82)
%
for j=1:s.ny+1;
for i=1:s.nx+1;
% Ucrs is weighted current (including waves)
if(ue(i,j)>ucrb(i,j) && hloc(i,j)>par.eps);
% Effective mobility parameter
me              =(ue(i,j)-ucrb(i,j)).^2 ./((delta.*par.g.*s.d50(jg)));
% 'normal' mobility parameter
m               =(ue(i,j)-0).^2 ./((delta.*par.g.*s.d50(jg)));
% bed load coefficent (without ps)
asb(i,j)        = 0.006d0 .* s.d50(jg) .* ws0(jg);
% calculate the main part
ceqb(i,j)       = asb(i,j) .*(m.^0.5).*(me.^0.7);
% goal is sediment concentration: divide by velocity
ceqb(i,j)       = ceqb(i,j) ./ used(i,j);
else;
ceqb(i,j)       = 0;
end;
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1+1);
%
% Suspended transport (-> solving the concentration profile with constant diffu
%
% orbital velocity (below 11.80)
uorb(:)        = par.px.*s.h./par.trep./sinh(min(max(s.k,0.01d0).*max(s.hh,par.delta.*s.h),10.0d0));
% peak orbital excursion (?)
a           = uorb.*par.trep ./(2d0.*par.px);
% wave-related grain friction friction (11.77)
fw          = exp(-6 + 5.2d0 .* min((a./par.d90(jg)).^-0.019d0, 1d-10));
% efficiency factor waves (11.75)
uw(:)          = max(0.063d0,(((1.5d0-(s.h./s.hh)).^2) ./ 8));
% bed shear stress due to waves (11.76)
taubw       = 0.25d0 .* par.rho .* fw .* uorb .^ 2;
% efficiency currents = 1, tau,cw is XBeach (eq. 11.70)
tautotal    =(s.taubx.^2 + s.tauby.^2).^0.5 + uw.*fw;
% critical bed shear stress (Soulsby, 1997 at Eq. 74)
taucr(:)       = shieldscrit(jg) .* par.g .*(par.rhos-par.rho) .* s.d50(jg);
%
for j=1:s.ny+1;
for i=1:s.nx+1;
if(tautotal(i,j)>taucr(i,j) && hloc(i,j)>par.eps);
% 'normal' mobility parameter
m       =(ue(i,j)-0).^2 ./((delta.*par.g.*s.d50(jg)));
% kscr (5c bed load Van Rijn)
kscr    = 0.5d0  .*(182.5-0.652.*m).*par.d50(jg);
% limits, min
kscr    = min(20.*par.d50(jg), kscr);
% limits, max
kscr    = max(150.*par.d50(jg), kscr);
% reference level based on 0.5 * kscr
refa    = 0.5.*kscr;
ca      = 0.0015 .*(s.d50(jg) ./ refa) .*((tautotal(i,j)-taucr(i,j))./taucr(i,j)).^1.5 ./ dster(jg).^0.3;
% near-bed sediment mixing coefficient (simple assumption)
e       = 0.004 .* dster(jg) .* uorb(i,j) .* 6 .* refa;
ceqs(i,j) = ca .* e .*( exp(-ws0(jg) .* refa ./ e) - exp(-ws0(jg) .* s.hh(i,j) ./ e) ) ./ ws0(jg);
else;
ceqs(i,j) = 0;
end;
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1+1);
end;
%
% maximum equilibrium bed concentration
ceqb = min(ceqb./hloc,par.cmax./2);
s.ceqbg(:,:,jg) =(1-par.bulk).*ceqb.*s.sedcal(jg).*s.wetz;
% maximum equilibrium suspended concentration
ceqs = min(ceqs./hloc,par.cmax./2);
s.ceqsg(:,:,jg) =(ceqs+par.bulk.*ceqb).*s.sedcal(jg).*s.wetz;
% Jaap: old brute method to prevent strong coastline erosion
% where (hloc<=par%hmin) ceqsg(:,:,jg) = 0.0d0
% end of grain size classes
end; jg =fix(par.ngd+1);
% end of grain size classes

end %subroutine sedtransform

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function [s,par]=waveturb(s,par,varargin);
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University !
% Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski,            !
% Jaap van Thiel de Vries, Robert McCall                                  !
%                                                                         !
% d.roelvink@unesco-ihe.org                                               !
% UNESCO-IHE Institute for Water Education                                !
% P.O. Box 3015                                                           !
% 2601 DA Delft                                                           !
% The Netherlands                                                         !
%                                                                         !
% This library is free software; you can redistribute it and/or           !
% modify it under the terms of the GNU Lesser General Public              !
% License as published by the Free Software Foundation; either            !
% version 2.1 of the License, or (at your option) any later version.      !
%                                                                         !
% This library is distributed in the hope that it will be useful,         !
% but WITHOUT ANY WARRANTY; without even the implied warranty of          !
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU        !
% Lesser General Public License for more details.                         !
%                                                                         !
% You should have received a copy of the GNU Lesser General Public        !
% License along with this library; if not, write to the Free Software     !
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307     !
% USA                                                                     !
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
persistent disturb dzsdt_cr i j ksource kturbu kturbv ml sturbu sturbv twothird ; 

 params;
 spaceparams;
 xmpi_module;



 if isempty(i), i=0; end;
 if isempty(j), j=0; end;
 if isempty(ml), ml=0; end;
 if isempty(disturb), disturb=0; end;
 if isempty(twothird), twothird=0; end;
 if isempty(ksource), ksource=[]; end;
 if isempty(kturbu), kturbu=[]; end;
 if isempty(kturbv), kturbv=[]; end;
 if isempty(sturbu), sturbu=[]; end;
 if isempty(sturbv), sturbv=[]; end;
 if isempty(dzsdt_cr), dzsdt_cr=[]; end;

%include 's.ind'
%include 's.inp'

if(~ ~isempty(kturbu));
ksource=zeros(s.nx+1,s.ny+1);
kturbu=zeros(s.nx+1,s.ny+1);
kturbv=zeros(s.nx+1,s.ny+1);
sturbu=zeros(s.nx+1,s.ny+1);
sturbv=zeros(s.nx+1,s.ny+1);
dzsdt_cr=zeros(s.nx+1,s.ny+1);
twothird=2.0d0./3.0d0;
end;
% use lagrangian velocities
%Jaap
kturbu(:)       = 0.0d0;
%Jaap
kturbv(:)       = 0.0d0;
%     dzsdt_cr=par%beta*s%c
dzsdt_cr(:)=par.beta.*sqrt(par.g.*s.hh);
% Update roller thickness
%s%rolthick=s%rolthick+par%dt*(abs(s%dzsdt)-dzsdt_cr) Dano: not abs!
s.rolthick=s.rolthick+par.dt.*(s.dzsdt-dzsdt_cr);
s.rolthick=max(s.rolthick,0.0d0);
s.rolthick=min(s.rolthick,s.hh);
fmask=find((s.wetz==0));NOTfmask=find(~(s.wetz==0));
s.rolthick=0.0d0;
%end where;

% Jaap compute sources and sinks
% long wave source
ksource(:) = 0;
if(par.lwt == 1);
%only important in shallow water, where s%c=sqrt(gh)
ksource(:) = par.g.*s.rolthick.*par.beta.*sqrt(par.g.*s.hh);
end;
if(par.turbadv == turbadv_none);
% See Battjes, 1975 / 1985
s.kturb =(s.dr./par.rho).^twothird;
else;
ksource = ksource + s.dr./par.rho;


% Jaap do long wave turb approach for both short waves and long waves
%
% Turbulence in uu-points
%
for j=1:s.ny+1;
for i=1:s.nx;
if(s.uu(i,j)>0.0d0);
kturbu(i,j)=par.thetanum.*s.kturb(i,j)+(1.0d0-par.thetanum).*s.kturb(min(i+1,s.nx),j);
elseif(s.uu(i,j)<0.0d0) ;
kturbu(i,j)=par.thetanum.*s.kturb(i+1,j)+(1.0d0-par.thetanum).*s.kturb(max(i,2),j);
else;
kturbu(i,j)=0.5d0.*(s.kturb(i,j)+s.kturb(i+1,j));
end;
end; i=fix(s.nx+1);
end; j=fix(s.ny+1+1);
if(xmpi_isbot);
 kturbu(s.nx+1,:) = s.kturb(s.nx+1,:);
end;
%
% Turbulence in vv-points
%
for j=1:s.ny;
for i=1:s.nx+1;
if(s.vv(i,j)>0);
kturbv(i,j)=par.thetanum.*s.kturb(i,j)+(1.0d0-par.thetanum).*s.kturb(i,min(j+1,s.ny));
elseif(s.vv(i,j)<0) ;
kturbv(i,j)=par.thetanum.*s.kturb(i,j+1)+(1.0d0-par.thetanum).*s.kturb(i,max(j,2));
else;
kturbv(i,j)=0.5d0.*(kturbv(i,j)+kturbv(i,j+1));
end;
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1);
%Robert
kturbv(:,s.ny+1) = s.kturb(:,s.ny+1);
%
% Turbulence advection in X and Y direction
%
if(par.turbadv == turbadv_lagrangian);
sturbu=kturbu.*s.uu.*s.hu.*s.wetu;
sturbv=kturbv.*s.vv.*s.hv.*s.wetv;
elseif(par.turbadv == turbadv_eulerian) ;
sturbu=kturbu.*s.ueu.*s.hu.*s.wetu;
sturbv=kturbv.*s.vev.*s.hv.*s.wetv;
end;
%
% Update turbulence
%
if(s.ny>0);
for j=2:s.ny+1;
for i=2:s.nx+1;
if(par.betad>0);
disturb=par.betad.*s.kturb(i,j).^1.5d0;
else;
ml=max(s.rolthick(i,j),0.07d0.*max(s.hh(i,j),par.hmin));
disturb=0.08d0.*max(s.hh(i,j),par.hmin)./ml.*s.kturb(i,j).^1.5d0;
end;
s.kturb(i,j) =(s.hold(i,j).*s.kturb(i,j)-par.dt.*((sturbu(i,j).*s.dnu(i,j)-sturbu(i-1,j).*s.dnu(i-1,j)+sturbv(i,j).*s.dsv(i,j)-sturbv(i,j-1).*s.dsv(i,j-1)).*s.dsdnzi(i,j)-(ksource(i,j)-disturb)))./max(s.hh(i,j),par.hmin);
s.kturb(i,j)=max(s.kturb(i,j),0.0d0);

end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1+1);
else;
j=1;
for i=2:s.nx+1;

s.kturb(i,j) =(s.hold(i,j).*s.kturb(i,j)-par.dt.*((sturbu(i,j).*s.dnu(i,j)-sturbu(i-1,j).*s.dnu(i-1,j)).*s.dsdnzi(i,j)-(ksource(i,j)-par.betad.*s.kturb(i,j).^1.5d0)))./max(s.hh(i,j),par.hmin);
s.kturb(i,j)=max(s.kturb(i,j),0.0d0);

end; i=fix(s.nx+1+1);
end;

s.kturb = s.kturb./max(s.hh,0.01d0);

% Jaap only required for advection mode?
if(xmpi_istop);
 s.kturb(1,:)=s.kturb(2,:);
end;
if(xmpi_isbot);
 s.kturb(s.nx+1,:)=s.kturb(s.nx+1-1,:);
end;
if(s.ny>0);
if(xmpi_isleft);
  s.kturb(:,1)=s.kturb(:,2);
end;
if(xmpi_isright);
 s.kturb(:,s.ny+1)=s.kturb(:,s.ny+1-1);
end;
end;

#ifdef usempi;
[s.kturb]=xmpi_shift_ee(s.kturb);
#end;

% turbadv == 'none'
end

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


function [s,par]=rvr(s,par,varargin);

persistent alpha_ml b1 beta bm m1 m2 m3 m4 m5 m6 urs ; 

 params;
 spaceparams;
 xmpi_module;



 if isempty(m1), m1=0; end;
 if isempty(m2), m2=0; end;
 if isempty(m3), m3=0; end;
 if isempty(m4), m4=0; end;
 if isempty(m5), m5=0; end;
 if isempty(m6), m6=0; end;
 if isempty(alpha_ml), alpha_ml=0; end;
 if isempty(beta), beta=0; end;

 if isempty(urs), urs=[]; end;
 if isempty(bm), bm=[]; end;
 if isempty(b1), b1=[]; end;

%include 's.ind'
%include 's.inp'

% only in first timestep..
if(~ ~isempty(urs));

urs=zeros(s.nx+1,s.ny+1);
bm=zeros(s.nx+1,s.ny+1);
b1=zeros(s.nx+1,s.ny+1);

% a = 0
m1 = 0;
% b = 0.79 +/- 0.023
m2 = 0.7939;
% c = -0.61 +/- 0.041
m3 = -0.6065;
% d = -0.35 +/- 0.032
m4 = 0.3539;
% e = 0.64 +/- 0.025
m5 = 0.6373;
% f = 0.60 +/- 0.043
m6 = 0.5995;
alpha_ml = -log10(exp(1.0d0))./m4;
beta  = exp(m3./m4);

end;

%Ursell number
urs(:) = 3.0d0./8.0d0.*sqrt(2.0d0).*s.h.*s.k./(s.k.*s.hh).^3;
urs = max(urs,0.000000000001d0);
%Boltzmann sigmoid (eq 6)
bm = m1 +(m2-m1)./(1.0d0+beta.*urs.^alpha_ml);
b1 =(-90.0d0+90.0d0.*tanh(m5./urs.^m6)).*par.px./180.0d0;
%Skewness (eq 8)
s.sk = bm.*cos(b1);
%Asymmetry(eq 9)
s.as = bm.*sin(b1);
s.ua = par.sws.*(par.facsk.*s.sk-par.facas.*s.as).*s.urms;

% multiply Sk and As with wetz to get zeros at h = 0 for output
s.sk = s.sk.*s.wetz;
s.as = s.as.*s.wetz;

end %subroutine rvr

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function [s,par]=vt(s,par,varargin);

persistent detadxmax detadxmean dh dt duddtmax duddtmean dudtmax dudtmean f0 f1 f2 f3 h0 i ih0 ih1 it0 it1 j nh nt p q siguref t0 t0fac ; 

 params;
 spaceparams;
 readkey_module;
 xmpi_module;



 if isempty(i), i=0; end;
 if isempty(j), j=0; end;
 if isempty(nh), nh=0; end;
 if isempty(nt), nt=0; end;
 if isempty(ih0), ih0=0; end;
 if isempty(it0), it0=0; end;
 if isempty(ih1), ih1=0; end;
 if isempty(it1), it1=0; end;
 if isempty(p), p=0; end;
 if isempty(q), q=0; end;
 if isempty(f0), f0=0; end;
 if isempty(f1), f1=0; end;
 if isempty(f2), f2=0; end;
 if isempty(f3), f3=0; end;
 if isempty(t0fac), t0fac=0; end;
 if isempty(siguref), siguref=0; end;
 if isempty(duddtmax), duddtmax=0; end;
 if isempty(dudtmax), dudtmax=0; end;
 if isempty(duddtmean), duddtmean=0; end;
 if isempty(dudtmean), dudtmean=0; end;
 if isempty(detadxmean), detadxmean=0; end;
 if isempty(dh), dh=0; end;
 if isempty(dt), dt=0; end;

 if isempty(h0), h0=[]; end;
 if isempty(t0), t0=[]; end;
 if isempty(detadxmax), detadxmax=[]; end;
% Robert: RF table now included in source code, rather than read from file
% Rienecker Fenton table with amongst others amplitudes non-linear components obtained with stream function theory

% Robert: 'RF.inc' contains definition of RF as real*8(5,33,40) with 'parameter' attribute
%   so RF values may not be modified! To save memory, only rows 13,14,15,16 and 18 of the
%   original matrix are stored. So new row 1 corresponds with old row 13, etc.

%include 's.ind'
%include 's.inp'


% only in first timestep..
if(~ ~isempty(h0));
h0=zeros(s.nx+1,s.ny+1);
t0=zeros(s.nx+1,s.ny+1);
detadxmax=zeros(s.nx+1,s.ny+1);
dh = 0.03d0;
dt = 1.25d0;
nh = fix(floor(0.99d0./dh));
nt = fix(floor(50.0d0./dt));
end;

% non-linearity of short waves is listed in table as function of dimensionless wave height h0 and dimensionless wave period t0

% compute dimensionless wave height and wave period in each grid point..
h0(:) = min(nh.*dh,max(dh,min(s.h,s.hh)./s.hh));
t0(:) = min(nt.*dt,max(dt,par.trep.*sqrt(par.g./s.hh)));

% estimate Sk, As and ua by interpolating table values
for j=1:s.ny+1;
for i=1:s.nx+1;
% interpolate table values....
ih0=fix(floor(h0(i,j)./dh));
it0=fix(floor(t0(i,j)./dt));
ih1=fix(min(ih0+1,nh));
it1=fix(min(it0+1,nt));
p=(h0(i,j)-ih0.*dh)./dh;
q=(t0(i,j)-it0.*dt)./dt;

f0=(1-p).*(1-q);
f1=p.*(1-q);
f2=q.*(1-p);
f3=p.*q;

% Skewness and assymetry
s.sk(i,j) = f0.*rf(1,ih0,it0)+f1.*rf(1,ih1,it0)+ f2.*rf(1,ih0,it1)+f3.*rf(1,ih1,it1);
s.as(i,j) = f0.*rf(2,ih0,it0)+f1.*rf(2,ih1,it0)+ f2.*rf(2,ih0,it1)+f3.*rf(2,ih1,it1);

% Sediment advection velocity from Skewness and Assymetry
% ua(i,j) = par%sws*par%facua*(Sk(i,j)-As(i,j))*urms(i,j)
s.ua(i,j) = par.sws.*(par.facsk.*s.sk(i,j)-par.facas.*s.as(i,j)).*s.urms(i,j);

% Estimate bore period Tbore and mean slope bore front to feeded back in roller energy balance

% correct slope in case 1.25>T0>50
if(t0(i,j)==50.0d0);
t0fac = 50.0d0./max((par.trep.*sqrt(par.g./s.hh(i,j))),50.0d0);
elseif(t0(i,j)==1.25);
t0fac = 1.25d0./min((par.trep.*sqrt(par.g./s.hh(i,j))),1.25d0);
else;
t0fac = 1.0d0;
end;

% detadxmax for Tbore...
% dimnesionless maximum acceleration under bore front
duddtmax = f0.*rf(3,ih0,it0)+f1.*rf(3,ih1,it0)+ f2.*rf(3,ih0,it1)+f3.*rf(3,ih1,it1);
siguref = f0.*rf(4,ih0,it0)+f1.*rf(4,ih1,it0)+ f2.*rf(4,ih0,it1)+f3.*rf(4,ih1,it1);
% translate dimensionless duddtmax to real world dudtmax
%         /scale with variance and go from [-] to [m/s^2]     /tableb./dimensionless dudtmax
dudtmax = s.urms(i,j)./max(par.eps,siguref).*sqrt(par.g./s.hh(i,j)).*t0fac.*duddtmax;
detadxmax(i,j) = dudtmax.*sinh(s.k(i,j).*s.hh(i,j))./max(max(s.c(i,j),sqrt(s.h(i,j).*par.g)),1d-10)./s.sigm(i,j);

% detadxmean for roller energy balance dissipation...
if(par.rfb==1);
duddtmean = f0.*rf(5,ih0,it0)+f1.*rf(5,ih1,it0)+ f2.*rf(5,ih0,it1)+f3.*rf(5,ih1,it1);
dudtmean = s.urms(i,j)./max(par.eps,siguref).*sqrt(par.g./s.hh(i,j)).*t0fac.*duddtmean;
detadxmean = dudtmean.*sinh(s.k(i,j).*s.hh(i,j))./max(s.c(i,j),sqrt(s.h(i,j).*par.g))./s.sigm(i,j);
s.br(i,j) = par.brfac.*sin(atan(detadxmean));
end;

end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1+1);

s.tbore = max(par.trep./25.0d0,min(par.trep./4.0d0,s.h./(max(max(s.c,sqrt(s.h.*par.g)),1d-10).*max(detadxmax,par.eps))));
s.tbore = par.tbfac.*s.tbore;

% multiply Sk and As with wetz to get zeros at h = 0 for output
s.sk = s.sk.*s.wetz;
s.as = s.as.*s.wetz;

end %subroutine vt

function [s, par]=setbathy_update(s, par,varargin);

persistent dummy i j zbnew ; 

 params;
 spaceparams;
 interp;



 if isempty(i), i=0; end;
 if isempty(j), j=0; end;
 if isempty(dummy), dummy=0; end;
 if isempty(zbnew), zbnew=zeros(s.nx+1,s.ny+1); end;

% interpolate from file
for j=1:s.ny+1;
for i=1:s.nx+1;
[s.tsetbathy,s.setbathy(i,j,:),par.nsetbathy,par.t,zbnew(i,j),dummy]=linear_interp(s.tsetbathy,s.setbathy(i,j,:),par.nsetbathy,par.t,zbnew(i,j),dummy);
end; i=fix(s.nx+1+1);
end; j=fix(s.ny+1+1);
% update water level
s.zs = s.zs+zbnew-s.zb;
% update bed level
s.zb = zbnew;
% update wet and dry cells
fmask=find((s.zs<s.zb+par.eps));NOTfmask=find(~(s.zs<s.zb+par.eps));
s.wetz=0;
s.zs = s.zb+par.eps;
s.hh = par.eps;
%end where;

end %subroutine setbathy_update

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function [s,par]=hybrid(s,par,varargin);

persistent first hav1d i indx irrb j j1 maxiter niter runup_old slopeind ; 

 params;
 interp;
 spaceparams;
 xmpi_module;



 if isempty(i), i=0; end;
 if isempty(j), j=0; end;
 if isempty(j1), j1=0; end;
 if isempty(indx), indx=0; end;
 if isempty(first), first=0; end;
 if isempty(niter), niter=0; end;
 if isempty(maxiter), maxiter=0; end;
 if isempty(slopeind), slopeind=[]; end;
 if isempty(hav1d), hav1d=[]; end;
 if isempty(irrb), irrb=0; end;
 if isempty(runup_old), runup_old=0; end;

%include 's.ind'
%include 's.inp'

if(~ ~isempty(hav1d));
hav1d=zeros(1,s.nx+1);
slopeind=(zeros(1,s.nx+1));
end;

for j=1:s.ny+1;
indx = fix(s.nx+1);
first = 0;
for i=1:s.nx;
% transition from wet to dry
if(s.wetz(i,j)-s.wetz(max(i-1,1),j)==-1 && first==0);
% only consider first dry point
first = 1;
% find wave height for runup at L1 meter from water line
[s.xz(:,j),s.h(:,j),dumvar3,dumvar4,s.hrunup(j),indx]=linear_interp(s.xz(:,j),s.h(:,j),s.nx+1,s.xz(i-1,j)-s.l1(i-1,j),s.hrunup(j),indx);
% Find toe of runup slope if present (dzbdx > 0.15).
% If not present Hrunup will converge to H at the water line (where H = 0 per definition)
for j1=indx:i-1;
% TODO: Strange condition. In case of no toe, or no steep slope,
%   there will still be extra turbulence at L1 meter from the water line...
% cross shore location structure toe
if(s.dzbdx(j1,j)<0.15d0 || s.structdepth(j1,j)>0.1d0);
indx = fix(j1);
end;
end; j1=fix(i-1+1);
% update Hrunup and runup x-location
% short wave height at revetment toe
s.hrunup(j) = s.h(indx,j);
% cross-shore location revetment toe
s.xhrunup(j) = s.xz(indx,j);
% cross-shore index revetment toe
s.istruct(j) = indx.*1.0d0;
% cross-shore location waterline (inlcuding lw-s%runup)
s.iwl(j) =(i-1).*1.0d0;

% now iteratively compute runup
hav1d = s.hh(:,j);
runup_old = realmax;
s.runup(j) = 0;
niter = 0;
maxiter = 50;
while(abs(s.runup(j)-runup_old)>0.01d0 & niter < maxiter);
niter = fix(niter +1);
runup_old = s.runup(j);
slopeind(:) = 0;
fmask=find((hav1d>par.eps & s.dzbdx(:,j)>0.15d0));NOTfmask=find(~(hav1d>par.eps & s.dzbdx(:,j)>0.15d0));
slopeind(fmask) = 1;
%end where;
s.strucslope(j) = sum(sum(s.dzbdx(indx:s.nx,j).*s.dsu(indx:s.nx,j).*slopeind([indx:s.nx])))./max(par.eps,sum(sum(s.dsu(indx:s.nx,j).*slopeind([indx:s.nx]))));
if(s.strucslope(j) > 0.0d0);
irrb = s.strucslope(j)./sqrt(2.*par.px.*max(s.hrunup(j),par.eps)./par.g./par.trep.^2);
s.runup(j) = par.facrun.*min(irrb,2.3d0).*s.hrunup(j).*cos(2.*par.px./par.trep.*par.t);
else;
s.runup(j) = 0.0d0;
end;
% This triggers s%runup to be calculated for the complete hinterland (also behind a dune).
% Only values calculated for the first dune front are used in the avalanching algorithm (morphevolution)
hav1d =  s.wetz(:,j).*max(par.eps,(s.hh(:,j) + s.runup(j))) +(1.0d0-s.wetz(:,j)).*max(par.eps,s.runup(j)+s.zs(i-1,j)-s.zb(:,j) );
end;
end;
end; i=fix(s.nx+1);
end; j=fix(s.ny+1+1);



end %subroutine hybrid

