The TELEMAC-MASCARET system  trunk
cvsp_rm_fraction.f
Go to the documentation of this file.
1 ! ***************************
2  SUBROUTINE cvsp_rm_fraction
3 ! ***************************
4 !
5  &(j,i,dzfcl)
6 !
7 !***********************************************************************
8 ! SISYPHE V7P2 16/05/2017
9 !***********************************************************************
10 !
11 !brief REMOVES (PARTS) OF A FRACTION AFTER EROSION
12 !+ FROM THE VERTICAL SORTING PROFILE;
13 !
14 !history UWE MERKEL
15 !+ 02/02/2012
16 !+ V6P2
17 !+
18 !
19 !history P. A. TASSI (EDF R&D, LNHE)
20 !+ 12/03/2013
21 !+ V6P3
22 !+ Cleaning, cosmetic
23 !
24 !
25 !history U. MERKEL (UHM), R. KOPMANN (BAW)
26 !+ 12/06/2016 / 2017
27 !+ V6P3 / V7P2
28 !+ Update
29 !
30 !history R. KOPMANN (BAW)
31 !+ 25/02/2019
32 !+ V7P2
33 !+ Removing 1/NSICLA
34 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
35 !| J |<--| INDEX OF A POINT IN MESH
36 !| I |<--| INDEX OF A FRACTION
37 !| DZFCL |<--| VALUE OF A FRACTION IN CM !
38 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
39 !
43 !
45  IMPLICIT NONE
46 !
47 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
48 !
49  INTEGER, INTENT(IN) :: J
50  INTEGER, INTENT(IN) :: I
51  DOUBLE PRECISION, INTENT(IN) :: DZFCL
52 !
53 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
54 !
55  DOUBLE PRECISION EROSUM,PROF,PROTH,PROV,REST1,DFRAC
56  DOUBLE PRECISION EROSTRENGTH,ERODEPTH
57  DOUBLE PRECISION PRO_D_LOW,PRO_F_LOW,PROV_TOTAL,F1,F2,D1,D2
58 !
59  INTEGER CNTR, CNTRB, II, KK, LOWPNT, JG, K
60  LOGICAL RET
61  DOUBLE PRECISION ST_S, ST_A, ST_M, ST_C
62  CHARACTER(LEN=7) OCSTR
63 !
64 !-----------------------------------------------------------------------
65 !
66  jg = j
67  IF (ncsize > 1) jg = mesh%KNOLG%I(j)
68 
69  WRITE(unit=ocstr, fmt='(A,I2)') 'RFTIN',i
70 
71  IF (cvsp_db(jg,0)) CALL cvsp_p('./',ocstr,jg)
72 !
73  erosum = - dzfcl !EROSION SUM POSITIV = EROSION
74  cntr = -1 !SECTION CNTR BELOW PRO_MAX
75 !
76 !-----------------------------------------------------------------------
77 ! LOOP THROUGH VSP SECTIONS UNTIL ALL EROSION IS DONE
78 ! COUNTING BACK EROSUM -> 0
79 ! IDENTIFY THE LAST SECTION POINT AFFECTED
80 !-----------------------------------------------------------------------
81 !
82  DO WHILE ((erosum > 0.d0.OR.cntr==-1).AND.(cntr.LT.pro_max(j)-2))
83  cntr = cntr + 1
84 !
85  f1 = pro_f(j,pro_max(j)-cntr,i)
86  f2 = pro_f(j,pro_max(j)-(cntr+1),i)
87 !
88  d1 = pro_d(j,pro_max(j)-cntr,i)
89  d2 = pro_d(j,pro_max(j)-(cntr+1),i)
90 !
91  prof = (f1 + f2) * 0.5d0 ! MEAN FRACTION
92  proth = d1 - d2 ! PROFILE SECTION THICKNESS
93  prov = prof * proth ! PROFILE SECTION VOLUME ( EQ. THICKNESS OF FRACTION)
94  dfrac = f1 - f2 ! DELTA FRACTION = CHANGE OF FRACTION OVER THIS SECTION
95 !
96  IF (proth.GT.0.d0) THEN
97 !
98 !-----------------------------------------------------------------------
99 ! CASE 2
100 ! THIS VSP SECTION HAS MORE THEN ENOUGH MATERIAL
101 ! THEN SPLIT IT AT THE DEPTH OF MAXIMUM EROSION,
102 ! TO KEEP ANTHING BELOW UNCHANGED
103 ! => MAKES CASE 2 TO CASE 1
104 !-----------------------------------------------------------------------
105 !
106  IF(erosum < prov) THEN
107  IF(cvsp_db(jg,0)) CALL cvsp_p('./VSP/','RFT_C2_',jg)
108 !
109 !-----------------------------------------------------------------------
110 ! MAXIMUM DEPTH OF EROSION: SOLVING QUADRATIC PROBLEM BY CROSS-MULTIPLICATION
111 ! LINEAR OR QUATRATIC ? PREVENTS RUNNING INTO FLOATING POINT CANCELLATION PROBLEMS
112 !-----------------------------------------------------------------------
113 !
114  st_m = abs(dfrac / proth)
115 !-QUADRATIC
116  IF(st_m.GT.1.d-7.AND.abs(dfrac).GT.1.d-12) THEN
117  st_c = f1 / st_m
118  st_a = erosum + 0.5d0 * st_m * st_c**2
119  st_s = sqrt(2.d0 * st_a / st_m)
120  erostrength = st_s - st_c
121 !
122 !-LINEAR (SAVE MODE) FOR ALMOST NO CHANGE IN FRACTION OVER DEPTH
123 ! = ACCEPTAPLE FRACTION SHIFT TO OTHER FRACTIONS (ILLEGAL!!)...
124  ELSE
125  erostrength = erosum / prov * proth
126  ENDIF
127 !
128  erodepth = d1 - erostrength
129 !
130 ! ATTENTION
131  IF((erodepth-d2).LT.0.d0) erodepth = d2
132 !
133 ! IF(CVSP_DB(JG,0)) THEN
134 ! WRITE(LU,*) ' '
135 ! ENDIF
136 !
137 !-----------------------------------------------------------------------
138 ! INSERT A NEW BREAKPOINT TO SPLIT THE VSP AND SHIFT THE REST
139 ! SHIFT
140 !-----------------------------------------------------------------------
141 !
142  DO kk = 0,cntr
143  DO ii=1,nsicla
144  pro_f(j,pro_max(j)-kk+1,ii) = pro_f(j,pro_max(j)-kk,ii)
145  pro_d(j,pro_max(j)-kk+1,ii) = pro_d(j,pro_max(j)-kk,ii)
146  ENDDO
147  ENDDO
148 !
149  DO ii=1,nsicla
150  pro_d(j,pro_max(j)-cntr,ii) = erodepth
151  pro_f(j,pro_max(j)-cntr,ii) = pro_f(j,pro_max(j)-cntr,ii)
152  & - ( pro_f(j,pro_max(j)-cntr,ii)
153  & -pro_f(j,pro_max(j)-(cntr+1),ii))/proth*erostrength
154  ENDDO
155 !
156  IF(pro_d(j,pro_max(j)-cntr,1).LT.
157  & pro_d(j,pro_max(j)-cntr-1,1)) THEN
158  WRITE(lu,*) 'DEPTHINVERSION!!!',jg,pro_max(j)-cntr,
159  & pro_max(j), proth, erostrength
160  CALL plante(1)
161  stop
162  ENDIF
163 !
164  pro_max(j) = pro_max(j) + 1
165 !
166  IF(cvsp_db(jg,0)) CALL cvsp_p('./VSP/','RFT_C3_',jg)
167 !
168  IF(cvsp_check_f(j,pro_max(j)-cntr,'RMF:CAS2 B')) THEN
169  ENDIF
170 !
171 !-----------------------------------------------------------------------
172 ! AFTER SPLITTING: UPDATE SECTION PROPERTIES
173 !-----------------------------------------------------------------------
174 !
175  prof = (pro_f(j,pro_max(j)-cntr,i) + !PROFILE SECTION FRACTION
176  & pro_f(j,pro_max(j)-(cntr+1),i)) * 0.5d0
177 
178  proth = (pro_d(j,pro_max(j)-cntr,i) -
179  & pro_d(j,pro_max(j)-(cntr+1),i)) !PROFILE SECTION THICKNESS
180 
181  prov = prof * proth !PROFILE SECTION VOLUME ( EQ. THICKNESS OF FRACTION)
182 !
183  !DELTA FRACTION = CHANGE OF FRACTION OVER THIS SECTION
184  dfrac = pro_f(j,pro_max(j)-cntr,i) -
185  & pro_f(j,pro_max(j)-(cntr+1),i)
186 
187  erosum = prov !NECESSARY! TO AVOID PROBLEMS WITH SLIVER POLYGONS
188 !
189  ENDIF ! CASE 2 CONVERSION
190 !
191 !-----------------------------------------------------------------------
192 ! CASE 1 (NOW ALL CASES, AS CASE II IS ALREADY CONVERTED)
193 ! THIS VSP SECTION HAS NOT ENOUGH OR EXACTLY ENOUGH MATERIAL TO SATISFY THE HUNGER
194 ! MEANS. IT IS REMOVED COMPLETLY
195 !-----------------------------------------------------------------------
196 !
197  erosum = - prov + erosum
198  ENDIF
199 ! END LOOP THROUGH VSP SECTION
200  END DO
201 !
202 ! DEEPEST/LAST TOUCHED SECTION POINT
203  cntrb = cntr + 1
204 !
205 !-----------------------------------------------------------------------
206 ! FIXING DEPTH AND FRACTION
207 ! AFTER DIGGING OUT 1 FRACTION, EVERYTHING ABOVE FALLS DOWN
208 ! TO FILL THE HOLE
209 !-----------------------------------------------------------------------
210 !
211  IF (cvsp_db(jg,0)) CALL cvsp_p('./VSP/','RFT_C4_',jg)
212 !
213 !-----------------------------------------------------------------------
214 ! SHIFT AND DUPLICATE DEEPEST POINT TO PRODUCE A CLEAR BREAK WHEN NORMALIZING THE FRACTIONS
215 ! IN THE NEXT STEP
216 !-----------------------------------------------------------------------
217 !
218  IF (cntrb.EQ.pro_max(j)) cntrb = pro_max(j) - 1 ! PREVENT RIGID BED DEMOLITION
219  DO kk = 0,cntrb
220 ! DEBUG
221  IF(pro_max(j)-kk+1 > pro_max_max-1) THEN
222  WRITE(lu,*) 'PRO_MAX_MAX_: ',j,pro_max(j)
223  CALL cvsp_p('./','MAX_I',jg)
224  CALL plante(1)
225  ENDIF
226 !
227  DO ii=1,nsicla
228  pro_f(j,pro_max(j)-kk+1,ii) = pro_f(j,pro_max(j)-kk,ii)
229  pro_d(j,pro_max(j)-kk+1,ii) = pro_d(j,pro_max(j)-kk,ii)
230  ENDDO
231  ENDDO
232 !
233 ! 1 POINT MORE / MEANS ONE POINT DEEPER AFFECTED
234  pro_max(j) = pro_max(j) + 1
235 !
236  IF(cvsp_db(jg,0)) CALL cvsp_p('./VSP/','TEST4_',jg)
237 !
238 !-----------------------------------------------------------------------
239 !ALL LAYERS ABOVE FALL DOWN = EVERY FRACTION FALLS, GOING FROM DEEPEST TO HIGHEST AFFECTED SECTION
240 !-----------------------------------------------------------------------
241 !
242  lowpnt = pro_max(j)-cntrb
243 !
244  pro_f_low = pro_f(j,lowpnt,i)
245  pro_d_low = pro_d(j,lowpnt,i) ! START AT LOWESET LEVEL
246  prov_total = 0.d0 ! ACCUMULATING THE EROSION THE HIGHER WE GO
247 !
248  DO kk = lowpnt+1,pro_max(j) ! LOOP OVER UPPER SECTION POINT
249  prof = ( pro_f(j,kk,i) + pro_f_low) * 0.5d0 ! PROFILE SECTION FRACTION
250  proth = ( pro_d(j,kk,i) - pro_d_low) ! PROFILE SECTION THICKNESS
251  prov = prof * proth ! PROFILE SECTION VOLUME (EQ. THICKNESS OF FRACTION)
252 !
253  IF (prov < zero) prov = 0.d0
254  prov_total = prov + prov_total
255  rest1 = ( - pro_f(j,kk,i) + 1.d0) ! SUM OF FRACTIONS AFTER EROSION
256  pro_d_low = pro_d(j,kk,i) ! REMEMBER THIS FOR NEXT SECTION STEP
257  ! CAUSE PRO_D(J,KK,II) IS NOT AVAILABEL ANY MORE
258  pro_f_low = pro_f(j,kk,i) ! KEEP THIS FOR THE NEXT SECTION !
259 !
260  DO ii=1,nsicla
261  pro_d(j,kk,ii) = -prov_total + pro_d(j,kk,ii) ! SHIFT DEPTH
262  IF (pro_d(j,kk,ii).LE.pro_d(j,kk-1,ii)) THEN ! CORRECTING 1D-15 ERRORS WHICH PRODUCE UNSTEADY PDFS
263  pro_d(j,kk,ii) = pro_d(j,kk-1,ii)
264  ENDIF
265  IF (rest1.GT.zero) THEN
266  pro_f(j,kk,ii) = pro_f(j,kk,ii)/ rest1 ! NORMALIZE FRACTION
267  ELSE
268  pro_f(j,kk,ii) = 0.d0 ! IN CASE OF ALMOST TOTAL LOSS
269  ENDIF
270  ENDDO
271 !
272  IF(rest1.GT.zero) THEN
273  pro_f(j,kk,i) = 0.d0 ! FRACTION I IS REMOVED
274  ELSE
275  pro_f(j,kk,i) = 0.d0 ! IN CASE OF ALMOST TOTAL LOSS
276  ENDIF
277 !
278  ENDDO
279 !
280  IF(cvsp_db(jg,0)) CALL cvsp_p('./VSP/','TEST5_',jg)
281 !
282 !-----------------------------------------------------------------------
283 ! SPECIAL TREATMENT LOWEST POINT
284 !-----------------------------------------------------------------------
285 !
286  rest1 = (-pro_f(j,lowpnt,i)+1.d0)
287 !
288  DO ii=1,nsicla
289  IF (rest1.GT.zero) THEN
290  pro_f(j,lowpnt,ii) =
291  & pro_f(j,lowpnt,ii) / rest1 ! NORMALIZE FRACTION
292  ELSE
293  pro_f(j,lowpnt,ii) = 0.d0 ! IN CASE OF ALMOST TOTAL LOSS
294  ENDIF
295  ENDDO
296 !
297  pro_f(j,lowpnt,i) = 0.d0 ! FRACTION I IS REMOVED
298  IF (rest1.LE.zero) pro_f(j,lowpnt,i) =
299  & 0.d0 ! IN CASE OF ALMOST TOTAL LOSS
300 !
301 ! ADDITIONAL CHECK FOR PRO_F
302  DO k = 1, pro_max(j)
303  ret = cvsp_check_f(j,k,'rm_fraction: ')
304  ENDDO
305 !
306  IF(cvsp_db(jg,0)) CALL cvsp_p('./','RFTLOW_',jg)
307  CALL cvsp_compress_clean(j)
308  IF(cvsp_db(jg,0)) CALL cvsp_p('./','RFTCCL_',jg)
309 !
310 !-----------------------------------------------------------------------
311 !
312  RETURN
313  END SUBROUTINE cvsp_rm_fraction
recursive logical function cvsp_check_f(J, K, SOMETEXT)
Definition: cvsp_check_f.f:7
subroutine cvsp_compress_clean(J)
double precision, dimension(:,:,:), allocatable, target pro_f
logical function cvsp_db(J_GLOBAL, TIMESTAMP)
Definition: cvsp_db.f:7
subroutine cvsp_p(PATH_PRE, FILE_PRE, JG)
Definition: cvsp_p.f:7
integer, dimension(:), allocatable pro_max
double precision, dimension(:,:,:), allocatable, target pro_d
subroutine cvsp_rm_fraction(J, I, DZFCL)
type(bief_mesh), target mesh