c fluvial_erosion
subroutine fluvial_erosion (xkf,xlf_BR,xlf_AL,width_c,
& x,y,h,h0,hi,ndon,nnode,
& surface,slope,length,
& water,sediment,
& ibucket,dt,fix,
& dhh,
& nb,nn,nbmax,
& iorder,itype_node,
& dhminfluv,dhmaxfluv,
& nlake,
& sea_level,outflux,ideposition)
c subroutine to calculate fluvial erosion
c INPUT: xkf = diffusivity
c xlf_BR = bedrock erosion length scale
c xlf_AL = alluvials erosion length scale
c width_c = coef controling channel width as fnct of discharge
c x,y = x- and y-nodal positions
c h = present topography
c h0 = bedrock-alluvion interface
c hi = initial topography (at time 0)
c ndon = donor array
c nnode = number of nodes
c surface = surface associated with each node
c slope = slope associated with each nodal link (stream)
c length = length associated with each nodal link (stream)
c ibucket = working array
c dt = time step length
c fix = boundary conditoin array
c dhh = amount of material eroded/deposited over the time step
c nb = number of natural neighbours per node
c nn = list of natural neighbours
c nbmax = maximum of natural neighbours
c iorder = working array containing the proper sequence in node
c number in which to perform the river erosion operations
c itype_node = type associated to each node (see later in code)
c dhminfluv = minimum amount removed by river erosion
c dhmaxfluv = maximum amount removed by river erosion
c nlake = determines whether a node is part of a lake or not
c sea_level = sea level
c ideposition = to prevent deposition by rivers (=0)
c OUTPUT: several arrays are updated:
c h = new current topography
c the contribution from river incision to outflux is calculated
c the following arrays are filled:
c water = water discharge at each point
c sediment = sediment load at each point
c subroutines called:
c - debug
c - find_order
common /vocal/ ivocal
real x(nnode),y(nnode),h(nnode),h0(nnode),hi(nnode)
real xkf(nnode),xlf_BR(nnode)
real*4 surface(nnode),dhh(nnode)
real*4 water(nnode),sediment(nnode)
real*4 slope(nnode),length(nnode)
integer ndon(nnode)
integer*4 ibucket(nnode)
real fix(nnode)
integer*4 nb(*),nn(nbmax,*)
integer*4 iorder(nnode),itype_node(nnode)
integer*4 nlake(nnode)
c sets all ibucket to 0
c sets all water to 1 and all sediment to 0
c note that water(i) can be different for every node; this is where a more
c "complex" orographic model should be included...
do i=1,nnode
ibucket(i)=0
sediment(i)=0.
dhh(i)=0.
enddo
dhminfluv=0.
dhmaxfluv=0.
c finds proper ordering (cascade algorithm)
if (ivocal.eq.1) call debug ('find_order$',0)
call find_order (ibucket,ndon,fix,iorder,
& nnode,norder)
if (ivocal.eq.1) call debug ('fluvial_erosion$',1)
c itype_node determines the type of node:
c -1: local minima next to (or on) a fixed boundary
c 0: diffusion only
c 1: channel because one of its parents was a channel
do i=1,nnode
itype_node(i)=0
enddo
do jorder=1,norder
i=iorder(jorder)
slop=-slope(i)
c special treatment for self donors and b.c.nodes
if (ndon(i).eq.i) then
dh=0.
itype_node(i)=-2
outflux=outflux+sediment(i)*(1.-fix(i))
sediment(i)=0.
water(i)=0.
c special treatment for nodes below sea level
elseif (h(i).lt.sea_level) then
dh=sediment(i)/surface(i)
h(i)=h(i)+dh
itype_node(i)=-1
sediment(i)=0.
else
if (itype_node(i).ne.0) itype_node(i)=1
c sedeqb is how much sediment the river can carry
sedeqb=xkf(i)*slop*water(i)*dt
c set channel width according to sqrt function of discharge
width_ch=width_c*sqrt(water(i))
c special treatment for lake nodes
if (nlake(i).eq.1) then
dh=sediment(i)/surface(i)
h(i)=h(i)+dh
sediment(i)=0.
c deposition
elseif (sediment(i).ge.sedeqb.and.ideposition.eq.1) then
dh=(sediment(i)-sedeqb)/surface(i)
c there is a maximum amount of sediment that can be dumped at one location.
c the maximum is given by the maximum height that the dumping may produce
c such that the slope it creates between this node and its donor is not
c greater than the erosional treshold. In other words, it is assumed that
c by deposition one cannot create a levee so big that at the next time step
c it is going to be eroded down
if (water(i).ne.0.) then
dhmax=sediment(i)*length(i)/xkf(i)/water(i)/dt+
& h(ndon(i))-h(i)
else
dhmax=dh
endif
if (dh.le.dhmax) then
sediment(i)=sedeqb
h(i)=h(i)+dh
else
dh=dhmax
sediment(i)=sediment(i)-dhmax*surface(i)
h(i)=h(i)+dh
endif
c erosion
else
c first case: eroding bedrock only
if (h(i).le.h0(i)) then
dsed_ch=(sedeqb)*(length(i)/xlf_BR(i))
dh=-(dsed_ch)/(width_ch*length(i))
dsediment=-dh*surface(i)
h(i)=h(i)+dh
sediment(i)=sediment(i)+dsediment
c second case: eroding alluvions only
else
dsed_ch=(sedeqb)*(length(i)/xlf_AL)
dh=-(dsed_ch)/(width_ch*length(i))
if (dh.ge.h0(i)-h(i)) then
dsediment=-dh*surface(i)
h(i)=h(i)+dh
sediment(i)=sediment(i)+dsediment
else
c third case: eroding alluvions and bedrock
dh1=h(i)-h0(i)
dsed_al_ch=dh1*length(i)*width_ch
dsed_ch=sedeqb*(length(i)/xlf_BR(i))
c calculate the height change of the bedrock interface.
dh=-(dsed_ch)/(width_ch*length(i))
dsediment=-dh*surface(i)
dsed_al_cell=dh1*surface(i)
h(i)=h0(i)+dh
sediment(i)=sediment(i)+dsediment+dsed_al_cell
endif
endif
endif
dhminfluv=amin1(dhminfluv,dh)
dhmaxfluv=amax1(dhmaxfluv,dh)
c from water+sediment+slope update height and sediment
if (ndon(i).ne.i) then
water(ndon(i))=water(ndon(i))+water(i)
sediment(ndon(i))=sediment(ndon(i))+sediment(i)
if (slop.ne.0.) itype_node(ndon(i))=1
endif
endif
c compute dhh
dhh(i)=fix(i)*dh
enddo
return
end