The TELEMAC-MASCARET system  trunk
bed1_init_sediment_gaia.f
Go to the documentation of this file.
1 ! **********************************
2  SUBROUTINE bed1_init_sediment_gaia
3 ! **********************************
4 !
5  &(nsicla,elay,zf,zr,npoin,xmvs0,es,nomblay,
6  & debu,volu2d,numstrat,maxvar)
7 !
8 !***********************************************************************
9 ! GAIA
10 !***********************************************************************
11 !
13 !
14 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
28 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
29 !
30  USE bief
31  USE interface_gaia, ex_bed1_init_sediment =>
34  & mass_mud_tot,mass_sand_tot,
35  & ratio_sand,ratio_mud,ratio_mud_sand,nsand,nmud,
36  & masstot,num_imud_icla,num_isand_icla,
37  & elay0,debu_mass,mass0tot,ava0,xkv0,
38  & conc_mud0,conc_mud,toce_mud0,toce_mud,
39  & trans_mass0,trans_mass,partheniades0,partheniades,hirano,
40  & flux_neg_mud_activ_layer,flux_pos_mud_activ_layer,sed_thick,
41  & flux_neg_sand_activ_layer,flux_pos_sand_activ_layer,bed_model,
42  & conc_mud_found,toce_mud_found,partheniades_found,mtrans_found
43 
45  USE interface_parallel, ONLY: p_sum
46  IMPLICIT NONE
47 !
48 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
49 !
50  INTEGER, INTENT(IN) :: NSICLA,NPOIN,NOMBLAY,MAXVAR
51  INTEGER, INTENT(IN) :: NUMSTRAT
52  TYPE(bief_obj), INTENT(INOUT) :: ELAY,ZF,ZR
53  TYPE(bief_obj), INTENT(INOUT) :: VOLU2D
54  DOUBLE PRECISION, INTENT(IN) :: XMVS0(nsicla)
55  LOGICAL, INTENT(IN) :: DEBU
56  DOUBLE PRECISION, INTENT(INOUT) :: ES(npoin,nomblay)
57 !
58 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
59 !
60  INTEGER :: K,ICLA,IPOIN,IERR
61  INTEGER :: ILAYER,IMUD,ISAND,ISTRAT
62  DOUBLE PRECISION :: CHECK_RS,CHECK_RM
63  DOUBLE PRECISION :: TERM,DISCR,TOT1,TOT2
64  DOUBLE PRECISION :: MASS_TOT,TOTSAND,TOTMUD
65  DOUBLE PRECISION :: XMVS_LAY
66  DOUBLE PRECISION, ALLOCATABLE :: ESTRATUM(:,:)
67  DOUBLE PRECISION, ALLOCATABLE :: RATIO_INIT(:,:,:)
68 !
69 !======================================================================!
70 !
71 
72  IF(debu) THEN ! This is checked ONLY
73  IF(nsicla.GT.1) hirano=.true. ! for pure a bedload case
74  ENDIF
75 !
76  IF(.NOT.debu) THEN
77 ! CASE NOT COMPUTATION CONTINUED
78 !
79 ! ALLOCATE ESTRATUM AND RATIO_INIT
80  ALLOCATE(estratum(numstrat,npoin), stat=ierr)
81  CALL check_allocate(ierr, 'ESTRATUM')
82  ALLOCATE(ratio_init(nsicla,numstrat,npoin), stat=ierr)
83  CALL check_allocate(ierr, 'RATIO_SAND')
84 !
85 ! BED COMPOSITION:
86 !
87 ! DEFAULT CASE : NO STRATIFICATION = ONLY ONE STRATUM 100 M DEEP
88 !
89 ! USER CAN CHANGE THE THICKNESS OF SEDIMENT HERE
90 ! (REPLACES SUBROUTINE NOEROD)
91 ! GRADED SEDIMENT: USER CAN DEFINE AN INITIAL STRATIFICATION
92 ! DEFINED BY LAYER THICKNESS AND COMPOSITION FOR EACH STRATUM
93 !
94 ! POROSITY IS DEFINED (KEYWORD) FOR STRATUMS.
95 ! THE VALUE FOR THE FIRST STRATUM IS COPIED IN
96 ! THE FIRST TWO NUMERICAL LAYERS (ACTIVE LAYER + FIRST STRATUM)
97  DO ipoin=1,npoin
98  DO istrat=1,numstrat
99 ! DEFAULT CASE: ALL STRATUMS HAVE SAME THICKNESS.
100 ! THIS CAN BE CHANGED BY USER IN USER_BED_INIT
101  estratum(istrat,ipoin) = sed_thick(istrat)
102 ! DEFAULT CASE: ALL STRATUMS HAVE SAME COMPOSITION.
103 ! THIS CAN BE CHANGED BY USER IN USER_BED_INIT
104  DO icla=1,nsicla
105  ratio_init(icla,istrat,ipoin) = ava0(icla)
106  ENDDO
107  ENDDO
108  ENDDO
109 !
110 ! USER CAN DO ADVANCED SETUP HERE
111  CALL user_bed_init(numstrat,npoin,nsicla,estratum,ratio_init)
112 !
113 ! INITIALISATION OF NUMERICAL THICKNESS (ES):
114  DO ipoin=1,npoin
115  DO ilayer=1,nomblay
116  es(ipoin,ilayer)=0.d0
117  ENDDO
118  ENDDO
119 !
120 ! INITIALISATION OF FLUX FOR ACTIVE LAYER IN CASE OF CONSOLIDATION
121 ! IF NOT HIRANO, FLUX HAVE TO STAY AT 0 BECAUSE THERE AS USE IN THE GENERAL CASE.
122  DO ipoin = 1,npoin
123  DO imud = 1,nmud
124  flux_neg_mud_activ_layer(imud,ipoin)= 0.d0
125  DO ilayer=1,nomblay
126  flux_pos_mud_activ_layer(imud,ilayer,ipoin)=0.d0
127  ENDDO
128  ENDDO
129  ENDDO
130  DO ipoin = 1,npoin
131  DO isand = 1,nsand
132  flux_neg_sand_activ_layer(isand,ipoin)=0.d0
133  DO ilayer=1,nomblay
134  flux_pos_sand_activ_layer(isand,ilayer,ipoin)=0.d0
135  ENDDO
136  ENDDO
137  ENDDO
138 !
139  IF(nsicla.EQ.1) THEN
140 !
141  hirano=.false.
142 !
143  IF(nomblay.GT.1.AND.nmud.EQ.0) THEN
144  WRITE(lu,*) 'ERROR: YOU HAVE ONLY 1 NON-COHESIVE CLASS '
145  WRITE(lu,*) 'AND SEVERAL LAYERS. THIS IS NOT ALLOWED!'
146  CALL plante(1)
147  stop
148  ENDIF
149 !
150 !
151  DO ipoin = 1,npoin
152  DO ilayer=1,nomblay
153  es(ipoin,ilayer) = estratum(ilayer,ipoin)
154  conc_mud(ilayer,ipoin) = conc_mud0(ilayer)
155  trans_mass(ilayer,ipoin) = trans_mass0(ilayer)
156  toce_mud(ilayer,ipoin) = toce_mud0(ilayer)
157  partheniades(ilayer,ipoin)=partheniades0(ilayer)
158  ENDDO
159  ENDDO
160 !
161  ELSE
162 !
163 ! MULTICLASS CASE
164 !
165  hirano=.true.
166 !
167  DO ipoin=1,npoin
168  elay%R(ipoin) = elay0
169 ! ELAY WILL BE CALCULATED IN BED1_UPDATE_ACTIVELAYER_HIRANO (AGAIN)
170 !
171 ! THE ACTIVE LAYER CORRESPONDS TO THE FIRST LAYER
172 ! ITS THICKNESS IS AUTOMATICALLY COMPUTED IN THE FIRST CALL TO
173 ! BED_UPDATE_HIRANO AND VALUE FOR MUD ARE UPDATE FOR THIS LAYER
174 ! BUT TRANS_MAS MUST BE EQUAL TO 0 FOR THIS LAYER
175 !
176  es(ipoin,1) = 0.d0 !!ES OF LAYER 1 COMPUTED AFTER IN HIRANO
177  conc_mud(1,ipoin) = conc_mud0(1) !!CONC_MUD OF LAYER 1 IS COMPUTED AFTER IN HIRANO
178  trans_mass(1,ipoin) = trans_mass0(1)!!TRANS_MASS OF LAYER 1 IS COMPUTED AFTER IN HIRANO
179  toce_mud(1,ipoin)= toce_mud0(1) !!TOCE_MUD OF LAYER 1 COMPUTED AFTER IN HIRANO
180  partheniades(1,ipoin)=partheniades0(1)!!TO CHANGE IF NECESSARY
181 
182  DO ilayer = 2,nomblay
183  istrat=ilayer-1
184  es(ipoin,ilayer) = estratum(istrat,ipoin)
185  conc_mud(ilayer,ipoin) = conc_mud0(istrat)
186  trans_mass(ilayer,ipoin) = trans_mass0(istrat)
187  toce_mud(ilayer,ipoin) = toce_mud0(istrat)
188  partheniades(ilayer,ipoin)=partheniades0(istrat)
189  ENDDO
190  ENDDO
191 !
192  ENDIF
193 ! CHECK PARAMETER VALUES FOR CONSOLIDATION
194  IF(bed_model.EQ.2)THEN
195  DO ipoin=1,npoin
196  DO ilayer = 2,nomblay
197  IF(conc_mud(ilayer,ipoin).LT.
198  & conc_mud(ilayer-1,ipoin)) THEN
199  WRITE(lu,*)'MUD CONCENTRATION OF ILAYER MUST BE',
200  & ' < TO MUD CONCNTRATION OF ILAYER +1 '
201  CALL plante(1)
202  stop
203  ENDIF
204  ENDDO
205  IF(trans_mass(nomblay,ipoin).LT.0.d0.OR.
206  & trans_mass(nomblay,ipoin).GT.1.d-8) THEN
207  WRITE(lu,*)'MASS TRANSFERT FOR LAST LAYER ',
208  & 'OF CONSOLIDATIONMUST BE EQUAL TO 0'
209  CALL plante(1)
210  stop
211  ENDIF
212  ENDDO
213  ENDIF
214 !
215 ! HERE THE NON ERODABLE BED IS FIXED
216 !
217  DO ipoin=1,npoin
218  zr%R(ipoin)=zf%R(ipoin)
219  DO ilayer=1,nomblay
220  zr%R(ipoin)=zr%R(ipoin)-es(ipoin,ilayer)
221  ENDDO
222  ENDDO
223 !
224 ! RATIO_SAND, RATIO_MUD AND RATIO_MUD_SAND ARE COMPUTED FROM RATIO_INIT(NCLA,NUMSTRAT,NPOIN)
225 !
226 ! 1: RATIO_SAND
227 !
228  ratio_sand = 0.d0
229  IF(nsand.GT.0) THEN
230  DO ipoin=1,npoin
231  DO istrat=1,numstrat
232  tot1=0.d0
233  tot2=0.d0
234  DO isand=1,nsand
235  tot1=tot1+ratio_init(num_isand_icla(isand),istrat,ipoin)
236  ENDDO
237  DO isand=1,nsand
238  IF(nsicla.GT.1)THEN
239  ilayer = istrat+1
240 ! ACTIVE LAYER IS FILLED IN BED_UPDATE_HIRANO BUT ITS
241 ! RATIO IS INIITIALISED HERE
242  IF(isand.LT.nsand)THEN
243  ratio_sand(isand,1,ipoin)=0.d0
244  ELSE
245  ratio_sand(nsand,1,ipoin)=1.d0
246  ENDIF
247  ELSE
248  ilayer=istrat
249  ENDIF
250  IF(isand.LT.nsand)THEN
251  IF(tot1.GT.0.d0)THEN
252  ratio_sand(isand,ilayer,ipoin)=
253  & ratio_init(num_isand_icla(isand),istrat,ipoin)/tot1
254  ELSE
255  ratio_sand(nsand,ilayer,ipoin)=0.d0
256  ENDIF
257  tot2=tot2+ratio_sand(isand,ilayer,ipoin)
258  ELSE
259  ratio_sand(nsand,ilayer,ipoin)=1.d0-tot2
260  ENDIF
261  ENDDO
262  ENDDO
263  ENDDO
264  ENDIF
265 !
266 ! 2: RATIO_MUD
267 !
268  ratio_mud = 0.d0
269  IF(nmud.GT.0) THEN
270  DO ipoin=1,npoin
271  DO istrat=1,numstrat
272  tot1=0.d0
273  tot2=0.d0
274  DO imud=1,nmud
275  tot1= tot1 + ratio_init(num_imud_icla(imud),istrat,ipoin)
276  ENDDO
277  DO imud=1,nmud
278  IF(nsicla.GT.1)THEN
279  ilayer = istrat+1
280 ! ACTIVE LAYER IS FILLED IN BED_UPDATE_HIRANO BUT ITS
281 ! RATIO IS INIITIALISED HERE
282  IF(imud.LT.nmud)THEN
283  ratio_mud(imud,1,ipoin)=0.d0
284  ELSE
285  ratio_mud(nmud,1,ipoin)=1.d0
286  ENDIF
287  ELSE
288  ilayer=istrat
289  ENDIF
290  IF(imud.LT.nmud)THEN
291  IF(tot1.GT.0.d0)THEN
292  ratio_mud(imud,ilayer,ipoin)=
293  & ratio_init(num_imud_icla(imud),istrat,ipoin)/tot1
294  ELSE
295  ratio_mud(nmud,ilayer,ipoin)=0.d0
296  ENDIF
297  tot2=tot2+ratio_mud(imud,ilayer,ipoin)
298  ELSE
299  ratio_mud(nmud,ilayer,ipoin)=1.d0-tot2
300  ENDIF
301  ENDDO
302  ENDDO
303  ENDDO
304  ENDIF
305 !
306 ! 3: RATIO_MUD_SAND
307 !
308  ratio_mud_sand = 0.d0
309  DO ipoin=1,npoin
310  DO istrat=1,numstrat
311  totmud=0.d0
312  totsand=0.d0
313  DO imud=1,nmud
314  totmud= totmud +
315  & ratio_init(num_imud_icla(imud),istrat,ipoin)
316  ENDDO
317  DO isand=1,nsand
318  totsand= totsand +
319  & ratio_init(num_isand_icla(isand),istrat,ipoin)
320  ENDDO
321  IF(nsicla.GT.1)THEN
322  ilayer = istrat+1
323 ! ACTIVE LAYER IS FILLED IN THE FIRST CALL TO BED_UPDATE_HIRANO
324  ELSE
325  ilayer=istrat
326  ENDIF
327  IF(totmud+totsand.GT.0.d0)THEN
328  ratio_mud_sand(ilayer,ipoin)=totmud/(totmud+totsand)
329  ELSE
330  ratio_mud_sand(ilayer,ipoin)=1.d0
331  ENDIF
332  ENDDO
333  ENDDO
334 !
335  ELSE
336 !
337 ! DEBU.EQ.TRUE
338 ! IN THIS PART WE COMPLETE MISSING DATA FROM PREV SED FILE
339 ! USING DATA READ IN THE STEERING FILE
340  IF(nsicla.EQ.1) THEN
341  DO ipoin=1,npoin
342  DO ilayer=1,nomblay
343  DO isand=1,nsand
344  ratio_sand(isand,ilayer,ipoin) =
345  & ava0(num_isand_icla(isand))
346  ENDDO
347  DO imud=1,nmud
348  ratio_mud(imud,ilayer,ipoin) =
349  & ava0(num_imud_icla(imud))
350  ENDDO
351  ENDDO
352  ENDDO
353 ! INITIALISE THE MUD ARRAYS IF THEY WERE NOT FOUND IN THE
354 ! PREVIOUS FILE
355  DO ipoin = 1,npoin
356  DO ilayer=1,nomblay
357  IF (.NOT.conc_mud_found) THEN
358  conc_mud(ilayer,ipoin) = conc_mud0(ilayer)
359  IF(ipoin==1.AND.ilayer==1)
360  & WRITE(lu,*) 'WARNING: MUD CONCENTRATION NOT FOUND IN ',
361  & 'THE PREVIOUS COMPUTATION FILE, IT IS SET ',
362  & 'FROM THE STEERING FILE'
363  ENDIF
364  IF (.NOT.mtrans_found) THEN
365  trans_mass(ilayer,ipoin) = trans_mass0(ilayer)
366  IF(ipoin==1.AND.ilayer==1)
367  & WRITE(lu,*) 'WARNING: MASS TRANSFER NOT FOUND IN ',
368  & 'THE PREVIOUS COMPUTATION FILE, IT IS SET ',
369  & 'FROM THE STEERING FILE'
370  ENDIF
371  IF (.NOT.toce_mud_found) THEN
372  toce_mud(ilayer,ipoin) = toce_mud0(ilayer)
373  IF(ipoin==1.AND.ilayer==1)
374  & WRITE(lu,*) 'WARNING: MUD TOCE NOT FOUND IN ',
375  & 'THE PREVIOUS COMPUTATION FILE, IT IS SET ',
376  & 'FROM THE STEERING FILE'
377  ENDIF
378  IF (.NOT.partheniades_found) THEN
379  partheniades(ilayer,ipoin) = partheniades0(ilayer)
380  IF(ipoin==1.AND.ilayer==1)
381  & WRITE(lu,*) 'WARNING: PARTHENIADES NOT FOUND IN ',
382  & 'THE PREVIOUS COMPUTATION FILE, IT IS SET ',
383  & 'FROM THE STEERING FILE'
384  ENDIF
385  ENDDO
386  ENDDO
387 !
388  ELSE
389 !
390  DO ipoin=1,npoin
391 ! THE ACTIVE LAYER THICKNESS IS ALWAYS READ FROM THE
392 ! STEERING FILE, AND NOT THE PREVIOUS COMPUTATION RESULTS
393 ! IT WILL BE CALCULATED IN BED1_UPDATE_ACTIVELAYER_HIRANO
394  elay%R(ipoin)= elay0
395  IF (.NOT.conc_mud_found) THEN
396  conc_mud(1,ipoin) = conc_mud0(1)
397  IF(ipoin==1)
398  & WRITE(lu,*) 'WARNING: MUD CONCENTRATION NOT FOUND IN ',
399  & 'THE PREVIOUS COMPUTATION FILE, IT IS SET ',
400  & 'FROM THE STEERING FILE'
401  ENDIF
402  IF (.NOT.mtrans_found) THEN
403  trans_mass(1,ipoin) = trans_mass0(1)
404  IF(ipoin==1)
405  & WRITE(lu,*) 'WARNING: MASS TRANSFER NOT FOUND IN ',
406  & 'THE PREVIOUS COMPUTATION FILE, IT IS SET ',
407  & 'FROM THE STEERING FILE'
408  ENDIF
409  IF (.NOT.toce_mud_found) THEN
410  toce_mud(1,ipoin) = toce_mud0(1)
411  IF(ipoin==1)
412  & WRITE(lu,*) 'WARNING: MUD TOCE NOT FOUND IN ',
413  & 'THE PREVIOUS COMPUTATION FILE, IT IS SET ',
414  & 'FROM THE STEERING FILE'
415  ENDIF
416  IF (.NOT.partheniades_found) THEN
417  partheniades(1,ipoin) = partheniades0(1)
418  IF(ipoin==1)
419  & WRITE(lu,*) 'WARNING: PARTHENIADES NOT FOUND IN ',
420  & 'THE PREVIOUS COMPUTATION FILE, IT IS SET ',
421  & 'FROM THE STEERING FILE'
422  ENDIF
423  DO ilayer = 2,nomblay
424  istrat=ilayer-1
425  IF (.NOT.conc_mud_found) THEN
426  conc_mud(ilayer,ipoin) = conc_mud0(istrat)
427  ENDIF
428  IF (.NOT.mtrans_found) THEN
429  trans_mass(ilayer,ipoin) = trans_mass0(istrat)
430  ENDIF
431  IF (.NOT.toce_mud_found) THEN
432  toce_mud(ilayer,ipoin) = toce_mud0(istrat)
433  ENDIF
434  IF (.NOT.partheniades_found) THEN
435  partheniades(ilayer,ipoin)=partheniades0(istrat)
436  ENDIF
437  ENDDO
438  ENDDO
439  ENDIF
440 !
441  IF(.NOT.debu_mass) THEN
442 ! COMPUTE RATIO_MUD_SAND FROM RATIO_SAND AND RATIO_MUD ISSUES FROM SEDIMENTOLOGICAL FILE
443  DO ipoin=1,npoin
444  DO ilayer=1,nomblay
445  totmud=0.d0
446  totsand=0.d0
447  DO imud=1,nmud
448  totmud= totmud + ratio_mud(imud,ilayer,ipoin)
449  ENDDO
450  DO isand=1,nsand
451  totsand= totsand + ratio_sand(isand,ilayer,ipoin)
452  ENDDO
453  IF(totmud+totsand.GT.0.d0)THEN
454  ratio_mud_sand(ilayer,ipoin)=totmud/(totmud+totsand)
455  ELSE
456  ratio_mud_sand(ilayer,ipoin)=1.d0
457  ENDIF
458  ENDDO
459  ENDDO
460  ENDIF
461 !
462  ENDIF
463 !
464  IF(.NOT.debu.OR..NOT.debu_mass) THEN
465 !
466 ! CHECK THE RATIOS
467 !
468  DO ipoin=1,npoin
469  DO ilayer=1,nomblay
470  check_rs=0.d0
471  check_rm=0.d0
472  IF(nsand.GT.0) THEN
473  DO isand=1,nsand
474  check_rs=check_rs+ratio_sand(isand,ilayer,ipoin)
475  ENDDO
476  IF(abs(check_rs-1.d0).GE.1.d-7) THEN
477  WRITE(lu,*)'SUM OF SAND RATE COEFF MUST BE EQUAL TO 1!'
478  WRITE(lu,*)'VERIFY YOUR MASS RATE OF SAND'
479  CALL plante(1)
480  stop
481  ENDIF
482  ENDIF
483  IF(nmud.GT.0) THEN
484  DO imud=1,nmud
485  check_rm=check_rm+ratio_mud(imud,ilayer,ipoin)
486  ENDDO
487  IF(abs(check_rm-1.d0).GE.1.d-7) THEN
488  WRITE(lu,*)'SUM OF MUD RATE COEFF MUST BE EQUAL TO 1!'
489  WRITE(lu,*)'VERIFY YOUR MASS RATE OF MUD'
490  CALL plante(1)
491  stop
492  ENDIF
493  ENDIF
494  ENDDO
495  ENDDO
496 !
497 ! MASS COMPUTATION : MASS_MUD;MASS_SAND;MASS_MUD_TOT;MASS_SAND_TOT
498 ! ALL MASSES IN [kg/m2]
499 !
500 ! FOR MASS COMPUTATION IS MANDATORY RATIO_MUD_SAND,
501 ! XMVS,CONC_MUD,ES(ILAYER)
502 !
503  IF(nsand.EQ.0)THEN
504  DO ipoin = 1,npoin
505  DO ilayer = 1,nomblay
506  mass_mud_tot(ilayer,ipoin) =
507  & es(ipoin,ilayer)*conc_mud(ilayer,ipoin)
508 ! COMPUTES MASS FOR EVERY MUD
509  DO imud = 1,nmud
510  mass_mud(imud,ilayer,ipoin)=mass_mud_tot(ilayer,ipoin)
511  & *ratio_mud(imud,ilayer,ipoin)
512  ENDDO
513  ENDDO
514  ENDDO
515  ELSE
516  DO ipoin = 1,npoin
517  DO ilayer = 1,nomblay
518  xmvs_lay=0.d0
519  DO isand = 1,nsand ! AVERAGE DENSITY OF SAND FOR THE LAYER
520  xmvs_lay=xmvs_lay+ratio_sand(isand,ilayer,ipoin)
521  & *xmvs0(num_isand_icla(isand))
522  ENDDO
523  term=0.d0
524  IF(nmud.GT.0) THEN
525  term =
526  & ratio_mud_sand(ilayer,ipoin)/conc_mud(ilayer,ipoin)-
527  & (xkv0(ilayer)*(1.d0-ratio_mud_sand(ilayer,ipoin)))/
528  & (xmvs_lay*(1.d0-xkv0(ilayer)))
529 ! TERM REPRESENTS THE DIFFERENCE BETWEEN MUD VOLUME AND VOID
530 ! VOLUME
531  ENDIF
532 ! DISCR IS POSITIVE WHEN MUD FILLS ALL THE SAND POROSITY
533 ! OTHERWISE IS ZERO
534  discr=max(0.d0,term)
535  mass_tot = es(ipoin,ilayer)/
536  & ((1.d0-ratio_mud_sand(ilayer,ipoin))/
537  & (xmvs_lay*(1.d0-xkv0(ilayer)))+discr)
538 !
539  mass_mud_tot(ilayer,ipoin) = ratio_mud_sand(ilayer,ipoin)
540  & *mass_tot
541  mass_sand_tot(ilayer,ipoin) =
542  & (1.d0-ratio_mud_sand(ilayer,ipoin))*mass_tot
543 ! COMPUTES MASS FOR EVERY MUD
544  DO imud = 1,nmud
545  mass_mud(imud,ilayer,ipoin) = mass_mud_tot(ilayer,ipoin)
546  & *ratio_mud(imud,ilayer,ipoin)
547  ENDDO
548 ! COMPUTES MASS FOR EVERY SAND
549  DO isand = 1,nsand
550  mass_sand(isand,ilayer,ipoin) =
551  & mass_sand_tot(ilayer,ipoin)*
552  & ratio_sand(isand,ilayer,ipoin)
553  ENDDO
554  ENDDO
555  ENDDO
556  ENDIF
557 !
558 !
559  ENDIF !(.NOT.DEBU.OR..NOT.DEBU_MASS)
560 !
561 ! REAL MASS COMPUTATION: FROM [kg/m2] to [kg]
562 ! USEFUL FOR FIRST LISTING AND MASS BALANCE
563 !
564  DO icla=1,nsicla
565  masstot(icla)=0.d0
566  ENDDO
567 !
568  IF(nsand.GT.0) THEN
569  DO ipoin=1,npoin
570  DO ilayer=1,nomblay
571  DO isand=1,nsand
572  k=num_isand_icla(isand)
573  masstot(k)=masstot(k)+mass_sand(isand,ilayer,ipoin)
574  & *volu2d%R(ipoin)
575  ENDDO
576  ENDDO
577  ENDDO
578  ENDIF
579 !
580  IF(nmud.GT.0) THEN
581  DO ipoin=1,npoin
582  DO ilayer=1,nomblay
583  DO imud=1,nmud
584  k=num_imud_icla(imud)
585  masstot(k)=masstot(k)+mass_mud(imud,ilayer,ipoin)
586  & *volu2d%R(ipoin)
587  ENDDO
588  ENDDO
589  ENDDO
590  ENDIF
591 !
592  IF(ncsize.GT.1) THEN
593  DO icla=1,nsicla
594  masstot(icla) = p_sum(masstot(icla))
595  ENDDO
596  ENDIF
597 !
598  DO icla=1,nsicla
599  mass0tot(icla)=masstot(icla)
600  ENDDO
601 !
602  WRITE(lu,*)
603  WRITE(lu,100)
604  DO icla=1,nsicla
605  WRITE(lu,111) icla,masstot(icla)
606  ENDDO
607 100 FORMAT(20x,'GAIA INITIAL MASS OF SEDIMENTS: ')
608 111 FORMAT(1x,'TOTAL MASS OF CLASS ',i2,' =',g16.7,'( KG )')
609 !
610 !-----------------------------------------------------------------------
611  IF(.NOT.debu) THEN
612  DEALLOCATE(estratum)
613  DEALLOCATE(ratio_init)
614  ENDIF
615 !-----------------------------------------------------------------------
616 !
617  RETURN
618  END
double precision, dimension(:,:,:), allocatable, target mass_mud
Surface mass of mud (kg/m2), for imud,ilayer,ipoin.
subroutine bed1_init_sediment_gaia(NSICLA, ELAY, ZF, ZR, NPOIN, XMVS0, ES, NOMBLAY, DEBU, VOLU2D, NUMSTRAT, MAXVAR)
double precision, dimension(:,:,:), allocatable, target mass_sand
Surface mass of sand (kg/m2), for isand,ilayer,ipoin.
subroutine user_bed_init(NUMSTRAT, NPOIN, NSICLA, ESTRATUM, RATIO_INIT)
Definition: user_bed_init.f:7
Definition: bief.f:3