cvsp_rm_fraction.f

Go to the documentation of this file.

!                   ***************************
                    SUBROUTINE cvsp_rm_fraction
!                   ***************************
!
     &(j,i,dzfcl)
!
!***********************************************************************
! SISYPHE   V7P2                                   16/05/2017
!***********************************************************************
!
!brief   REMOVES (PARTS) OF A FRACTION AFTER EROSION
!+       FROM THE VERTICAL SORTING PROFILE;
!
!history UWE MERKEL
!+        02/02/2012
!+        V6P2
!+
!
!history  P. A. TASSI (EDF R&D, LNHE)
!+        12/03/2013
!+        V6P3
!+   Cleaning, cosmetic
!
!
!history  U. MERKEL (UHM), R. KOPMANN (BAW)
!+        12/06/2016 / 2017
!+        V6P3 / V7P2
!+   Update
!
!history  R. KOPMANN (BAW)
!+        25/02/2019
!+        V7P2
!+   Removing 1/NSICLA
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| J              |<--| INDEX OF A POINT IN MESH
!| I              |<--| INDEX OF A FRACTION
!| DZFCL          |<--| VALUE OF A FRACTION IN CM !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe
      USE declarations_sisyphe
      USE cvsp_outputfiles
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER,          INTENT(IN)    :: J
      INTEGER,          INTENT(IN)    :: I
      DOUBLE PRECISION, INTENT(IN)    :: DZFCL
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION EROSUM,PROF,PROTH,PROV,REST1,DFRAC
      DOUBLE PRECISION EROSTRENGTH,ERODEPTH
      DOUBLE PRECISION PRO_D_LOW,PRO_F_LOW,PROV_TOTAL,F1,F2,D1,D2
!
      INTEGER CNTR, CNTRB, II, KK,  LOWPNT, JG, K
      LOGICAL RET
      DOUBLE PRECISION  ST_S, ST_A, ST_M, ST_C
      CHARACTER(LEN=7) OCSTR
!
!-----------------------------------------------------------------------
!
      jg = j
      IF (ncsize > 1) jg = mesh%KNOLG%I(j)

        WRITE(unit=ocstr, fmt='(A,I2)') 'RFTIN',i

      IF (cvsp_db(jg,0)) CALL cvsp_p('./',ocstr,jg)
!
      erosum = - dzfcl    !EROSION SUM POSITIV = EROSION
      cntr = -1                 !SECTION CNTR BELOW PRO_MAX
!
!-----------------------------------------------------------------------
! LOOP THROUGH VSP SECTIONS UNTIL ALL EROSION IS DONE
! COUNTING BACK EROSUM -> 0
! IDENTIFY THE LAST SECTION POINT AFFECTED
!-----------------------------------------------------------------------
!
      DO WHILE ((erosum > 0.d0.OR.cntr==-1).AND.(cntr.LT.pro_max(j)-2))
        cntr  =  cntr + 1
!
        f1 = pro_f(j,pro_max(j)-cntr,i)
        f2 = pro_f(j,pro_max(j)-(cntr+1),i)
!
        d1 = pro_d(j,pro_max(j)-cntr,i)
        d2 = pro_d(j,pro_max(j)-(cntr+1),i)
!
        prof  = (f1 + f2) * 0.5d0 ! MEAN FRACTION
        proth = d1 - d2           ! PROFILE SECTION THICKNESS
        prov  =  prof * proth     ! PROFILE SECTION VOLUME ( EQ. THICKNESS OF FRACTION)
        dfrac =  f1 - f2          ! DELTA FRACTION = CHANGE OF FRACTION OVER THIS SECTION
!
        IF (proth.GT.0.d0) THEN
!
!-----------------------------------------------------------------------
! CASE 2
! THIS VSP SECTION HAS MORE THEN ENOUGH MATERIAL
! THEN SPLIT IT AT THE DEPTH OF MAXIMUM EROSION,
! TO KEEP ANTHING BELOW UNCHANGED
! => MAKES CASE 2 TO CASE 1
!-----------------------------------------------------------------------
!
          IF(erosum < prov) THEN
            IF(cvsp_db(jg,0)) CALL cvsp_p('./VSP/','RFT_C2_',jg)
!
!-----------------------------------------------------------------------
! MAXIMUM DEPTH OF EROSION: SOLVING QUADRATIC PROBLEM BY CROSS-MULTIPLICATION
! LINEAR OR QUATRATIC ? PREVENTS RUNNING INTO FLOATING POINT CANCELLATION PROBLEMS
!-----------------------------------------------------------------------
!
            st_m = abs(dfrac / proth)
!-QUADRATIC
            IF(st_m.GT.1.d-7.AND.abs(dfrac).GT.1.d-12) THEN
              st_c = f1 / st_m
              st_a = erosum + 0.5d0 * st_m * st_c**2
              st_s = sqrt(2.d0 * st_a / st_m)
              erostrength = st_s - st_c
!
!-LINEAR (SAVE MODE) FOR ALMOST NO CHANGE IN FRACTION OVER DEPTH
! = ACCEPTAPLE FRACTION SHIFT TO OTHER FRACTIONS (ILLEGAL!!)...
            ELSE
              erostrength = erosum / prov * proth
            ENDIF
!
            erodepth = d1 - erostrength
!
!     ATTENTION
            IF((erodepth-d2).LT.0.d0) erodepth = d2
!
!           IF(CVSP_DB(JG,0)) THEN
!             WRITE(LU,*) '  '
!           ENDIF
!
!-----------------------------------------------------------------------
! INSERT A NEW BREAKPOINT TO SPLIT THE VSP AND SHIFT THE REST
! SHIFT
!-----------------------------------------------------------------------
!
            DO kk = 0,cntr
              DO ii=1,nsicla
                pro_f(j,pro_max(j)-kk+1,ii) = pro_f(j,pro_max(j)-kk,ii)
                pro_d(j,pro_max(j)-kk+1,ii) = pro_d(j,pro_max(j)-kk,ii)
              ENDDO
            ENDDO
!
            DO ii=1,nsicla
              pro_d(j,pro_max(j)-cntr,ii) = erodepth
              pro_f(j,pro_max(j)-cntr,ii) = pro_f(j,pro_max(j)-cntr,ii)
     &             - ( pro_f(j,pro_max(j)-cntr,ii)
     &             -pro_f(j,pro_max(j)-(cntr+1),ii))/proth*erostrength
            ENDDO
!
            IF(pro_d(j,pro_max(j)-cntr,1).LT.
     &         pro_d(j,pro_max(j)-cntr-1,1)) THEN
              WRITE(lu,*) 'DEPTHINVERSION!!!',jg,pro_max(j)-cntr,
     &             pro_max(j), proth, erostrength
              CALL plante(1)
              stop
            ENDIF
!
            pro_max(j) = pro_max(j) + 1
!
            IF(cvsp_db(jg,0)) CALL cvsp_p('./VSP/','RFT_C3_',jg)
!
            IF(cvsp_check_f(j,pro_max(j)-cntr,'RMF:CAS2 B')) THEN
            ENDIF
!
!-----------------------------------------------------------------------
! AFTER SPLITTING: UPDATE SECTION PROPERTIES
!-----------------------------------------------------------------------
!
            prof  = (pro_f(j,pro_max(j)-cntr,i) + !PROFILE SECTION FRACTION
     &           pro_f(j,pro_max(j)-(cntr+1),i)) * 0.5d0

            proth = (pro_d(j,pro_max(j)-cntr,i) -
     &           pro_d(j,pro_max(j)-(cntr+1),i)) !PROFILE SECTION THICKNESS

            prov  =  prof * proth    !PROFILE SECTION VOLUME ( EQ. THICKNESS OF FRACTION)
!
            !DELTA FRACTION = CHANGE OF FRACTION OVER THIS SECTION
            dfrac = pro_f(j,pro_max(j)-cntr,i) -
     &              pro_f(j,pro_max(j)-(cntr+1),i)

            erosum = prov            !NECESSARY! TO AVOID PROBLEMS WITH SLIVER POLYGONS
!
          ENDIF                     ! CASE 2 CONVERSION
!
!-----------------------------------------------------------------------
! CASE 1 (NOW ALL CASES, AS CASE II IS ALREADY CONVERTED)
! THIS VSP SECTION HAS NOT ENOUGH OR EXACTLY ENOUGH MATERIAL TO SATISFY THE HUNGER
! MEANS. IT IS REMOVED COMPLETLY
!-----------------------------------------------------------------------
!
        erosum = - prov + erosum
        ENDIF
!     END LOOP THROUGH VSP SECTION
      END DO
!
!     DEEPEST/LAST TOUCHED SECTION POINT
      cntrb = cntr + 1
!
!-----------------------------------------------------------------------
! FIXING DEPTH AND FRACTION
! AFTER DIGGING OUT 1 FRACTION, EVERYTHING ABOVE FALLS DOWN
! TO FILL THE HOLE
!-----------------------------------------------------------------------
!
      IF (cvsp_db(jg,0)) CALL cvsp_p('./VSP/','RFT_C4_',jg)
!
!-----------------------------------------------------------------------
! SHIFT AND DUPLICATE DEEPEST POINT TO PRODUCE A CLEAR BREAK WHEN NORMALIZING THE FRACTIONS
! IN THE NEXT STEP
!-----------------------------------------------------------------------
!
      IF (cntrb.EQ.pro_max(j)) cntrb = pro_max(j) - 1 ! PREVENT RIGID BED DEMOLITION
      DO kk = 0,cntrb
!       DEBUG
        IF(pro_max(j)-kk+1 > pro_max_max-1) THEN
          WRITE(lu,*) 'PRO_MAX_MAX_: ',j,pro_max(j)
          CALL cvsp_p('./','MAX_I',jg)
          CALL plante(1)
        ENDIF
!
        DO ii=1,nsicla
          pro_f(j,pro_max(j)-kk+1,ii) = pro_f(j,pro_max(j)-kk,ii)
          pro_d(j,pro_max(j)-kk+1,ii) = pro_d(j,pro_max(j)-kk,ii)
        ENDDO
      ENDDO
!
!     1 POINT MORE / MEANS ONE POINT DEEPER AFFECTED
      pro_max(j) = pro_max(j) + 1
!
      IF(cvsp_db(jg,0)) CALL cvsp_p('./VSP/','TEST4_',jg)
!
!-----------------------------------------------------------------------
!ALL LAYERS ABOVE FALL DOWN = EVERY FRACTION FALLS, GOING FROM DEEPEST TO HIGHEST AFFECTED SECTION
!-----------------------------------------------------------------------
!
      lowpnt = pro_max(j)-cntrb
!
      pro_f_low = pro_f(j,lowpnt,i)
      pro_d_low = pro_d(j,lowpnt,i) ! START AT LOWESET LEVEL
      prov_total = 0.d0             ! ACCUMULATING THE EROSION THE HIGHER WE GO
!
      DO kk = lowpnt+1,pro_max(j)                     ! LOOP OVER UPPER SECTION POINT
        prof  = ( pro_f(j,kk,i) + pro_f_low) * 0.5d0 ! PROFILE SECTION FRACTION
        proth = ( pro_d(j,kk,i) - pro_d_low)         ! PROFILE SECTION THICKNESS
        prov  =  prof * proth                        ! PROFILE SECTION VOLUME (EQ. THICKNESS OF FRACTION)
!
        IF (prov < zero) prov = 0.d0
        prov_total = prov + prov_total
        rest1 = ( - pro_f(j,kk,i) + 1.d0) ! SUM OF FRACTIONS AFTER EROSION
        pro_d_low = pro_d(j,kk,i)         ! REMEMBER THIS FOR NEXT SECTION STEP
                                          ! CAUSE PRO_D(J,KK,II) IS NOT AVAILABEL ANY MORE
        pro_f_low = pro_f(j,kk,i)         ! KEEP THIS FOR THE NEXT SECTION !
!
        DO ii=1,nsicla
          pro_d(j,kk,ii) = -prov_total + pro_d(j,kk,ii) ! SHIFT DEPTH
          IF (pro_d(j,kk,ii).LE.pro_d(j,kk-1,ii)) THEN  ! CORRECTING 1D-15 ERRORS WHICH PRODUCE UNSTEADY PDFS
            pro_d(j,kk,ii) = pro_d(j,kk-1,ii)
          ENDIF
          IF (rest1.GT.zero) THEN
            pro_f(j,kk,ii) = pro_f(j,kk,ii)/ rest1 ! NORMALIZE FRACTION
          ELSE
            pro_f(j,kk,ii) = 0.d0       ! IN CASE OF ALMOST TOTAL LOSS
          ENDIF
        ENDDO
!
        IF(rest1.GT.zero) THEN
          pro_f(j,kk,i) = 0.d0          ! FRACTION I IS REMOVED
        ELSE
          pro_f(j,kk,i) = 0.d0          ! IN CASE OF ALMOST TOTAL LOSS
        ENDIF
!
      ENDDO
!
      IF(cvsp_db(jg,0)) CALL cvsp_p('./VSP/','TEST5_',jg)
!
!-----------------------------------------------------------------------
! SPECIAL TREATMENT LOWEST POINT
!-----------------------------------------------------------------------
!
      rest1 = (-pro_f(j,lowpnt,i)+1.d0)
!
      DO ii=1,nsicla
        IF (rest1.GT.zero) THEN
          pro_f(j,lowpnt,ii) =
     &         pro_f(j,lowpnt,ii) / rest1    ! NORMALIZE FRACTION
        ELSE
          pro_f(j,lowpnt,ii) = 0.d0          ! IN CASE OF ALMOST TOTAL LOSS
        ENDIF
      ENDDO
!
      pro_f(j,lowpnt,i) = 0.d0                 ! FRACTION I IS REMOVED
      IF (rest1.LE.zero) pro_f(j,lowpnt,i) =
     &     0.d0                                ! IN CASE OF ALMOST TOTAL LOSS
!
!     ADDITIONAL CHECK FOR PRO_F
      DO k = 1, pro_max(j)
        ret = cvsp_check_f(j,k,'rm_fraction: ')
      ENDDO
!
      IF(cvsp_db(jg,0)) CALL cvsp_p('./','RFTLOW_',jg)
      CALL cvsp_compress_clean(j)
      IF(cvsp_db(jg,0)) CALL cvsp_p('./','RFTCCL_',jg)
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE cvsp_rm_fraction