Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
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)
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
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