The TELEMAC-MASCARET system  trunk
bedload_main_gaia.f
Go to the documentation of this file.
1 ! ****************************
2  SUBROUTINE bedload_main_gaia
3 ! ****************************
4 !
5  &(acladm,ksp,ksr, v2dpar,unsv2d,cf,ebor,fw,hn,liqbor,
6  & mask, maskel, maskpt, qbor, u2d,
7  & v2d, s,unladm,uw,thetaw,mu,tob,tobw,tw,zf,
8  & debug, hidfac, icf, ielmt, kddl, kdir,
9  & kent, klog, kneu, ksort,
10  & npoin, nptfr, nsicla, optban, beta, dcla,
11  & grav, hidi, hmin, vce, xmve, xmvs0, xwc,
12  & pi, karman, zero, karim_holly_yang,msk, susp, vf,
13  & mesh, liebor, limtec, masktr,
14  & it1, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11,
15  & t12,t13,unorm,ac, dt,
16  & breach, calfa_cl, coefpn,
17  & hiding, qscl_c, qscl_s, qs_c,
18  & qsclxc, qsxc, qsclyc, qsyc, salfa_cl,
19  & entets, seccurrent, slopeff,
20  & phised, devia, beta2, bijk,sedco,houle,
21  & u3d,v3d,code,flbcla,maxadv,ratio_sand,h_tel,
22  & hw, thetac, tobcw_mean, tobcw_max, cstaeq)
23 !
24 !***********************************************************************
25 ! GAIA
26 !***********************************************************************
27 !
30 !
31 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
141 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
142 !
143  USE bief
144  USE interface_gaia, ex_bedload_main => bedload_main_gaia
146  & num_icla_isand,evcl_mb,evol_mb,min_sed_mass_comp,
147  & mass_sand_tot,ratio_mud_sand,
148  & num_icla_imud,mass_mud,mass_sand_active_layer,
149  & mass_sand_masked, evcl_m_tot_sand,
150  & ratio_evol_tot_sand,mofac_bed,mudb,f_mudb,
151  & nestor,xkv0,volu2d,massnestor,
152  & nmud,num_isand_icla,num_imud_icla,
153  & sanfra
155  USE interface_parallel, ONLY: p_sum
156  IMPLICIT NONE
157 !
158 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
159 !
160  TYPE(bief_obj), INTENT(IN) :: ACLADM, KSR,V2DPAR,UNSV2D
161  TYPE(bief_obj), INTENT(IN) :: CF,FW,KSP,HN,LIQBOR
162  TYPE(bief_obj), INTENT(IN) :: MASK, MASKEL, MASKPT
163  TYPE(bief_obj), INTENT(IN) :: QBOR
164  TYPE(bief_obj), INTENT(INOUT) :: U2D,V2D,TOB,MU,UNORM,EBOR
165  TYPE(bief_obj), INTENT(IN) :: S,UNLADM
166  TYPE(bief_obj), INTENT(IN) :: UW, THETAW, TOBW, TW
167  TYPE(bief_obj), INTENT(IN) :: HW, THETAC
168  TYPE(bief_obj), INTENT(INOUT) :: TOBCW_MEAN, TOBCW_MAX
169  TYPE(bief_obj), INTENT(IN) :: ZF
170  INTEGER, INTENT(IN) :: DEBUG, HIDFAC, ICF,MAXADV
171  INTEGER, INTENT(IN) :: IELMT, KDDL, KDIR, KENT
172  INTEGER, INTENT(IN) :: KLOG, KNEU, KSORT
173  INTEGER, INTENT(IN) :: NPOIN, NPTFR
174  INTEGER, INTENT(IN) :: NSICLA, OPTBAN
175  DOUBLE PRECISION, INTENT(IN) :: BETA
176  DOUBLE PRECISION, INTENT(IN) :: DCLA(nsicla),GRAV
177  DOUBLE PRECISION, INTENT(IN) :: HIDI(nsicla),HMIN,VCE
178  DOUBLE PRECISION, INTENT(IN) :: XMVE,XMVS0(nsicla),XWC(nsicla)
179  DOUBLE PRECISION, INTENT(IN) :: PI,KARMAN,ZERO
180  DOUBLE PRECISION, INTENT(IN) :: KARIM_HOLLY_YANG
181  LOGICAL, INTENT(IN) :: MSK, SUSP, VF
182  LOGICAL, INTENT(IN) :: SECCURRENT
183  LOGICAL, INTENT(IN) :: SEDCO(nsicla),HOULE
184  TYPE(bief_mesh), INTENT(INOUT) :: MESH
185  TYPE(bief_obj), INTENT(INOUT) :: FLBCLA
186  TYPE(bief_obj), INTENT(INOUT) :: LIEBOR, LIMTEC, MASKTR
187  TYPE(bief_obj), INTENT(INOUT) :: IT1,T1,T2,T3,T4,T5,T6,T7
188  TYPE(bief_obj), INTENT(INOUT) :: T8,T9,T10,T11,T12,T13
189  DOUBLE PRECISION, INTENT(INOUT) :: AC(nsicla), DT
190  DOUBLE PRECISION, INTENT(INOUT) :: RATIO_SAND(nsand,nomblay,npoin)
191  TYPE(bief_obj), INTENT(INOUT) :: BREACH, CALFA_CL, COEFPN
192  TYPE(bief_obj), INTENT(INOUT) :: HIDING
193  TYPE(bief_obj), INTENT(INOUT) :: QSCL_C,QSCL_S
194  TYPE(bief_obj), INTENT(INOUT) :: QS_C, QSCLXC, QSXC, QSCLYC
195  TYPE(bief_obj), INTENT(INOUT) :: QSYC, SALFA_CL
196  LOGICAL, INTENT(INOUT) :: ENTETS
197  DOUBLE PRECISION, INTENT(IN) :: BETA2, PHISED
198  INTEGER, INTENT (IN) :: SLOPEFF, DEVIA
199  DOUBLE PRECISION, INTENT(IN) :: BIJK
200  TYPE(bief_obj), INTENT(IN) :: U3D,V3D
201  CHARACTER(LEN=24), INTENT(IN) :: CODE
202  TYPE(bief_obj), INTENT(IN) :: H_TEL
203  TYPE(bief_obj), INTENT(INOUT) :: CSTAEQ
204 !
205 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
206 !
207  INTEGER I,K,IPOIN,ICLA,IMUD,ISAND
208  DOUBLE PRECISION, ALLOCATABLE :: TMP_RATIO(:)
209 !
210 !======================================================================!
211 !======================================================================!
212 ! PROGRAM !
213 !======================================================================!
214 !======================================================================!
215 !
216 !
217 !
218 ! Calculation of sand fraction content at each node (Wilcock and Crowe, 2003)
219  IF (icf == 10) THEN
220  DO k = 1, npoin
221  sanfra(k) = 0.0d0
222  DO i = 1, nsand
223  IF (dcla(i).LT.2.d-3) THEN
224  sanfra(k) = sanfra(k) + ratio_sand(i,1,k)
225  ENDIF
226  ENDDO
227  ENDDO
228  ENDIF
229 !
230 
231 ! INITIALISES TECHNICAL BOUNDARY CONDITIONS
232 !
233  IF (debug > 0) WRITE(lu,*) 'BEDLOAD_DIFFIN_GAIA'
235  & (u2d, v2d, mesh%NBOR, mesh%XNEBOR, mesh%YNEBOR,
236  & maskel, mesh%NELBOR, nptfr, kent, ksort, klog,
237  & kdir, kddl, kneu, msk, it1, liebor, masktr,limtec,
238  & mesh%IKLBOR%I,mesh%NELEB,mesh%NELEBX)
239  IF (debug > 0) WRITE(lu,*) 'END_BEDLOAD_DIFFIN'
240 !
241  ! Memory optimisation (from intel debug)
242  ALLOCATE(tmp_ratio(npoin))
243  DO i = 1, nsicla
244  k=num_icla_isand(i)
245 !
246 ! FOR SAND
247  IF(.NOT.sedco(i)) THEN
248  IF (debug > 0) WRITE(lu,*)
249  & 'BEDLOAD_SOLIDISCHARGE_GAIA : ',i,'/',nsicla
250  tmp_ratio = ratio_sand(k,1,1:npoin)
252  & (mesh, u2d, v2d, unorm,hn, tw, uw, mu,tob,cf,
253  & tobw,fw,thetaw,
254  & tmp_ratio,
255  & maskpt, maskel, acladm,
256  & unladm,ksp,ksr, liqbor, debug, npoin,
257  & nptfr, ielmt, icf, kent, optban, hidfac, grav,
258  & dcla(i), xwc(i), xmve, xmvs0(i), vce, hmin,
259  & hidi(i),karman,zero,pi,
260  & karim_holly_yang,susp,msk,t1,t2,
261  & t3, t4, t5, t6, t7, t8, t9, t10, t11,t12, ac(i),
262  & hiding,qscl_c%ADR(i)%P,qscl_s%ADR(i)%P,
263  & slopeff,coefpn,phised,
264  & calfa_cl%ADR(i)%P,salfa_cl%ADR(i)%P,
265  & beta,zf,s,
266  & devia, beta2 , seccurrent, bijk,houle,unsv2d,
267  & u3d,v3d,code, h_tel,
268  & hw,thetac,tobcw_mean,tobcw_max,cstaeq%ADR(i)%P,
269  & sanfra)
270  IF(debug > 0) WRITE(lu,*) 'END_BEDLOAD_SOLIDISCHARGE'
271 !
272  ELSE
273 ! FOR COHESIVE SEDIMENT: ZERO BEDLOAD TRANSPORT RATE
274  CALL os('X=0 ',x=qscl_c%ADR(i)%P)
275  CALL os('X=0 ',x=qsclxc%ADR(i)%P)
276  CALL os('X=0 ',x=qsclyc%ADR(i)%P)
277  ENDIF
278 !
279 ! MORPHOLOGICAL FACTOR on evolution bed is applicated
280 ! on QSC directly to take in account rigid bed
281 ! After evolution of bed, QSC take the good value
282  DO ipoin=1,npoin
283  qscl_c%ADR(i)%P%R(ipoin)=
284  & qscl_c%ADR(i)%P%R(ipoin)*mofac_bed
285  ENDDO
286 !
287  ENDDO
288 !
289 !
290  DO ipoin=1,npoin
291  evcl_m_tot_sand(ipoin)= 0.d0
292  ENDDO
293 !
294 ! COMPUTES THE EVOLUTION FOR EACH CLASS
295 !
296  DO i = 1, nsicla
297 !
298  k=num_icla_isand(i)
299 !
300  IF(.NOT.sedco(i)) THEN
301 !
302  DO ipoin=1,npoin
303 ! IF MORE THAN 30% OF MUD, BEDLOAD IS NOT ACTIVE
304  IF(ratio_mud_sand(1,ipoin).LE.0.3d0)THEN
305  mass_sand_active_layer(ipoin)= mass_sand(k,1,ipoin)
306  mass_sand_masked(ipoin)=0.d0
307  ELSE
308  mass_sand_active_layer(ipoin)= 0.d0
309  mass_sand_masked(ipoin)= mass_sand(k,1,ipoin)
310  ENDIF
311  ENDDO
312 !
313  IF (debug>0) WRITE(lu,*) 'BEDLOAD_EVOL_GAIA : ',i,'/',nsicla
314  tmp_ratio = ratio_sand(k,1,1:npoin)
315  CALL bedload_evol_gaia(s,coefpn,calfa_cl%ADR(i)%P,
316  & salfa_cl%ADR(i)%P, limtec,
317  & ebor%ADR(k)%P,maskel,mask,
318  & v2dpar,unsv2d,debug,npoin,nptfr,ielmt,
319  & kent,kdir,kddl,dt,xmvs0(i),
320  & vf,entets,msk,
321  & mesh,qscl_c%ADR(i)%P,
322  & t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,
323  & t13,breach,qsclxc%ADR(i)%P,
324  & qsclyc%ADR(i)%P,slopeff,
325  & i,flbcla,liqbor,qbor%ADR(i)%P,maxadv,
326  & mass_sand_active_layer,
327  & tmp_ratio,evcl_mb%ADR(i)%P)
328  IF(debug.GT.0) WRITE(lu,*) 'END_BEDLOAD_EVOL'
329 !
330  DO ipoin=1,npoin
331  mass_sand(k,1,ipoin)= mass_sand_active_layer(ipoin)
332  & + mass_sand_masked(ipoin)
333  ENDDO
334 !
335  DO ipoin=1,npoin
336  evcl_m_tot_sand(ipoin)= evcl_m_tot_sand(ipoin)+
337  & evcl_mb%ADR(i)%P%R(ipoin)
338  ENDDO
339 !
340  ENDIF
341  ENDDO
342  DEALLOCATE(tmp_ratio)
343 !
344 !============ run Nestor ==============================================!
345 !
346  IF(nestor) THEN
347  IF(nsand.GT.0) THEN
348  massnestor = 0.d0
349  DO ipoin=1,npoin
350  DO isand=1,nsand
351  k = num_isand_icla(isand)
352  massnestor(k) = massnestor(k) +
353  & evcl_mb%ADR(k)%P%R(ipoin) *volu2d%R(ipoin)
354  ENDDO
355  ENDDO
356  ENDIF
357  IF(nmud.GT.0) THEN
358  DO ipoin=1,npoin
359  DO imud=1,nmud
360  k = num_imud_icla(imud)
361  massnestor(k) = massnestor(k) +
362  & evcl_mb%ADR(k)%P%R(ipoin) *volu2d%R(ipoin)
363  ENDDO
364  ENDDO
365  ENDIF
366  IF(ncsize.GT.1) THEN
367  DO icla=1,nsicla
368  massnestor(icla) = p_sum(massnestor(icla))
369  ENDDO
370  ENDIF
371 !
372  CALL nestor_interface_gaia(2,123,xmvs0,xkv0(1),volu2d)
373 !
374  IF(nsand.GT.0) THEN
375  t8%R = 0.d0
376  DO ipoin=1,npoin
377  DO isand=1,nsand
378  k = num_isand_icla(isand)
379  t8%R(k) = t8%R(k) +
380  & evcl_mb%ADR(k)%P%R(ipoin) *volu2d%R(ipoin)
381  ENDDO
382  ENDDO
383  ENDIF
384  IF(nmud.GT.0) THEN
385  DO ipoin=1,npoin
386  DO imud=1,nmud
387  k = num_imud_icla(imud)
388  t8%R(k) = t8%R(k) +
389  & evcl_mb%ADR(k)%P%R(ipoin) *volu2d%R(ipoin)
390  ENDDO
391  ENDDO
392  ENDIF
393  IF(ncsize.GT.0) THEN
394  DO icla = 1,nsicla
395  t8%R(icla) = p_sum(t8%R(icla))
396  massnestor(icla) = t8%R(icla) - massnestor(icla)
397  END DO
398  ENDIF
399  ENDIF
400 !
401 !======================================================================!
402 !
403 ! EVOLUTIONS AND QS FOR EACH CLASS ARE ADDED
404 !
405 ! INITIALISES
406 !
407  CALL os('X=0 ',x=qs_c)
408  CALL os('X=0 ',x=qsxc)
409  CALL os('X=0 ',x=qsyc)
410  CALL os('X=0 ',x=mudb)
411  CALL os('X=0 ',x=f_mudb)
412 ! FIXME: COMPUTE EVOL_MB IS NOT REALLY USEFULL BECAUSE MASS BALANCE
413 ! IS THEN COMPUTED FOR EACH CLASS AND THE SUM IS LOCALLY COMPUTED
414  CALL os('X=0 ',x=evol_mb)
415 !
416 !
417  DO i=1,nsicla
418 ! After evolution of bed, QSC take the good value
419 ! see before for multiplcation by MOFAC_BED
420  DO ipoin=1,npoin
421  qscl_c%ADR(i)%P%R(ipoin)=
422  & qscl_c%ADR(i)%P%R(ipoin)/mofac_bed
423  ENDDO
424 !
425  k=num_icla_isand(i)
426  IF(.NOT.sedco(i)) THEN
427  CALL os('X=X+Y ', x=qs_c, y=qscl_c%ADR(i)%P)
428  CALL os('X=X+YZ ', x=qsxc, y=qscl_c%ADR(i)%P,
429  & z=calfa_cl%ADR(i)%P)
430  CALL os('X=X+YZ ', x=qsyc, y=qscl_c%ADR(i)%P,
431  & z=salfa_cl%ADR(i)%P)
432  CALL os('X=X+Y ', x=evol_mb, y=evcl_mb%ADR(i)%P)
433  ENDIF
434  ENDDO
435 !
436 ! UPDATE OF SAND MASS WITH EVOLUTION DUE TO BEDLOAD
437 ! (OF THE FIRST LAYER)
438 !
439  DO i=1,nsicla
440  k=num_icla_isand(i)
441  IF(.NOT.sedco(i)) THEN
442  DO ipoin=1,npoin
443  mass_sand(k,1,ipoin) = mass_sand(k,1,ipoin) +
444  & evcl_mb%ADR(i)%P%R(ipoin)
445  ENDDO
446  ENDIF
447  ENDDO
448 !
449 ! POURCENTAGE OF MASS SAND IN THE ACTIVE LAYER THAT HAVE MOVED
450  DO ipoin=1,npoin
451 ! MASS_SAND_TOT IS COMING FROM BED UPDATE ACTIVE HIRANO BEFORE BEDLOAD
452  IF(mass_sand_tot(1,ipoin).GE.min_sed_mass_comp)THEN
453  ratio_evol_tot_sand(ipoin)= evcl_m_tot_sand(ipoin)/
454  & mass_sand_tot(1,ipoin)
455  ratio_evol_tot_sand(ipoin)=
456  & min(ratio_evol_tot_sand(ipoin),1.d0)
457  ratio_evol_tot_sand(ipoin)=
458  & max(ratio_evol_tot_sand(ipoin),-1.d0)
459  ELSE
460  ratio_evol_tot_sand(ipoin)= 0.d0
461  ENDIF
462  ENDDO
463 ! THE SAME RATIO OF MUD MOVED AND IS ADDED IN SUSPENSION FLUX
464  DO i = 1, nsicla
465  k=num_icla_imud(i)
466  IF(sedco(i)) THEN
467  DO ipoin =1,npoin
468  IF(ratio_evol_tot_sand(ipoin).LE.0.d0)THEN
469  mudb%ADR(i)%P%R(ipoin)= ratio_evol_tot_sand(ipoin)*
470  & mass_mud(k,1,ipoin)
471  ELSE
472  mudb%ADR(i)%P%R(ipoin) = 0.d0
473  ENDIF
474  mass_mud(k,1,ipoin)= mass_mud(k,1,ipoin)+
475  & mudb%ADR(i)%P%R(ipoin)
476  f_mudb%ADR(i)%P%R(ipoin)=
477  & - mudb%ADR(i)%P%R(ipoin)/dt
478 ! F_MUDB IS ADDED FOR SUPENSION FLUX AFTER
479 ! IN BED1_SUSPENSION_ERODE
480  ENDDO
481  ENDIF
482  ENDDO
483 !
484 ! TIDAL FLATS WITH MASKING
485 !
486  IF(optban.EQ.2) CALL os('X=XY ',x=evol_mb,y=maskpt)
487 !
488 !======================================================================!
489 !======================================================================!
490 !
491  RETURN
492  END
subroutine nestor_interface_gaia(OPTION, GRAFCOUNT, XMVS0, XKV01, VOLU2D)
integer, target nomblay
Number of bed load model layers = NUMSTRAT+1 to take the active layer into account.
subroutine bedload_diffin_gaia(U2D, V2D, NBOR, XNEBOR, YNEBOR, MASKEL, NELBOR, NPTFR, KENT, KSORT, KLOG, KDIR, KDDL, KNEU, MSK, CLT, LITBOR, MASKTR, LIMTRA, IKLBOR, NELEB, NELEBX)
double precision, dimension(:,:,:), allocatable, target mass_sand
Surface mass of sand (kg/m2), for isand,ilayer,ipoin.
subroutine breach
Definition: breach.f:4
integer nsand
Total number of sand.
subroutine bedload_main_gaia(ACLADM, KSP, KSR, V2DPAR, UNSV2D, CF, EBOR, FW, HN, LIQBOR, MASK, MASKEL, MASKPT, QBOR, U2D, V2D, S, UNLADM, UW, THETAW, MU, TOB, TOBW, TW, ZF, DEBUG, HIDFAC, ICF, IELMT, KDDL, KDIR, KENT, KLOG, KNEU, KSORT, NPOIN, NPTFR, NSICLA, OPTBAN, BETA, DCLA, GRAV, HIDI, HMIN, VCE, XMVE, XMVS0, XWC, PI, KARMAN, ZERO, KARIM_HOLLY_YANG, MSK, SUSP, VF, MESH, LIEBOR, LIMTEC, MASKTR, IT1, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12, T13, UNORM, AC, DT, BREACH, CALFA_CL, COEFPN, HIDING, QSCL_C, QSCL_S, QS_C, QSCLXC, QSXC, QSCLYC, QSYC, SALFA_CL, ENTETS, SECCURRENT, SLOPEFF, PHISED, DEVIA, BETA2, BIJK, SEDCO, HOULE, U3D, V3D, CODE, FLBCLA, MAXADV, RATIO_SAND, H_TEL, HW, THETAC, TOBCW_MEAN, TOBCW_MAX, CSTAEQ)
subroutine bedload_evol_gaia(S, COEFPN, CALFA, SALFA, LIMTEC, EBOR, MASKEL, MASK, V2DPAR, UNSV2D, DEBUG, NPOIN, NPTFR, IELMT, KENT, KDIR, KDDL, DT, XMVS, VF, ENTETS, MSK, MESH, QS, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12, T13, BREACH, QSX, QSY, SLOPEFF, ICLA, FLBCLA, LIQBOR, QBOR, MAXADV, MASS_SAND, RATIO_SAND, EVCL_MB)
subroutine bedload_solidischarge_gaia(MESH, U2D, V2D, UNORM, HN, TW, UW, MU, TOB, CF, TOBW, FW, THETAW, RATIO_SAND, MASKPT, MASKEL, ACLADM, UNLADM, KSP, KSR, LIQBOR, DEBUG, NPOIN, NPTFR, IELMT, ICF, KENT, OPTBAN, HIDFAC, GRAV, DCLA, XWC, XMVE, XMVS, VCE, HMIN, HIDI, KARMAN, ZERO, PI, K_H_Y, SUSP, MSK, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12, AC, HIDING, QSC, QSS, SLOPEFF, COEFPN, PHISED, CALFA, SALFA, BETA, ZF, S, DEVIA, BETA2, SECCURRENT, BIJK, HOULE, UNSV2D, U3D, V3D, CODE, H_TEL, HW, THETAC, TOBCW_MEAN, TOBCW_MAX, CSTAEQ, SANFRA)
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
Definition: bief.f:3