cvsp_rm_fraction.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\sisyphe\cvsp_rm_fraction.f
00002 !
00065                      SUBROUTINE CVSP_RM_FRACTION
00066 !                    ***************************
00067 !
00068      &(J,I,DZFCL,EVL)
00069 !
00070 !***********************************************************************
00071 ! SISYPHE   V6P3                                   14/03/2013
00072 !***********************************************************************
00073 !
00074 !
00075 !
00076 !
00077 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00078 !| J              |<--| INDEX OF A POINT IN MESH
00079 !| I              |<--| INDEX OF A FRACTION
00080 !| DZFCL          |<--| VALUE OF A FRACTION IN CM !
00081 !| EVL            |<--| SUM OF ALL FRACTIONS
00082 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00083 !
00084       USE DECLARATIONS_SISYPHE
00085       USE CVSP_OUTPUTFILES
00086 !
00087       IMPLICIT NONE
00088       INTEGER LNG,LU
00089       COMMON/INFO/LNG,LU
00090 !
00091 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00092 !
00093       INTEGER,          INTENT(IN)    :: J
00094       INTEGER,          INTENT(IN)    :: I
00095       DOUBLE PRECISION, INTENT(IN)    :: DZFCL
00096       DOUBLE PRECISION, INTENT(IN)    :: EVL
00097 !
00098 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00099 !
00100       LOGICAL, EXTERNAL :: CVSP_CHECK_F, DB
00101 !
00102       DOUBLE PRECISION EROSUM,PROF,PROTH,PROV,REST1,DFRAC
00103       DOUBLE PRECISION EROSTRENGTH,ERODEPTH
00104       DOUBLE PRECISION PRO_D_LOW,PRO_F_LOW,PROV_TOTAL,F1,F2,D1,D2
00105 !
00106       INTEGER CNTR, CNTRB, II, KK,  LOWPNT, JG
00107       DOUBLE PRECISION  ST_S, ST_A, ST_M, ST_C
00108 !
00109 !-----------------------------------------------------------------------
00110 !
00111       JG = J
00112       IF (NCSIZE > 1) JG = MESH%KNOLG%I(J)
00113       IF (DB(JG,0)) CALL CVSP_P('./VSP/','RFT_IN_',JG)
00114 !
00115       EROSUM = - DZFCL    !EROSION SUM POSITIV = EROSION
00116       CNTR = -1                 !SECTION CNTR BELOW PRO_MAX
00117 !
00118 !-----------------------------------------------------------------------
00119 ! LOOP THROUGH VSP SECTIONS UNTIL ALL EROSION IS DONE
00120 ! COUNTING BACK EROSUM -> 0
00121 ! IDENTIFY THE LAST SECTION POINT AFFECTED
00122 !-----------------------------------------------------------------------
00123 !
00124       DO WHILE ((EROSUM > 0.D0.OR.CNTR==-1).AND.(CNTR.LT.PRO_MAX(J)-2))
00125         CNTR  =  CNTR + 1
00126 !
00127         F1 = PRO_F(J,PRO_MAX(J)-CNTR,I)
00128         F2 = PRO_F(J,PRO_MAX(J)-(CNTR+1),I)
00129 !
00130         D1 = PRO_D(J,PRO_MAX(J)-CNTR,I)
00131         D2 = PRO_D(J,PRO_MAX(J)-(CNTR+1),I)
00132 !
00133         PROF  = (F1 + F2) * 0.5D0 ! MEAN FRACTION
00134         PROTH = D1 - D2           ! PROFILE SECTION THICKNESS
00135         PROV  =  PROF * PROTH     ! PROFILE SECTION VOLUME ( EQ. THICKNESS OF FRACTION)
00136         DFRAC =  F1 - F2          ! DELTA FRACTION = CHANGE OF FRACTION OVER THIS SECTION
00137 !
00138         IF (PROTH.GT.0.D0) THEN
00139 !
00140 !-----------------------------------------------------------------------
00141 ! CASE 2
00142 ! THIS VSP SECTION HAS MORE THEN ENOUGH MATERIAL
00143 ! THEN SPLIT IT AT THE DEPTH OF MAXIMUM EROSION,
00144 ! TO KEEP ANTHING BELOW UNCHANGED
00145 ! => MAKES CASE 2 TO CASE 1
00146 !-----------------------------------------------------------------------
00147 !
00148           IF(EROSUM < PROV) THEN
00149             IF(DB(JG,0)) CALL CVSP_P('./VSP/','RFT_C2_',JG)
00150 !
00151 !-----------------------------------------------------------------------
00152 ! MAXIMUM DEPTH OF EROSION: SOLVING QUADRATIC PROBLEM BY CROSS-MULTIPLICATION
00153 ! LINEAR OR QUATRATIC ? PREVENTS RUNNING INTO FLOATING POINT CANCELLATION PROBLEMS
00154 !-----------------------------------------------------------------------
00155 !
00156             ST_M = ABS(DFRAC / PROTH)
00157 !-QUADRATIC
00158             IF(ST_M.GT.1.D-7.AND.ABS(DFRAC).GT.1.D-12) THEN
00159               ST_C = F1 / ST_M
00160               ST_A = EROSUM + 0.5D0 * ST_M * ST_C**2
00161               ST_S = SQRT(2.D0 * ST_A / ST_M)
00162               EROSTRENGTH = ST_S - ST_C
00163 !
00164 !-LINEAR (SAVE MODE) FOR ALMOST NO CHANGE IN FRACTION OVER DEPTH
00165 ! = ACCEPTAPLE FRACTION SHIFT TO OTHER FRACTIONS (ILLEGAL!!)...
00166             ELSE
00167               EROSTRENGTH = EROSUM / PROV * PROTH
00168             ENDIF
00169 !
00170             ERODEPTH = D1 - EROSTRENGTH
00171 !
00172 !     ATTENTION
00173             IF((ERODEPTH-D2).LT.0.D0) ERODEPTH = D2
00174 !
00175             IF(DB(JG,0)) THEN
00176 !             PRINT*, '  '
00177             ENDIF
00178 !
00179 !-----------------------------------------------------------------------
00180 ! INSERT A NEW BREAKPOINT TO SPLIT THE VSP AND SHIFT THE REST
00181 ! SHIFT
00182 !-----------------------------------------------------------------------
00183 !
00184             DO KK = 0,CNTR
00185               DO II=1,NSICLA
00186                 PRO_F(J,PRO_MAX(J)-KK+1,II) = PRO_F(J,PRO_MAX(J)-KK,II)
00187                 PRO_D(J,PRO_MAX(J)-KK+1,II) = PRO_D(J,PRO_MAX(J)-KK,II)
00188               ENDDO
00189             ENDDO
00190 !
00191             DO II=1,NSICLA
00192               PRO_D(J,PRO_MAX(J)-CNTR,II) = ERODEPTH
00193               PRO_F(J,PRO_MAX(J)-CNTR,II) = PRO_F(J,PRO_MAX(J)-CNTR,II)
00194      &             - ( PRO_F(J,PRO_MAX(J)-CNTR,II)
00195      &             -PRO_F(J,PRO_MAX(J)-(CNTR+1),II))/PROTH*EROSTRENGTH
00196             ENDDO
00197 !
00198             IF(PRO_D(J,PRO_MAX(J)-CNTR,1).LT.
00199      &         PRO_D(J,PRO_MAX(J)-CNTR-1,1)) THEN
00200               WRITE(LU,*) 'DEPTHINVERSION!!!',JG,PRO_MAX(J)-CNTR,
00201      &             PRO_MAX(J), PROTH, EROSTRENGTH
00202               CALL PLANTE(1)
00203               STOP
00204             ENDIF
00205 !
00206             PRO_MAX(J) = PRO_MAX(J) + 1
00207 !
00208             IF(DB(JG,0)) CALL CVSP_P('./VSP/','RFT_C3_',JG)
00209 !
00210             IF(CVSP_CHECK_F(J,PRO_MAX(J)-CNTR,'RMF:CAS2 B')) THEN
00211             ENDIF
00212 !
00213 !-----------------------------------------------------------------------
00214 ! AFTER SPLITTING: UPDATE SECTION PROPERTIES
00215 !-----------------------------------------------------------------------
00216 !
00217             PROF  = (PRO_F(J,PRO_MAX(J)-CNTR,I) + !PROFILE SECTION FRACTION
00218      &           PRO_F(J,PRO_MAX(J)-(CNTR+1),I)) * 0.5D0
00219 
00220             PROTH = (PRO_D(J,PRO_MAX(J)-CNTR,I) -
00221      &           PRO_D(J,PRO_MAX(J)-(CNTR+1),I)) !PROFILE SECTION THICKNESS
00222 
00223             PROV  =  PROF * PROTH    !PROFILE SECTION VOLUME ( EQ. THICKNESS OF FRACTION)
00224 !
00225             !DELTA FRACTION = CHANGE OF FRACTION OVER THIS SECTION
00226             DFRAC = PRO_F(J,PRO_MAX(J)-CNTR,I) -
00227      &              PRO_F(J,PRO_MAX(J)-(CNTR+1),I)
00228 
00229             EROSUM = PROV            !NECESSARY! TO AVOID PROBLEMS WITH SLIVER POLYGONS
00230 !
00231           ENDIF                     ! CASE 2 CONVERSION
00232 !
00233 !-----------------------------------------------------------------------
00234 ! CASE 1 (NOW ALL CASES, AS CASE II IS ALREADY CONVERTED)
00235 ! THIS VSP SECTION HAS NOT ENOUGH OR EXACTLY ENOUGH MATERIAL TO SATISFY THE HUNGER
00236 ! MEANS. IT IS REMOVED COMPLETLY
00237 !-----------------------------------------------------------------------
00238 !
00239         EROSUM = - PROV + EROSUM
00240         ENDIF
00241 !     END LOOP THROUGH VSP SECTION
00242       END DO
00243 !
00244 !     DEEPEST/LAST TOUCHED SECTION POINT
00245       CNTRB = CNTR + 1
00246 !
00247 !-----------------------------------------------------------------------
00248 ! FIXING DEPTH AND FRACTION
00249 ! AFTER DIGGING OUT 1 FRACTION, EVERYTHING ABOVE FALLS DOWN
00250 ! TO FILL THE HOLE
00251 !-----------------------------------------------------------------------
00252 !
00253       IF (DB(JG,0)) CALL CVSP_P('./VSP/','RFT_C4_',JG)
00254 !
00255 !-----------------------------------------------------------------------
00256 ! SHIFT AND DUPLICATE DEEPEST POINT TO PRODUCE A CLEAR BREAK WHEN NORMALIZING THE FRACTIONS
00257 ! IN THE NEXT STEP
00258 !-----------------------------------------------------------------------
00259 !
00260       IF (CNTRB.EQ.PRO_MAX(J)) CNTRB = PRO_MAX(J) - 1 ! PREVENT RIGID BED DEMOLITION
00261       DO KK = 0,CNTRB
00262 !       DEBUG
00263         IF(PRO_MAX(J)-KK+1 > PRO_MAX_MAX-1) THEN
00264           WRITE(LU,*) 'PRO_MAX_MAX_: ',J,PRO_MAX(J)
00265           CALL CVSP_P('./ERR/','MAX_I',JG)
00266           CALL PLANTE(1)
00267         ENDIF
00268 !
00269         DO II=1,NSICLA
00270           PRO_F(J,PRO_MAX(J)-KK+1,II) = PRO_F(J,PRO_MAX(J)-KK,II)
00271           PRO_D(J,PRO_MAX(J)-KK+1,II) = PRO_D(J,PRO_MAX(J)-KK,II)
00272         ENDDO
00273       ENDDO
00274 !
00275 !     1 POINT MORE / MEANS ONE POINT DEEPER AFFECTED
00276       PRO_MAX(J) = PRO_MAX(J) + 1
00277 !
00278       IF(DB(JG,0)) CALL CVSP_P('./VSP/','TEST4_',JG)
00279 !
00280 !-----------------------------------------------------------------------
00281 !ALL LAYERS ABOVE FALL DOWN = EVERY FRACTION FALLS, GOING FROM DEEPEST TO HIGHEST AFFECTED SECTION
00282 !-----------------------------------------------------------------------
00283 !
00284       LOWPNT = PRO_MAX(J)-CNTRB
00285 !
00286       PRO_F_LOW = PRO_F(J,LOWPNT,I)
00287       PRO_D_LOW = PRO_D(J,LOWPNT,I) ! START AT LOWESET LEVEL
00288       PROV_TOTAL = 0.D0             ! ACCUMULATING THE EROSION THE HIGHER WE GO
00289 !
00290       DO KK = LOWPNT+1,PRO_MAX(J)                     ! LOOP OVER UPPER SECTION POINT
00291         PROF  = ( PRO_F(J,KK,I) + PRO_F_LOW) * 0.5D0 ! PROFILE SECTION FRACTION
00292         PROTH = ( PRO_D(J,KK,I) - PRO_D_LOW)         ! PROFILE SECTION THICKNESS
00293         PROV  =  PROF * PROTH                        ! PROFILE SECTION VOLUME (EQ. THICKNESS OF FRACTION)
00294 !
00295         IF (PROV < ZERO) PROV = 0.D0
00296         PROV_TOTAL = PROV + PROV_TOTAL
00297         REST1 = ( - PRO_F(J,KK,I) + 1.D0) ! SUM OF FRACTIONS AFTER EROSION
00298         PRO_D_LOW = PRO_D(J,KK,I)         ! REMEMBER THIS FOR NEXT SECTION STEP
00299                                           ! CAUSE PRO_D(J,KK,II) IS NOT AVAILABEL ANY MORE
00300         PRO_F_LOW = PRO_F(J,KK,I)         ! KEEP THIS FOR THE NEXT SECTION !
00301 !
00302         DO II=1,NSICLA
00303           PRO_D(J,KK,II) = -PROV_TOTAL + PRO_D(J,KK,II) ! SHIFT DEPTH
00304           IF (PRO_D(J,KK,II).LE.PRO_D(J,KK-1,II)) THEN  ! CORRECTING 1D-15 ERRORS WHICH PRODUCE UNSTEADY PDFS
00305             PRO_D(J,KK,II) = PRO_D(J,KK-1,II)
00306           ENDIF
00307           IF (REST1.GT.ZERO) THEN
00308             PRO_F(J,KK,II) = PRO_F(J,KK,II)/ REST1 ! NORMALIZE FRACTION
00309           ELSE
00310             PRO_F(J,KK,II) = 1.D0 / NSICLA         ! IN CASE OF ALMOST TOTAL LOSS
00311           ENDIF
00312         ENDDO
00313 !
00314         IF(REST1.GT.ZERO) THEN
00315           PRO_F(J,KK,I) = 0.D0          ! FRACTION I IS REMOVED
00316         ELSE
00317           PRO_F(J,KK,I) = 1.D0 / NSICLA ! IN CASE OF ALMOST TOTAL LOSS
00318         ENDIF
00319 !
00320       ENDDO
00321 !
00322       IF(DB(JG,0)) CALL CVSP_P('./VSP/','TEST5_',JG)
00323 !
00324 !-----------------------------------------------------------------------
00325 ! SPECIAL TREATMENT LOWEST POINT
00326 !-----------------------------------------------------------------------
00327 !
00328       REST1 = (-PRO_F(J,LOWPNT,I)+1.D0)
00329 !
00330       DO II=1,NSICLA
00331         IF (REST1.GT.ZERO) THEN
00332           PRO_F(J,LOWPNT,II) =
00333      &         PRO_F(J,LOWPNT,II) / REST1    ! NORMALIZE FRACTION
00334         ELSE
00335           PRO_F(J,LOWPNT,II) = 1.D0 / NSICLA ! IN CASE OF ALMOST TOTAL LOSS
00336         ENDIF
00337       ENDDO
00338 !
00339       PRO_F(J,LOWPNT,I) = 0.D0                 ! FRACTION I IS REMOVED
00340       IF (REST1.LE.ZERO) PRO_F(J,LOWPNT,I) =
00341      &     1.D0 / NSICLA                       ! IN CASE OF ALMOST TOTAL LOSS
00342 !
00343       IF(DB(JG,0)) CALL CVSP_P('./VSP/','RFT_OU_',JG)
00344 !
00345       CALL CVSP_COMPRESS_CLEAN(J)
00346 !
00347 !-----------------------------------------------------------------------
00348 !
00349       RETURN
00350       END SUBROUTINE CVSP_RM_FRACTION

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0