The TELEMAC-MASCARET system  trunk
mass_balance.f
Go to the documentation of this file.
1 ! ***********************
2  SUBROUTINE mass_balance
3 ! ***********************
4 !
5  &(dt,nptfr,entet,nsicla,numliq,nfrliq,flbcla,
6  & lt,nit,npoin,volu2d,charr,susp,evcl_mb,evcl_ms,masstot,
7  & mass0tot)
8 !
9 !***********************************************************************
10 ! GAIA
11 !***********************************************************************
12 !
18 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
37 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
38 !
39  USE bief
40  USE interface_gaia, ex_mass_balance => mass_balance
42  & fluer,
43  & fludp, summcumucl, sumrmascl,
44  & sum_erosion,sum_deposition,
45  & massnestor,summassnestor,
46  & mass0act,nestor,bedload_b_flux,
47  & sumbedload_b_flux,sumbedload_b,
48  & mcumucla
49 !
51  USE interface_parallel, ONLY: p_sum
52  IMPLICIT NONE
53 !
54 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
55 !
56  INTEGER, INTENT(IN) :: NPTFR,NFRLIQ,NSICLA,LT,NIT
57  INTEGER, INTENT(IN) :: NPOIN,NUMLIQ(nptfr)
58  DOUBLE PRECISION, INTENT(IN) :: DT
59  LOGICAL, INTENT(IN) :: ENTET,SUSP,CHARR
60  DOUBLE PRECISION, INTENT(IN) :: MASSTOT(nsicla)
61  DOUBLE PRECISION, INTENT(IN) :: MASS0TOT(nsicla)
62 !
63 !-----------------------------------------------------------------------
64 !
65 ! VECTOR STRUCTURES
66 !
67  TYPE(bief_obj), INTENT(IN) :: VOLU2D
68  TYPE(bief_obj), INTENT(IN) :: EVCL_MB,EVCL_MS
69  TYPE(bief_obj), INTENT(INOUT) :: FLBCLA
70 !
71 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
72 !
73  INTEGER I,IFRLQ,IPTFR,ICLA,K
74  DOUBLE PRECISION LOST,RELATI,DENOM
75  DOUBLE PRECISION FLUXTCLA
76  DOUBLE PRECISION RMASCLA(nsiclm)
77  DOUBLE PRECISION EROSION_FLUX(nsiclm),DEPOSITION_FLUX(nsiclm)
78  DOUBLE PRECISION SUMMASSTOT,SUMMASS0TOT
79  DOUBLE PRECISION SUMRMASCLA,SUMMCUMUCLA
80  DOUBLE PRECISION SUMEROSTOT,SUMDEPOTOT
81  DOUBLE PRECISION SUMLOST
82  DOUBLE PRECISION SUMNESTOR,SUMMASS0ACT
83 !
84 !-----------------------------------------------------------------------
85 !
86 ! COMPUTES EROSION AND DEPOSITION FLUX PER CLASS EVOLUTION
87 ! FOR THIS TIME STEP
88 ! NOTE: EROSION AND DEPOSITION FLUX MUST BE 0 IF THERE IS
89 ! NO SUSPENSION
90 !
91  IF(nsicla.GT.0) THEN
92  DO icla=1,nsicla
93  erosion_flux(icla) = 0.d0
94  deposition_flux(icla) = 0.d0
95  ENDDO
96  ENDIF
97  IF(nsicla.GT.0) THEN
98  IF (susp) THEN
99  DO icla=1,nsicla
100  DO i=1,npoin
101  erosion_flux(icla) = erosion_flux(icla)
102  & + fluer%ADR(icla)%P%R(i)*volu2d%R(i)
103  ENDDO
104  IF(ncsize.GT.1) erosion_flux(icla) =
105  & p_sum(erosion_flux(icla))
106  DO i=1,npoin
107  deposition_flux(icla) = deposition_flux(icla)
108  & + fludp%ADR(icla)%P%R(i)*volu2d%R(i)
109  ENDDO
110  IF(ncsize.GT.1) deposition_flux(icla) =
111  & p_sum(deposition_flux(icla))
112  ENDDO
113  ENDIF
114  ENDIF
115 !
116 !=======================================================================
117 !
118 ! PREPARES TERMS FOR BALANCE IN EXTENDED GRANULOMETRY
119 !
120  IF(nsicla.GT.0) THEN
121 !
122  DO icla=1,nsicla
123 !
124  k = num_icla_isand(icla)
125 ! COMPUTES THE EVOLUTION PER CLASS
126 !
127  rmascla(icla) = 0.d0
128  IF(susp.AND.charr) THEN
129  DO i=1,npoin
130  rmascla(icla) = rmascla(icla)+volu2d%R(i)*
131  & (evcl_mb%ADR(icla)%P%R(i) + evcl_ms%ADR(icla)%P%R(i))
132  ENDDO
133  ELSEIF(susp) THEN
134  DO i=1,npoin
135  rmascla(icla) = rmascla(icla)
136  & + evcl_ms%ADR(icla)%P%R(i)*volu2d%R(i)
137  ENDDO
138  ELSEIF(charr) THEN
139  DO i=1,npoin
140  rmascla(icla) = rmascla(icla)
141  & + evcl_mb%ADR(icla)%P%R(i)*volu2d%R(i)
142  ENDDO
143  ENDIF
144  IF(ncsize.GT.1) rmascla(icla) = p_sum(rmascla(icla))
145 !
146 ! COMPUTES THE FREE FLUXES BY CLASS, FOR EVERY BOUNDARY:
147 ! BEDLOAD_B_FLUX
148 ! FLUXTCLA CONTAINS THE SUM OF BEDLOAD_B_FLUX OVER BOUNDARIES
149 !
150  fluxtcla = 0.d0
151  IF(nfrliq.GT.0.AND.charr) THEN
152  DO ifrlq=1,nfrliq
153  bedload_b_flux(ifrlq,icla) = 0.d0
154  ENDDO
155  IF(nptfr.GT.0) THEN
156  DO iptfr=1,nptfr
157  ifrlq = numliq(iptfr)
158  IF(ifrlq.GT.0) THEN
159  bedload_b_flux(ifrlq,icla) =
160  & bedload_b_flux(ifrlq,icla)+flbcla%ADR(icla)%P%R(iptfr)
161  fluxtcla = fluxtcla + flbcla%ADR(icla)%P%R(iptfr)
162  ENDIF
163  ENDDO
164  ENDIF
165  IF(ncsize.GT.1) THEN
166  fluxtcla = p_sum(fluxtcla)
167  DO i=1,nfrliq
168  bedload_b_flux(i,icla) = p_sum(bedload_b_flux(i,icla))
169  ENDDO
170  ENDIF
171  ENDIF
172 !
173 ! COMPUTES THE CUMULATED FLUXES PER CLASS
174 !
175  mcumucla(icla) = -fluxtcla
176  summcumucl(icla) = summcumucl(icla) + mcumucla(icla)*dt
177  DO ifrlq=1,nfrliq
178  sumbedload_b(ifrlq,icla) = sumbedload_b(ifrlq,icla) +
179  & bedload_b_flux(ifrlq,icla)*dt
180  ENDDO
181  sumrmascl(icla) = sumrmascl(icla) + rmascla(icla)
182  sum_deposition(icla) = sum_deposition(icla)
183  & + deposition_flux(icla)*dt
184  sum_erosion(icla) = sum_erosion(icla)
185  & + erosion_flux(icla)*dt
186  summassnestor(icla) = summassnestor(icla) + massnestor(icla)
187  ENDDO
188 !
189  ENDIF
190 !
191 !=======================================================================
192 !
193 ! FIXME: IF(NESTOR)...
194 !
195 ! WRITES OUT THE BALANCE
196 !
197  IF(entet) THEN
198 !
199  IF(nsicla.GT.0) THEN
200 !
201 ! INITIALIZE FOR BALANCE OVER ALL CLASSES
202 !
203  summass0tot = 0.d0
204  summasstot = 0.d0
205  summass0act = 0.d0
206  sumrmascla = 0.d0
207  summcumucla = 0.d0
208  sumerostot = 0.d0
209  sumdepotot = 0.d0
210  sumlost = 0.d0
211  sumnestor = 0.d0
212  IF(nfrliq.GT.0) THEN
213  DO ifrlq=1,nfrliq
214  sumbedload_b_flux(ifrlq) = 0.d0
215  ENDDO
216  ENDIF
217 !
218 ! BALANCE FOR EACH CLASS
219 !
220  DO i=1,nsicla
221 ! EROSION AND DEPOSITION FLUX MUST BE 0 IF THERE IS
222 ! NO SUSPENSION
223  lost = rmascla(i)+(-mcumucla(i) + erosion_flux(i)
224  & - deposition_flux(i))*dt - massnestor(i)
225  WRITE(lu,*)
226  WRITE(lu,2000)
227  WRITE(lu,2001) i
228  WRITE(lu,2010) rmascla(i)
229  WRITE(lu,3031) mcumucla(i)
230  IF(nfrliq.GT.0.AND.charr) THEN
231  DO ifrlq=1,nfrliq
232  WRITE(lu,3033) ifrlq,-bedload_b_flux(ifrlq,i)
233 ! SUMS FLUXES OVER ALL CLASSES
234  sumbedload_b_flux(ifrlq) = sumbedload_b_flux(ifrlq)
235  & + bedload_b_flux(ifrlq,i)
236  ENDDO
237  ENDIF
238  WRITE(lu,2112) erosion_flux(i)
239  WRITE(lu,2113) deposition_flux(i)
240  WRITE(lu,2910) masstot(i)
241  IF(nestor) WRITE(lu,2022) massnestor(i)
242  WRITE(lu,2033) lost
243  IF(mass0act(i).GT.1.d-8) THEN
244  relati = lost/mass0act(i)
245  WRITE(lu,2130) relati
246  ENDIF
247 ! FIXME: RELATIVE ERROR SEEMS STRANGE.. KEEP IT FOR THE MOMENT
248 ! RK: IS RELATIVE TO THE INITIAL MASS IN ACTIVE LAYER
249  denom = max(rmascla(i),(mcumucla(i) + erosion_flux(i)
250  & -deposition_flux(i))*dt)
251 !
252 ! UPDATES VARIABLES FOR BALANCE OVER ALL CLASSES
253 !
254  summasstot = summasstot + masstot(i)
255  summass0tot = summass0tot + mass0tot(i)
256  summass0act = summass0act + mass0act(i)
257  sumrmascla = sumrmascla + rmascla(i)
258  summcumucla = summcumucla + mcumucla(i)
259  sumerostot = sumerostot + erosion_flux(i)
260  sumdepotot = sumdepotot + deposition_flux(i)
261  sumnestor = sumnestor + massnestor(i)
262  sumlost = sumlost + lost
263 !
264  IF(denom.GT.1.d-8) THEN
265  lost = lost / denom
266  WRITE(lu,2240) lost
267  ENDIF
268  ENDDO
269 !
270 ! BALANCE OVER ALL CLASSES
271 !
272  IF(nsicla.GT.1) THEN
273  WRITE(lu,*)
274  WRITE(lu,2002)
275  WRITE(lu,3010) sumrmascla
276  WRITE(lu,3031) summcumucla
277  IF(nfrliq.GT.0.AND.charr) THEN
278  DO ifrlq=1,nfrliq
279  WRITE(lu,3033) ifrlq,-sumbedload_b_flux(ifrlq)
280  ENDDO
281  ENDIF
282  WRITE(lu,2112) sumerostot
283  WRITE(lu,2113) sumdepotot
284  WRITE(lu,2910) summasstot
285  IF(nestor) WRITE(lu,2022) sumnestor
286  WRITE(lu,2033) sumlost
287  IF(summass0act.GT.1.d-8) THEN
288  relati = sumlost / summass0act
289  WRITE(lu,2130) relati
290  ENDIF
291  ENDIF
292 !
293 ! FINAL GLOBAL BALANCE
294 !
295  IF(lt.EQ.nit) THEN
296 !
297 ! INITIALIZE FOR BALANCE OVER ALL CLASSES
298 !
299  summass0tot = 0.d0
300  summasstot = 0.d0
301  summass0act = 0.d0
302  sumrmascla = 0.d0
303  summcumucla = 0.d0
304  sumerostot = 0.d0
305  sumdepotot = 0.d0
306  sumlost = 0.d0
307  sumnestor = 0.d0
308  IF(nfrliq.GT.0) THEN
309  DO ifrlq=1,nfrliq
310  sumbedload_b_flux(ifrlq) = 0.d0
311  ENDDO
312  ENDIF
313 !
314 ! BALANCE FOR EACH CLASS
315 !
316  WRITE(lu,*)
317  WRITE(lu,2150)
318  DO i=1,nsicla
319  lost = sumrmascl(i) - summcumucl(i) + sum_erosion(i)
320  & - sum_deposition(i) - summassnestor(i)
321  WRITE(lu,*)
322  WRITE(lu,2000)
323  WRITE(lu,2001) i
324  WRITE(lu,3010) sumrmascl(i)
325  WRITE(lu,3032) summcumucl(i)
326  IF(nfrliq.GT.0.AND.charr) THEN
327  DO ifrlq=1,nfrliq
328  WRITE(lu,3034) ifrlq,-sumbedload_b(ifrlq,i)
329  sumbedload_b_flux(ifrlq) = sumbedload_b_flux(ifrlq)
330  & + sumbedload_b(ifrlq,i)
331  ENDDO
332  ENDIF
333  WRITE(lu,2114) sum_erosion(i)
334  WRITE(lu,2115) sum_deposition(i)
335  WRITE(lu,2028) mass0tot(i)
336  WRITE(lu,2029) mass0act(i)
337  WRITE(lu,2910) masstot(i)
338  IF(nestor) WRITE(lu,2022) summassnestor(i)
339  WRITE(lu,2034) lost
340  IF(mass0act(i).GT.1.d-8) THEN
341  relati = lost / mass0act(i)
342  WRITE(lu,2130) relati
343  ENDIF
344 !
345  summass0tot = summass0tot + mass0tot(i)
346  summasstot = summasstot + masstot(i)
347  summass0act = summass0act + mass0act(i)
348  sumrmascla = sumrmascla + sumrmascl(i)
349  summcumucla = summcumucla + summcumucl(i)
350  sumerostot = sumerostot + sum_erosion(i)
351  sumdepotot = sumdepotot + sum_deposition(i)
352  sumnestor = sumnestor + summassnestor(i)
353  sumlost = sumlost + lost
354 !
355  IF(denom.GT.1.d-8) THEN
356  lost = lost / denom
357  WRITE(lu,2160) lost
358  ENDIF
359  lost = summasstot - summass0tot - sumnestor
360  & - summcumucla + sumerostot - sumdepotot
361  WRITE(lu,2031) lost
362  IF(summass0tot.GT.1.d-8) THEN
363  relati = lost / summass0tot
364  WRITE(lu,2131) relati
365  ENDIF
366 
367 
368  ENDDO
369 !
370 ! BALANCE OVER ALL CLASSES
371 !
372  IF(nsicla.GT.1) THEN
373  WRITE(lu,*)
374  WRITE(lu,2002)
375  WRITE(lu,3010) sumrmascla
376  WRITE(lu,3032) summcumucla
377  IF(nfrliq.GT.0.AND.charr) THEN
378  DO ifrlq=1,nfrliq
379  WRITE(lu,3034) ifrlq,-sumbedload_b_flux(ifrlq)
380  ENDDO
381  ENDIF
382  WRITE(lu,2114) sumerostot
383  WRITE(lu,2115) sumdepotot
384  WRITE(lu,2028) summass0tot
385  WRITE(lu,2029) summass0act
386  WRITE(lu,2910) summasstot
387  IF(nestor) WRITE(lu,2022) sumnestor
388  WRITE(lu,2034) sumlost
389  IF(summass0act.GT.1.d-8) THEN
390  relati = sumlost / summass0act
391  WRITE(lu,2130) relati
392  ENDIF
393  lost = summasstot - summass0tot - sumnestor
394  & - summcumucla + sumerostot - sumdepotot
395  WRITE(lu,2031) lost
396  IF(summass0tot.GT.1.d-8) THEN
397  relati = lost / summass0tot
398  WRITE(lu,2131) relati
399  ENDIF
400  ENDIF
401 !
402  ENDIF
403  ENDIF ! ENDIF NSICLA
404 !
405  ENDIF ! ENDIF ENTET
406 !
407 2000 FORMAT(20x,'GAIA MASS-BALANCE OF SEDIMENTS PER CLASS: ')
408 2002 FORMAT(20x,'GAIA MASS-BALANCE OF SEDIMENTS OVER ALL CLASSES: ')
409 2001 FORMAT(5x,'SEDIMENT CLASS NUMBER = ',1i8)
410 2010 FORMAT(5x,'TOTAL BED EVOLUTIONS = ',g16.7,
411  & ' ( KG )')
412 2028 FORMAT(5x,'INITIAL MASS = ',g16.7,
413  & ' ( KG )')
414 2029 FORMAT(5x,'INITIAL MASS ACTIVE LAYER = ',g16.7,
415  & ' ( KG )')
416 2033 FORMAT(5x,'LOST MASS = ',g16.7,
417  & ' ( KG )')
418 2034 FORMAT(5x,'CUMULATED LOST MASS = ',g16.7,
419  & ' ( KG )')
420 2112 FORMAT(5x,'EROSION FLUX = ',g16.7,
421  & ' ( KG/S )')
422 2113 FORMAT(5x,'DEPOSITION FLUX = ',g16.7,
423  & ' ( KG/S )')
424 2114 FORMAT(5x,'CUMULATED EROSION = ',g16.7,
425  & ' ( KG )')
426 2115 FORMAT(5x,'CUMULATED DEPOSITION = ',g16.7,
427  & ' ( KG )')
428 2150 FORMAT(20x,'FINAL MASS-BALANCE OF SEDIMENTS: ')
429 2130 FORMAT(5x,'RELATIVE ERROR TO INITIAL ACT LAYER MASS = ',g16.7)
430 2131 FORMAT(5x,'RELATIVE ERROR TO TOTAL INITIAL MASS = ',g16.7)
431 2160 FORMAT(5x,'CUMULATED RELATIVE ERROR ON MASS = ',g16.7)
432 2240 FORMAT(5x,'RELATIVE ERROR 0N MASS = ',g16.7)
433 2022 FORMAT(5x,'NESTOR MASS = ',g16.7,
434  & ' ( KG )')
435 2910 FORMAT(5x,'TOTAL MASS = ',g16.7,
436  & ' ( KG )')
437 3010 FORMAT(5x,'CUMULATED BED EVOLUTIONS = ',g16.7,
438  & ' ( KG )')
439 3031 FORMAT(5x,'BOUNDARIES BEDLOAD FLUX = '
440  & ,g16.7,' ( KG/S >0 = ENTERING )')
441 3032 FORMAT(5x,'CUMULATED BOUNDARIES BEDLOAD MASS = ',g16.7,
442  & ' ( KG )')
443 3033 FORMAT(5x,'BEDLOAD FLUX BOUNDARY ',i4,' = ', g16.7
444  & ,' ( KG/S >0 = ENTERING )')
445 3034 FORMAT(5x,'CUMULATED BEDLOAD BOUNDARY ',i4,' = ', g16.7
446  & ,' ( KG )')
447 2031 FORMAT(5x,'CUMULATED LOST MASS (INI-FINAL+FLUXES) = ',
448  & g16.7,' ( KG )')
449 
450 !
451 !-----------------------------------------------------------------------
452 !
453  RETURN
454  END
subroutine mass_balance(DT, NPTFR, ENTET, NSICLA, NUMLIQ, NFRLIQ, FLBCLA, LT, NIT, NPOIN, VOLU2D, CHARR, SUSP, EVCL_MB, EVCL_MS, MASSTOT, MASS0TOT)
Definition: mass_balance.f:9
integer, parameter nsiclm
Maximum number of sediment classes.
integer, dimension(:), allocatable num_icla_isand
Definition: bief.f:3