The TELEMAC-MASCARET system  trunk
cvsp_rm_fraction_gaia.f
Go to the documentation of this file.
1 ! ********************************
2  SUBROUTINE cvsp_rm_fraction_gaia
3 ! ********************************
4 !
5  &(j,i,dzfcl)
6 !
7 !***********************************************************************
8 ! GAIA V8P1 16/05/2017
9 !***********************************************************************
10 !
13 !
17 !
22 !
23 !
28 !
33 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
37 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
38 !
42 !
44  IMPLICIT NONE
45 !
46 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
47 !
48  INTEGER, INTENT(IN) :: J
49  INTEGER, INTENT(IN) :: I
50  DOUBLE PRECISION, INTENT(IN) :: DZFCL
51 !
52 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
53 !
54  DOUBLE PRECISION EROSUM,PROF,PROTH,PROV,REST1,DFRAC
55  DOUBLE PRECISION EROSTRENGTH,ERODEPTH
56  DOUBLE PRECISION PRO_D_LOW,PRO_F_LOW,PROV_TOTAL,F1,F2,D1,D2
57 !
58  INTEGER CNTR, CNTRB, II, KK, LOWPNT, JG, K
59  LOGICAL RET
60  DOUBLE PRECISION ST_S, ST_A, ST_M, ST_C
61  CHARACTER(LEN=7) OCSTR
62 !
63 !-----------------------------------------------------------------------
64 !
65  jg = j
66  IF (ncsize > 1) jg = mesh%KNOLG%I(j)
67 
68  WRITE(unit=ocstr, fmt='(A,I2)') 'RFTIN',i
69 
70  IF (cvsp_db_gaia(jg,0)) CALL cvsp_p_gaia('./',ocstr,jg)
71 !
72  erosum = - dzfcl !EROSION SUM POSITIV = EROSION
73  cntr = -1 !SECTION CNTR BELOW PRO_MAX
74 !
75 !-----------------------------------------------------------------------
76 ! LOOP THROUGH VSP SECTIONS UNTIL ALL EROSION IS DONE
77 ! COUNTING BACK EROSUM -> 0
78 ! IDENTIFY THE LAST SECTION POINT AFFECTED
79 !-----------------------------------------------------------------------
80 !
81  DO WHILE ((erosum > 0.d0.OR.cntr==-1).AND.(cntr.LT.pro_max(j)-2))
82  cntr = cntr + 1
83 !
84  f1 = pro_f(j,pro_max(j)-cntr,i)
85  f2 = pro_f(j,pro_max(j)-(cntr+1),i)
86 !
87  d1 = pro_d(j,pro_max(j)-cntr,i)
88  d2 = pro_d(j,pro_max(j)-(cntr+1),i)
89 !
90  prof = (f1 + f2) * 0.5d0 ! MEAN FRACTION
91  proth = d1 - d2 ! PROFILE SECTION THICKNESS
92  prov = prof * proth ! PROFILE SECTION VOLUME ( EQ. THICKNESS OF FRACTION)
93  dfrac = f1 - f2 ! DELTA FRACTION = CHANGE OF FRACTION OVER THIS SECTION
94 !
95  IF (proth.GT.0.d0) THEN
96 !
97 !-----------------------------------------------------------------------
98 ! CASE 2
99 ! THIS VSP SECTION HAS MORE THEN ENOUGH MATERIAL
100 ! THEN SPLIT IT AT THE DEPTH OF MAXIMUM EROSION,
101 ! TO KEEP ANTHING BELOW UNCHANGED
102 ! => MAKES CASE 2 TO CASE 1
103 !-----------------------------------------------------------------------
104 !
105  IF(erosum < prov) THEN
106  IF(cvsp_db_gaia(jg,0))
107  & CALL cvsp_p_gaia('./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_GAIA(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_gaia(jg,0))
167  & CALL cvsp_p_gaia('./VSP/','RFT_C3_',jg)
168 !
169  IF(cvsp_check_f_gaia(j,pro_max(j)-cntr,'RMF:CAS2 B')) THEN
170  ENDIF
171 !
172 !-----------------------------------------------------------------------
173 ! AFTER SPLITTING: UPDATE SECTION PROPERTIES
174 !-----------------------------------------------------------------------
175 !
176  prof = (pro_f(j,pro_max(j)-cntr,i) + !PROFILE SECTION FRACTION
177  & pro_f(j,pro_max(j)-(cntr+1),i)) * 0.5d0
178 
179  proth = (pro_d(j,pro_max(j)-cntr,i) -
180  & pro_d(j,pro_max(j)-(cntr+1),i)) !PROFILE SECTION THICKNESS
181 
182  prov = prof * proth !PROFILE SECTION VOLUME ( EQ. THICKNESS OF FRACTION)
183 !
184  !DELTA FRACTION = CHANGE OF FRACTION OVER THIS SECTION
185  dfrac = pro_f(j,pro_max(j)-cntr,i) -
186  & pro_f(j,pro_max(j)-(cntr+1),i)
187 
188  erosum = prov !NECESSARY! TO AVOID PROBLEMS WITH SLIVER POLYGONS
189 !
190  ENDIF ! CASE 2 CONVERSION
191 !
192 !-----------------------------------------------------------------------
193 ! CASE 1 (NOW ALL CASES, AS CASE II IS ALREADY CONVERTED)
194 ! THIS VSP SECTION HAS NOT ENOUGH OR EXACTLY ENOUGH MATERIAL TO SATISFY THE HUNGER
195 ! MEANS. IT IS REMOVED COMPLETLY
196 !-----------------------------------------------------------------------
197 !
198  erosum = - prov + erosum
199  ENDIF
200 ! END LOOP THROUGH VSP SECTION
201  END DO
202 !
203 ! DEEPEST/LAST TOUCHED SECTION POINT
204  cntrb = cntr + 1
205 !
206 !-----------------------------------------------------------------------
207 ! FIXING DEPTH AND FRACTION
208 ! AFTER DIGGING OUT 1 FRACTION, EVERYTHING ABOVE FALLS DOWN
209 ! TO FILL THE HOLE
210 !-----------------------------------------------------------------------
211 !
212  IF (cvsp_db_gaia(jg,0)) CALL cvsp_p_gaia('./VSP/','RFT_C4_',jg)
213 !
214 !-----------------------------------------------------------------------
215 ! SHIFT AND DUPLICATE DEEPEST POINT TO PRODUCE A CLEAR BREAK WHEN NORMALIZING THE FRACTIONS
216 ! IN THE NEXT STEP
217 !-----------------------------------------------------------------------
218 !
219  IF (cntrb.EQ.pro_max(j)) cntrb = pro_max(j) - 1 ! PREVENT RIGID BED DEMOLITION
220  DO kk = 0,cntrb
221 ! DEBUG
222  IF(pro_max(j)-kk+1 > pro_max_max-1) THEN
223  WRITE(lu,*) 'PRO_MAX_MAX_: ',j,pro_max(j)
224  CALL cvsp_p_gaia('./','MAX_I',jg)
225  CALL plante(1)
226  ENDIF
227 !
228  DO ii=1,nsicla
229  pro_f(j,pro_max(j)-kk+1,ii) = pro_f(j,pro_max(j)-kk,ii)
230  pro_d(j,pro_max(j)-kk+1,ii) = pro_d(j,pro_max(j)-kk,ii)
231  ENDDO
232  ENDDO
233 !
234 ! 1 POINT MORE / MEANS ONE POINT DEEPER AFFECTED
235  pro_max(j) = pro_max(j) + 1
236 !
237  IF(cvsp_db_gaia(jg,0)) CALL cvsp_p_gaia('./VSP/','TEST4_',jg)
238 !
239 !-----------------------------------------------------------------------
240 !ALL LAYERS ABOVE FALL DOWN = EVERY FRACTION FALLS, GOING FROM DEEPEST TO HIGHEST AFFECTED SECTION
241 !-----------------------------------------------------------------------
242 !
243  lowpnt = pro_max(j)-cntrb
244 !
245  pro_f_low = pro_f(j,lowpnt,i)
246  pro_d_low = pro_d(j,lowpnt,i) ! START AT LOWESET LEVEL
247  prov_total = 0.d0 ! ACCUMULATING THE EROSION THE HIGHER WE GO
248 !
249  DO kk = lowpnt+1,pro_max(j) ! LOOP OVER UPPER SECTION POINT
250  prof = ( pro_f(j,kk,i) + pro_f_low) * 0.5d0 ! PROFILE SECTION FRACTION
251  proth = ( pro_d(j,kk,i) - pro_d_low) ! PROFILE SECTION THICKNESS
252  prov = prof * proth ! PROFILE SECTION VOLUME (EQ. THICKNESS OF FRACTION)
253 !
254  IF (prov < zero) prov = 0.d0
255  prov_total = prov + prov_total
256  rest1 = ( - pro_f(j,kk,i) + 1.d0) ! SUM OF FRACTIONS AFTER EROSION
257  pro_d_low = pro_d(j,kk,i) ! REMEMBER THIS FOR NEXT SECTION STEP
258  ! CAUSE PRO_D(J,KK,II) IS NOT AVAILABEL ANY MORE
259  pro_f_low = pro_f(j,kk,i) ! KEEP THIS FOR THE NEXT SECTION !
260 !
261  DO ii=1,nsicla
262  pro_d(j,kk,ii) = -prov_total + pro_d(j,kk,ii) ! SHIFT DEPTH
263  IF (pro_d(j,kk,ii).LE.pro_d(j,kk-1,ii)) THEN ! CORRECTING 1D-15 ERRORS WHICH PRODUCE UNSTEADY PDFS
264  pro_d(j,kk,ii) = pro_d(j,kk-1,ii)
265  ENDIF
266  IF (rest1.GT.zero) THEN
267  pro_f(j,kk,ii) = pro_f(j,kk,ii)/ rest1 ! NORMALIZE FRACTION
268  ELSE
269  pro_f(j,kk,ii) = 0.d0 ! IN CASE OF ALMOST TOTAL LOSS
270  ENDIF
271  ENDDO
272 !
273  IF(rest1.GT.zero) THEN
274  pro_f(j,kk,i) = 0.d0 ! FRACTION I IS REMOVED
275  ELSE
276  pro_f(j,kk,i) = 0.d0 ! IN CASE OF ALMOST TOTAL LOSS
277  ENDIF
278 !
279  ENDDO
280 !
281  IF(cvsp_db_gaia(jg,0)) CALL cvsp_p_gaia('./VSP/','TEST5_',jg)
282 !
283 !-----------------------------------------------------------------------
284 ! SPECIAL TREATMENT LOWEST POINT
285 !-----------------------------------------------------------------------
286 !
287  rest1 = (-pro_f(j,lowpnt,i)+1.d0)
288 !
289  DO ii=1,nsicla
290  IF (rest1.GT.zero) THEN
291  pro_f(j,lowpnt,ii) =
292  & pro_f(j,lowpnt,ii) / rest1 ! NORMALIZE FRACTION
293  ELSE
294  pro_f(j,lowpnt,ii) = 0.d0 ! IN CASE OF ALMOST TOTAL LOSS
295  ENDIF
296  ENDDO
297 !
298  pro_f(j,lowpnt,i) = 0.d0 ! FRACTION I IS REMOVED
299  IF (rest1.LE.zero) pro_f(j,lowpnt,i) =
300  & 0.d0 ! IN CASE OF ALMOST TOTAL LOSS
301 !
302 ! ADDITIONAL CHECK FOR PRO_F
303  DO k = 1, pro_max(j)
304  ret = cvsp_check_f_gaia(j,k,'rm_fraction: ')
305  ENDDO
306 !
307  IF(cvsp_db_gaia(jg,0)) CALL cvsp_p_gaia('./','RFTLOW_',jg)
309  IF(cvsp_db_gaia(jg,0)) CALL cvsp_p_gaia('./','RFTCCL_',jg)
310 !
311 !-----------------------------------------------------------------------
312 !
313  RETURN
314  END SUBROUTINE cvsp_rm_fraction_gaia
subroutine cvsp_compress_clean_gaia(J)
subroutine cvsp_p_gaia(PATH_PRE, FILE_PRE, JG)
Definition: cvsp_p_gaia.f:7
subroutine cvsp_rm_fraction_gaia(J, I, DZFCL)
double precision zero
Parameter used for clipping variables or testing values against zero.
integer, target nsicla
Number of sediment classes of bed material (less than NISCLM)
double precision, dimension(:,:,:), allocatable, target pro_f
Vertical sorting profile: fraction for each layer, class, point.
double precision, dimension(:,:,:), allocatable, target pro_d
Vertical sorting profile: depth for each layer, class, point.
recursive logical function cvsp_check_f_gaia(J, K, SOMETEXT)
integer pro_max_max
Maximum Number of Profile SECTIONS.
logical function cvsp_db_gaia(J_GLOBAL, TIMESTAMP)
Definition: cvsp_db_gaia.f:7
type(bief_mesh), target mesh
Mesh structure.
integer, dimension(:), allocatable pro_max
Maximum layer number in a vertical sorting profile for each point.