The TELEMAC-MASCARET system  trunk
bed1_suspension_erode.f
Go to the documentation of this file.
1 ! ********************************
2  SUBROUTINE bed1_suspension_erode
3 ! ********************************
4 !
5 !***********************************************************************
6 ! GAIA
7 !***********************************************************************
8 !
12 !
13 !
14 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
15 !
16  USE bief
19  IMPLICIT NONE
20 !
21 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
22 !
23  INTEGER IPOIN
24  INTEGER ILAYER,ISAND,IMUD,ICLA
25  INTEGER NSAND_VIRTUAL
26  DOUBLE PRECISION CHECK_POSITIF_MUD_TOT1,CHECK_POSITIF_MUD_TOT2
27  DOUBLE PRECISION CHECK_POSITIF_MUD,CHECK_POSITIF_SAND
28  DOUBLE PRECISION FLUER_MIX_TOT
29  DOUBLE PRECISION INTERPOL
30 !
31  IF(nsand.EQ.0) THEN
32  nsand_virtual = 1
33 ! ON FAIT CELA POUR PASSER UNE FOIS DANS LA BOUCLE SABLE MEME S IL N Y EN PAS
34 ! CELA NE POSE PAS DE PB CAR RATIO_MUD_SAND=1
35  DO ipoin = 1,npoin
36  DO ilayer = 1,nomblay
37  ratio_sand(nsand_virtual,ilayer,ipoin)= 1.d0
38  ENDDO
39  ENDDO
40  ELSE
41  nsand_virtual = nsand
42  ENDIF
43 ! INITIALIZATION OF PHYSICAL PARAMETERS OF THE ACTIVE LAYER IF MUD
44  IF(hirano.AND.nmud.GE.1) THEN
45  DO ipoin = 1,npoin
46 
47 ! UPDATE CRITICAL SHEAR STRESS OF MUD IN ACTIV LAYER IN FUNCTION OF CONCENTRATION OF MUD IN THE ACTIVE LAYER
48 ! SIMPLE INTERPOLATION BETWEEN TOCE_MUD DISCRETIZATION GIVEN BY USER
49  IF(conc_mud(1,ipoin).LE.conc_mud(2,ipoin))THEN
50  conc_mud(1,ipoin)=conc_mud(2,ipoin)
51  toce_mud(1,ipoin)=toce_mud(2,ipoin)
52  partheniades(1,ipoin)=partheniades(2,ipoin)
53  ELSEIF(conc_mud(1,ipoin).GE.
54  & conc_mud(nomblay,ipoin))THEN
55  conc_mud(1,ipoin)=conc_mud(nomblay,ipoin)
56  toce_mud(1,ipoin)=toce_mud(nomblay,ipoin)
57  partheniades(1,ipoin)=partheniades(nomblay,ipoin)
58  ELSE
59  DO ilayer=2,nomblay
60  IF(conc_mud_activ_layer(ipoin).LE.
61  & conc_mud(ilayer,ipoin))THEN
62  interpol=(conc_mud_activ_layer(ipoin)-
63  & conc_mud(ilayer-1,ipoin))/
64  & (conc_mud(ilayer,ipoin)-
65  & conc_mud(ilayer-1,ipoin))
66 !
67  toce_mud(1,ipoin)=interpol*
68  & (toce_mud(ilayer,ipoin)
69  & -toce_mud(ilayer-1,ipoin))
70  & +toce_mud(ilayer-1,ipoin)
71 !
72  partheniades(1,ipoin)=interpol*
73  & (partheniades(ilayer,ipoin)
74  & -partheniades(ilayer-1,ipoin))
75  & +partheniades(ilayer-1,ipoin)
76  GOTO 200
77  ENDIF
78  ENDDO !END LOOP NOMBLAY
79 200 CONTINUE
80  ENDIF
81 
82  ENDDO !END LOOP NPOIN
83  ENDIF !ENDIF HIRANO.AND.NMUD.GE.1
84 !
85 ! INITIALIZATION OF SHEAR STRESS OF MIXTURE AND EQUILIBRIUM CONCENTRATION FOR EACH LAYER
86 !
87  DO ilayer = 1,nomblay
88 !
89  DO isand = 1,nsand_virtual
90 ! TOCE_SAND(ISAND,IPOIN) IS COMPUTED IN INIT_SEDIMENT_GAIA
92 ! CE N'EST PAS FAIT ACTUELLMENT DANS SISYPHE POUR LA SUSPENSION
93  DO ipoin=1,npoin
94  IF(ratio_mud_sand(ilayer,ipoin).LE.0.3d0)THEN
95  toce_mix(isand,ilayer,ipoin)= toce_sand(isand,ipoin)
96  ratio_toce%ADR(num_isand_icla(isand))%P%R(ipoin)= 1.d0
97  ELSEIF(ratio_mud_sand(ilayer,ipoin).GE.0.5d0)THEN
98  toce_mix(isand,ilayer,ipoin) = toce_mud(ilayer,ipoin)
99 !
100  IF(nsand.NE.0)THEN
101  ratio_toce%ADR(num_isand_icla(isand))%P%R(ipoin)=
102  & 1.d0
103 !-----------not use in this case and Fix at 1
104  ENDIF
105 !
106  ELSE
107  toce_mix(isand,ilayer,ipoin) =
108  & (ratio_mud_sand(ilayer,ipoin)-0.3d0)/(0.5d0-0.3d0)
109  & *(toce_mud(ilayer,ipoin) - toce_sand(isand,ipoin))
110  & + toce_sand(isand,ipoin)
111 !
112  ratio_toce%ADR(num_isand_icla(isand))%P%R(ipoin)=
113  & toce_mix(isand,ilayer,ipoin)/toce_sand(isand,ipoin)
114  ENDIF
115 !!!!!!!!!!!!warning RATIO_TOCE IS NOT SAVE WITH NOMBLAY
116  ENDDO !END LOOP NPOIN
117  ENDDO !END LOOP NSAND_VIRTUAL
118 !
119  DO isand = 1,nsand
120  IF(nsand.NE.0)THEN
121  IF(susp_sand)THEN
122 ! COMPUTE EQUILIBRIUM CONCENTRATION FOR EACH LAYER!! IT IS NOT NECESSARY FOR EACH POINT
123 ! CAN BE OPTIMIZED IF NO LOOP ON NPOIN IN COMPUTE CAE
125  & dcla(num_isand_icla(isand)),npoin,charr,
126  & xmve,xmvs0(num_isand_icla(isand)),vce,grav,
127  & zero,zref,ac(num_isand_icla(isand)),
128  & cstaeq%ADR(num_isand_icla(isand))%P,
129  & qscl_c%ADR(num_isand_icla(isand))%P,
130  & icq,u2d,v2d,csratio,debug,
131  & ratio_toce%ADR(num_isand_icla(isand))%P)
132  ELSE
133  DO ipoin=1,npoin
134  cstaeq%ADR(num_isand_icla(isand))%P%R(ipoin)=0.d0
135  ENDDO
136  ENDIF
137 !
138  ENDIF
139 !!!!!!!!!!!!!warning CSTAEQ IS NOT SAVE WITH NOMBLAY
140 !!!!!!!!!!!!!for this time:
141  DO ipoin=1,npoin
142  cae_ilay(isand,ilayer,ipoin) =
143  & cstaeq%ADR(num_isand_icla(isand))%P%R(ipoin)
144  ENDDO
145  ENDDO !END LOOP NSAND
146  ENDDO !END LOOP NOMBLAY
147 !
148 !!!!!!!STRAT EROSION
149  DO ipoin = 1,npoin
150 !
151  DO imud = 1,nmud
152  qer_mud(imud) = 0.d0
153  ENDDO
154  DO isand = 1,nsand_virtual
155  qer_sand(isand) = 0.d0
156  time(isand) = dt
157  ENDDO
158 !
159  DO ilayer = 1,nomblay
160 !
161 ! COMPUTE PUR SAND FLUX
162 !
163  DO isand = 1,nsand
164  IF(nsand.NE.0)THEN
165 ! COMPUTE CAE MUST BE CALL HERE IF NO LOOP ON NPOIN
166  fluer_pur_sand(isand,ipoin) =
167  & cae_ilay(isand,ilayer,ipoin)*xwc(num_isand_icla(isand))
168  ELSE
169  cae_ilay(isand,ilayer,ipoin) = 0.d0
170  fluer_pur_sand(isand,ipoin) = 0.d0
171  ENDIF
172  ENDDO
173 ! COMPUTE PUR MUD FLUX
174  IF(nmud.GT.0)THEN
175  DO isand = 1,nsand_virtual
176  IF(taup%R(ipoin).GT.toce_mix(isand,ilayer,ipoin)) THEN
177  fluer_pur_mud(isand,ipoin) = partheniades(ilayer,ipoin)
178  & *(taup%R(ipoin)/toce_mix(isand,ilayer,ipoin) - 1.d0)
179  ELSE
180  fluer_pur_mud(isand,ipoin) = 0.d0
181  ENDIF
182  ENDDO
183  ELSE
184  DO isand = 1,nsand_virtual
185  fluer_pur_mud(isand,ipoin) = 0.d0
186  ENDDO
187  ENDIF
188 !
189 ! COMPUTE MIX FLUXES SAND-MUD
190  fluer_mix_tot = 0.d0
191 !
192  DO isand = 1,nsand_virtual
193 !
194  IF(taup%R(ipoin).GT.toce_mix(isand,ilayer,ipoin))THEN
195  IF(ratio_mud_sand(ilayer,ipoin).LE.0.3d0)THEN
196 ! MUD RATIO < 30%, (PURE SAND FLUX)
197  fluer_mix(isand,ipoin)= fluer_pur_sand(isand,ipoin)
198 ! MUD RATIO > 50%, (PURE MUD FLUX)
199  ELSEIF(ratio_mud_sand(ilayer,ipoin).GE.0.5d0)THEN
200  fluer_mix(isand,ipoin) = fluer_pur_mud(isand,ipoin)
201 ! MUD RATIO >30% AND <50%, (INTERPOLATION)
202  ELSE
203  fluer_mix(isand,ipoin) = (ratio_mud_sand(ilayer,ipoin)
204  & - 0.3d0)/(0.5d0-0.3d0)
205  & *(fluer_pur_mud(isand,ipoin)
206  & -fluer_pur_sand(isand,ipoin))
207  & + fluer_pur_sand(isand,ipoin)
208  ENDIF
209  ELSE
210  fluer_mix(isand,ipoin) = 0.d0
211  ENDIF
212 !
213  fluer_mix_tot = fluer_mix_tot
214  & + fluer_mix(isand,ipoin)*ratio_sand(isand,ilayer,ipoin)
215 !
216  ENDDO ! END LOOP NSAND_VIRTUAL
217 !
218  IF(fluer_mix_tot.LE.0.d0.AND.mass_mix_tot(ilayer,ipoin).GE.
219  & min_sed_mass_comp) GOTO 10
220 ! des valeurs infimes de masse ne peuvent pas faire du pavage, on laisse la possibilite d eroder la couche en dessous
221 ! EXIT LAYER LOOP
222  DO isand = 1,nsand_virtual
223 ! COMPUTATION OF ERODABLE QUANTITY
224  qe_moy(isand) = fluer_mix(isand,ipoin)*
225  & ratio_sand(isand,ilayer,ipoin)*time(isand)*mofac_bed
226 !
227  IF(qe_moy(isand).GT.0.d0) THEN
228 !
229  IF(qe_moy(isand).LT.
230  & (mass_mud_tot(ilayer,ipoin)
231  & *ratio_sand(isand,ilayer,ipoin)
232  & +mass_sand(isand,ilayer,ipoin))) THEN
233 !
234 ! IF ERODABLE QUANTITY < MASS IN LAYER , THE LAYER IS PARTIALLY ERODED
235 ! AND NO TIME TO ERODE THE NEXT LAYER
236 !
237 ! QUANTITY OF MUD ERODED RESPECT RATIO_MUD_SAND
238  check_positif_mud_tot1 =
239  & min(mass_mud_tot(ilayer,ipoin),
240  & qe_moy(isand)*ratio_mud_sand(ilayer,ipoin))
241  check_positif_mud_tot2 = 0.d0
242 !
243  DO imud = 1,nmud
244 ! QUANTITY OF EACH CLASS OF MUD ERODED RESPECT RATIO_MUD
245  check_positif_mud =
246  & min(mass_mud(imud,ilayer,ipoin),
247  & check_positif_mud_tot1
248  & *ratio_mud(imud,ilayer,ipoin))
249  mass_mud(imud,ilayer,ipoin) =
250  & mass_mud(imud,ilayer,ipoin) - check_positif_mud
251  qer_mud(imud) = qer_mud(imud) + check_positif_mud
252  check_positif_mud_tot2 =
253  & check_positif_mud_tot2 + check_positif_mud
254  ENDDO
255 !
256 ! QUANTITY OF SAND ERODED
257  check_positif_sand = min(mass_sand(isand,ilayer,ipoin),
258  & (1.d0-ratio_mud_sand(ilayer,ipoin))*qe_moy(isand))
259 
260  qer_sand(isand) = qer_sand(isand) + check_positif_sand
261  mass_sand(isand,ilayer,ipoin) =
262  & mass_sand(isand,ilayer,ipoin) - check_positif_sand
263  time(isand) = 0.d0
264  !
265  ELSE
266 !
267 ! IF ERODABLE QUANTITY > MASS IN LAYER , THE LAYER IS FULLY ERODED
268 ! AND THE TIME REMAINING FOR EROSION OF NEXT LAYER IS COMPUTED
269 !
270 ! QUANTITY OF MUD ERODED RESPECT RATIO_SAND
271  check_positif_mud_tot1 = min(mass_mud_tot(ilayer,ipoin),
272  & mass_mud_tot(ilayer,ipoin)
273  & *ratio_sand(isand,ilayer,ipoin))
274 !
275  check_positif_mud_tot2 = 0.d0
276  DO imud = 1,nmud
277 ! QUANTITY OF EACH CLASS OF MUD ERODED RESPECT RATIO_MUD
278  check_positif_mud =
279  & min(mass_mud(imud,ilayer,ipoin),
280  & check_positif_mud_tot1
281  & *ratio_mud(imud,ilayer,ipoin))
282 !
283  mass_mud(imud,ilayer,ipoin)=
284  & mass_mud(imud,ilayer,ipoin) - check_positif_mud
285 !
286  qer_mud(imud) = qer_mud(imud) + check_positif_mud
287 !
288  check_positif_mud_tot2 =
289  & check_positif_mud_tot2 + check_positif_mud
290  ENDDO
291 ! QUANTITY OF SAND ERODED
292  check_positif_sand =
293  & max(mass_sand(isand,ilayer,ipoin),0.d0)
294  mass_sand(isand,ilayer,ipoin) =
295  & mass_sand(isand,ilayer,ipoin) - check_positif_sand
296  qer_sand(isand) = qer_sand(isand) + check_positif_sand
297  ENDIF
298 ! TIME REMAINING FOR EROSION OF NEXT LAYER IS COMPUTED
299  time(isand) = time(isand) -
300  & ((check_positif_mud_tot2 + check_positif_sand)
301  & /(qe_moy(isand)/dt))
302 !
303  ENDIF
304 !
305  ENDDO ! END LOOP ISAND_VIRTUAL
306 !
307  DO isand = 1,nsand_virtual
308 ! IF REMAINING TIME IS NUL FOR ONE CLASS OF SAND, NO EROSION OF NEXT LAYER
309  IF(time(isand).LE.0.d0) GOTO 10
310 ! EXIT LAYER LOOP
311  ENDDO
312 !
313  ENDDO ! END LOOP NOMBLAY
314 10 CONTINUE
315 !
316  DO imud = 1,nmud
317  fluer%ADR(num_imud_icla(imud))%P%R(ipoin) =
318  & max(qer_mud(imud)/dt,0.d0)/mofac_bed
319 ! ADD MUD PART INCLUDED IN BEDLOAD
320  fluer%ADR(num_imud_icla(imud))%P%R(ipoin) =
321  & fluer%ADR(num_imud_icla(imud))%P%R(ipoin)+
322  & f_mudb%ADR(num_imud_icla(imud))%P%R(ipoin)
323  & /mofac_bed
324  ENDDO
325 !
326  DO isand = 1,nsand
327  fluer%ADR(num_isand_icla(isand))%P%R(ipoin) =
328  & max(qer_sand(isand)/dt,0.d0)/mofac_bed
329  ENDDO
330 !
331  ENDDO ! END LOOP NPOIN
332 !
333 ! MASS EVOLUTION DUE TO SUSPENSION: NEGATIVE SINCE EROSION
334  DO ipoin=1,npoin
335  DO icla=1,nsicla
336  evcl_ms%ADR(icla)%P%R(ipoin)=
337  & -fluer%ADR(icla)%P%R(ipoin)*dt
338  ENDDO
339  ENDDO
340 !
341  RETURN
342  END SUBROUTINE bed1_suspension_erode
double precision, dimension(:,:,:), allocatable cae_ilay
integer icq
Reference concentration formula.
double precision, dimension(:,:), allocatable fluer_mix
double precision, dimension(:,:), allocatable ratio_mud_sand
Ratio of mud to sand, for ilayer,ipoin.
type(bief_obj), target u2d
Components of depth-averaged velocity.
double precision, dimension(:,:), allocatable, target toce_mud
Critical erosion shear stress of the mud, for each bed layer, for each point. It is variable in space...
double precision, dimension(:,:,:), allocatable, target mass_mud
Surface mass of mud (kg/m2), for imud,ilayer,ipoin.
type(bief_obj), target v2d
integer, target nomblay
Number of bed load model layers = NUMSTRAT+1 to take the active layer into account.
double precision, dimension(nsiclm), target xwc
Settling velocities.
double precision, dimension(:), allocatable qer_sand
double precision, target dt
Time step It may be different from the one in TELEMAC because of the morphological factor...
double precision xmve
Water density (from steering file of T2D or T3D)
double precision, dimension(:,:), allocatable fluer_pur_mud
double precision, parameter min_sed_mass_comp
Minimum value to detect sediment mass.
type(bief_obj), target evcl_ms
Mass evolution for class (due to suspension)
logical hirano
Hirano model used.
type(bief_obj), target ratio_toce
Ratio between critical shear stress of pure sediment and mixed sediment in the same layer...
double precision, dimension(:,:,:), allocatable, target mass_sand
Surface mass of sand (kg/m2), for isand,ilayer,ipoin.
double precision vce
Water viscosity: it is defined here because the viscosity set in TELEMAC2D or TELEMAC3D may not b the...
subroutine suspension_compute_cae(TAUP, HN, DCLA, NPOIN, CHARR, XMVE, XMVS, VCE, GRAV, ZERO, ZREF, AC, CSTAEQ, QSC, ICQ, U2D, V2D, CSRATIO, DEBUG, RATIO_TOCE)
double precision, dimension(nsiclm) xmvs0
Sand density.
double precision, dimension(:,:), allocatable, target conc_mud
Mud concentration for each bed layer, at each point. It is variable in space for the active layer (la...
double precision, dimension(:,:), allocatable mass_mud_tot
Surface total mass of mud (kg/m2), for ilayer,ipoin.
double precision, dimension(:), allocatable qe_moy
integer nsand
Total number of sand.
double precision zero
Parameter used for clipping variables or testing values against zero.
integer, dimension(:), allocatable num_isand_icla
integer, dimension(:), allocatable num_imud_icla
Tables to switch from mud number to class number and from sand number to class number.
type(bief_obj), target hn
Water depth.
double precision, dimension(nsiclm), target dcla
Sediment diameter for each class It is only relevant for non-cohesive sediments. For the bedload...
double precision mofac_bed
Morphological factor on the bed: distorts the evolution of the morphodynamics with respect to the hyd...
subroutine bed1_suspension_erode
double precision, dimension(:,:), allocatable fluer_pur_sand
double precision, dimension(:,:), allocatable, target partheniades
Partheniades erosion coefficient, for each bed layer, for each point. It is variable in space for the...
integer, target nsicla
Number of sediment classes of bed material (less than NISCLM)
type(bief_obj), target zref
Reference elevation.
double precision grav
Gravity acceleration.
double precision, dimension(:,:,:), allocatable, target ratio_mud
Ratio of mud to all muds, for imud,ilayer,ipoin.
type(bief_obj), target csratio
Ratio between bottom concentration and average concentration.
integer nmud
Total number of muds.
logical, target susp_sand
Suspension for all sands (mud is assumed to be suspended)
double precision, dimension(:,:,:), allocatable, target ratio_sand
Ratio of sand to all sands, for isand,ilayer,ipoin.
double precision, dimension(:), allocatable qer_mud
double precision, dimension(nsiclm), target ac
Critical shields parameter.
double precision, dimension(:,:), allocatable toce_sand
Critical erosion shear stress of the sand, for each sand, for each point.
type(bief_obj), target taup
Shear stress modified by skin friction.
type(bief_obj), target qscl_c
Bedload transport rate for a sediment class [kg*(m-1*s-1)].
integer debug
Debugger.
type(bief_obj), target cstaeq
Sediment equilibrium concentration.
logical, target charr
Include bedload in the simulation.
type(bief_obj), target fluer
Erosion flux.
double precision, dimension(:), allocatable time
double precision, dimension(:,:,:), allocatable toce_mix
type(bief_obj), target f_mudb
Flux of mud in bedload added in suspension.
double precision, dimension(:,:), allocatable mass_mix_tot
Surface total mass of sediments (kg/m2), for ilayer,ipoin.
integer, pointer npoin
Number of 2d points in the mesh.
double precision, dimension(:), allocatable conc_mud_activ_layer
Concentation of mud in active layer (array defined for temporary work in some subroutines) ...
Definition: bief.f:3