The TELEMAC-MASCARET system  trunk
bedload_nerbed_vf_gaia.f
Go to the documentation of this file.
1 ! ******************************
2  SUBROUTINE bedload_nerbed_vf_gaia
3 ! ******************************
4 !
5  &(mesh,liebor,ksort,v2dpar,qsx,qsy,npoin,nseg,
6  & nptfr,dt,qs,t1,t2,t3,breach,nubo,vnoin,mass_sand)
7 !
8 !***********************************************************************
9 ! GAIA
10 !***********************************************************************
11 !
13 !
14 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
33 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
34 !
35  USE interface_gaia, ex_bedload_nerbed_vf => bedload_nerbed_vf_gaia
36  USE bief
38  IMPLICIT NONE
39 !
40 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
41 !
42  TYPE(bief_mesh), INTENT(INOUT) :: MESH
43  TYPE(bief_obj), INTENT(IN) :: LIEBOR
44  TYPE(bief_obj), INTENT(IN) :: QSX, QSY
45  INTEGER, INTENT(IN) :: NPOIN, NSEG, NPTFR,KSORT
46  DOUBLE PRECISION, INTENT(IN) :: DT
47  TYPE(bief_obj), INTENT(INOUT) :: QS, T1, T2, T3
48  TYPE(bief_obj), INTENT(INOUT) :: BREACH
49  DOUBLE PRECISION, INTENT(IN) :: V2DPAR(npoin)
50  INTEGER, INTENT(IN) :: NUBO(2,nseg)
51  DOUBLE PRECISION, INTENT(IN) :: VNOIN(3,nseg)
52  DOUBLE PRECISION, INTENT(IN) :: MASS_SAND(npoin)
53 !
54 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
55 !
56  INTEGER :: I, K
57  INTEGER :: IEL, IEL1, IEL2, ISEGIN
58  DOUBLE PRECISION :: QSP1, QSP2, QSPC
59  DOUBLE PRECISION :: XN, YN, TEMP,PROD_SCAL
60  DOUBLE PRECISION :: VNOIN1, VNOIN2, RNORM
61 !
62 !======================================================================!
63 !======================================================================!
64 ! PROGRAM !
65 !======================================================================!
66 !======================================================================!
67 !
68  ! ****************** !
69  ! I - INITIALISES !
70  ! ****************** !
71 !
72  ! BREACH INDICATES IF NON ERODABLE BED WILL BE REACHED
73  ! DURING TIME STEP FOR THIS POINT
74  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
75 !
76 ! GIVE T1 AND T2 THE SAME STRUCTURE AS QS
77 !
78  CALL cpstvc(qs,t1)
79  CALL cpstvc(qs,t2)
80 !
81  DO iel = 1, npoin
82  breach%I(iel) = 0
83  t1%R(iel)=0.d0
84  t2%R(iel)=0.d0
85  ENDDO
86 !
87  ! ************************************************* !
88  ! II - DETERMINES THE OUTGOING FLUX FOR EACH CELL ! (_IMP_)
89  ! ************************************************* !
90  ! THE PRINCIPLE IS THAT QS IS CALCULATED FOR EACH SEGMENT AS
91  ! HALF THE SUM OF THE QS AT THE CENTERS OF THE ELEMENTS WHICH
92  ! SEGMENT FORMS THE BOUNDARY. IT IS THEN PROJECTED ON THE NORMAL
93  ! TO THE SEGMENT, MULTIPLIED BY THE LENGTH OF THE SEGMENT, AND
94  ! THIS FLUX TERM IS ADDED (OR SUBTRACTED) TO THE TWO ELEMENTS.
95 !
96  DO isegin = 1, nseg
97 !
98  iel1 = nubo(1,isegin)
99  iel2 = nubo(2,isegin)
100 !
101  ! II.1 - SEGMENT LENGTH (RNORM)
102  ! ----------------------------------
103  vnoin1 = vnoin(1,isegin)
104  vnoin2 = vnoin(2,isegin)
105  rnorm = vnoin(3,isegin)
106  prod_scal= (mesh%X%R(iel2)-mesh%X%R(iel1))*vnoin1+
107  & (mesh%Y%R(iel2)-mesh%Y%R(iel1))*vnoin2
108  IF(prod_scal.LT.0.d0)THEN
109  iel1 = nubo(2,isegin)
110  iel2 = nubo(1,isegin)
111  ENDIF
112 !
113  ! II.2 - PROJECTS QS FOR THE SEGMENT ONTO THE SEGMENT NORMAL
114  ! ------------------------------------------------------------
115  qsp1 = vnoin1*qsx%R(iel1) + vnoin2*qsy%R(iel1)
116  qsp2 = vnoin1*qsx%R(iel2) + vnoin2*qsy%R(iel2)
117  qspc = (qsp1+qsp2)*0.5d0
118 !
119  ! II.3 - QS SUCH AS THE OUTGOING FLUX IS MAXIMUM
120  ! ----------------------------------------------
121  t1%R(iel1) = t1%R(iel1) + rnorm*max(qspc,qsp1,0.d0)
122  t1%R(iel2) = t1%R(iel2) - rnorm*min(qspc,qsp2,0.d0)
123 !
124  IF(qspc > 0.d0) THEN
125  t2%R(iel1) = t2%R(iel1) + rnorm*qsp1
126  ELSEIF(qspc < 0.d0) THEN
127  t2%R(iel2) = t2%R(iel2) - rnorm*qsp2
128  ENDIF
129 !
130  ENDDO
131 !
132  ! ************************************** !
133  ! III - LOOP ON THE BOUNDARY NODES !
134  ! ************************************** !
135 !
136  DO k = 1, nptfr
137  iel = mesh%NBOR%I(k)
138 !
139  ! III.1 - FREE EVOLUTION: SEDIMENTS ARE FREE TO LEAVE
140  ! ---------------------------------------------------------
141  IF (liebor%I(k) == ksort) THEN
142 !
143  ! XNEBOR (*+NPTFR) AND YNEBOR (*+NPTFR)
144  ! CONTAIN THE VECTOR NORMAL TO A BOUNDARY NODE
145  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
146  xn = mesh%XNEBOR%R(k+nptfr)
147  yn = mesh%YNEBOR%R(k+nptfr)
148  temp = qsx%R(iel)*xn + qsy%R(iel)*yn
149  IF (temp > 0.d0) THEN
150  t1%R(iel) = t1%R(iel) + temp
151  t2%R(iel) = t2%R(iel) + temp
152  ENDIF
153 !
154  ENDIF
155 !
156  ! III.2 - FOR A SOLID BOUNDARY: NOTHING TO PROGRAM
157  ! BECAUSE THE SEDIMENT FLUX IS ZERO HERE
158  ! FOR IMPOSED EVOLUTION : SEE BEDLOAD_SOLVS_VF_GAIA.F
159  ! --------------------------------------------------------
160  ENDDO
161 !
162  IF(ncsize > 1) THEN
163  CALL parcom(t1, 2, mesh)
164  CALL parcom(t2, 2, mesh)
165  ENDIF
166 !
167  ! ************************************************ !
168  ! IV - COMPUTES THE MAXIMUM FLUX AUTHORISED PER CELL!
169  ! ************************************************ !
170 !
171  DO i = 1, npoin
172 !
173  t3%R(i)=mass_sand(i)*v2dpar(i)/dt
174  IF (t3%R(i) < 0.d0) t3%R(i) = 0.d0
175 !
176  ! IF THE OUTGOING FLUX IS TOO LARGE, QS IS CAPPED AT THE NODE
177  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
178  IF(t1%R(i) > t3%R(i)) THEN
179  breach%I(i) = 1
180  IF(t2%R(i) > t3%R(i)) THEN
181  qs%R(i) = qs%R(i)*t3%R(i)/t2%R(i)
182  ENDIF
183  ENDIF
184 !
185  ENDDO
186 !
187 !======================================================================!
188 !======================================================================!
189 !
190  RETURN
191  END
192 
subroutine bedload_nerbed_vf_gaia(MESH, LIEBOR, KSORT, V2DPAR, QSX, QSY, NPOIN, NSEG, NPTFR, DT, QS, T1, T2, T3, BREACH, NUBO, VNOIN, MASS_SAND)
subroutine breach
Definition: breach.f:4
subroutine cpstvc(X, Y)
Definition: cpstvc.f:7
subroutine parcom(X, ICOM, MESH)
Definition: parcom.f:7
Definition: bief.f:3