cvsp_alt.f

Go to the documentation of this file.

!                 **********************************
                  DOUBLE PRECISION FUNCTION cvsp_alt
!                 **********************************
!
     &(j, formula)
!
!***********************************************************************
! SISYPHE   V7P2                                   05/12/2017
!***********************************************************************
!
!brief   CALCULATES A DYNAMIC ACTIVE LAYER THICKNESS
!+        ACCORDING TO 1 OF A COUPLE OF FORMULAS
!
!
!history  UWE MERKEL
!+        20/07/2011
!+       V6P2
!+
!
!history  P. A. TASSI (EDF R&D, LNHE)
!+        12/03/2013
!+        V6P3
!+   cleaning, cosmetic
!
!history  LEOPOLD STADLER (BAW) & J-M HERVOUET (EDF LAB, LNHE)
!+        28/07/2014
!+        V7P0
!+   Computation of D90 and DIAMAX secured to avoid divisions by 0.
!
!history  UWE MERKEL (UHM) R. KOPMANN (BAW)
!+        22/11/2016 / 2017
!+        V6P3 / V7P2
!+   Many changes
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| J              |<--| INDEX OF A POINT IN MESH
!| FORMULA        |<--| WHICH FORMULA TO USE TO CALCULATE THE ALT
!| VCE            |-->| WATER VISCOSITY
!| XMVE           |-->| FLUID DENSITY
!| XMVS           |-->| SEDIMENT DENSITY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_sisyphe
!
      USE declarations_special
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: J
      INTEGER, INTENT(IN) :: FORMULA
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER  I
      DOUBLE PRECISION RHO, RHO_S, G, D50, D90, DMAX, FMAX, TAUC, TAUB
      DOUBLE PRECISION PON, SUMME
      DOUBLE PRECISION DSTAR, DENS
!
!-----------------------------------------------------------------------
!
! ATTENTION !
!
!-----------------------------------------------------------------------
!
! EXPECTS GRAIN CLASSES TO BE SORTED IN ASCENDING ORDER
! UNLIKE SISYPHE!!!
! IMPROVE IT!!!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
! CHECK ASCENDING ORDER OF CLASSES D50
!-----------------------------------------------------------------------
      DO i=1,nsicla-1
        IF(fdm(i).GE.fdm(i+1)) THEN
          WRITE(lu,*) 'STOPPING!!!! GRAIN CLASSES HAVE TO BE',
     &                ' IN ASCENDING ORDER!!! FOR DYNAMIC ALT'
          CALL plante(1)
          stop
        ENDIF
      ENDDO
!-----------------------------------------------------------------------
! BASICS
!-----------------------------------------------------------------------
      g = grav
      rho = xmve
      rho_s = xmvs
      pon = xkv
      d50 = acladm%R(j)
!-----------------------------------------------------------------------
! DMAX - NEW APPROXIMATION: DMAX = D99 / BEFORE DMAX = MAX(FDM)
!-----------------------------------------------------------------------
      fmax = 0.99d0
      summe = avail(j,1,1)
      dmax = 0.d0
      IF (summe.GE.fmax) THEN
        dmax = fdm(1)
      ELSE
        DO i=2,nsicla
          IF(avail(j,1,i).GT.0.d0) THEN
            summe = avail(j,1,i) + summe
            IF(summe.GE.fmax.AND.dmax.EQ.0.d0) THEN
              dmax = fdm(i-1) + ( (fdm(i)-fdm(i-1)) /
     &        avail(j,1,i)) * (fmax-(summe-avail(j,1,i)))
            ENDIF
          ENDIF
        ENDDO
      ENDIF
!DMAX = MAX(FDM(1), DMAX)
      IF (dmax.LE.0.d0.AND.summe.GT.zero) THEN
        WRITE(lu,*)'UHM DMAXERROR'
        DO i=1,nsicla
          WRITE(lu,*)'DMAXERROR',j,i,fdm(i),avail(j,1,i),summe
          WRITE(lu,*)zf%R(j)-zr%R(j)
        ENDDO
        stop
      ENDIF
!-----------------------------------------------------------------------
! D90 - ONLY FIRST APPROXIMATION!
!-----------------------------------------------------------------------
      fmax = 0.90d0
      summe = avail(j,1,1)
      d90 = 0.d0
      IF (summe.GE.fmax) THEN
        d90 = fdm(1)
      ELSE
        DO i=2,nsicla
          IF(avail(j,1,i).GT.0.d0) THEN
            summe = avail(j,1,i) + summe
            IF(summe.GE.fmax.AND.d90.EQ.0.d0) THEN
              d90 = fdm(i-1) + ( (fdm(i)-fdm(i-1)) /
     &              avail(j,1,i)) * (fmax-(summe-avail(j,1,i)))
            ENDIF
          ENDIF
        ENDDO
      ENDIF
!D90 = MAX(FDM(1), D90)
      IF (d90.LE.0.d0.AND.summe.GT.zero) THEN
        WRITE(lu,*)'UHM D90ERROR', d90, fmax, summe
        DO i=1,nsicla
          WRITE(lu,*)'D90ERROR',j,i,fdm(i),avail(j,1,i),summe
        ENDDO
        WRITE(lu,*)'D90ERROR'
      ENDIF
!-----------------------------------------------------------------------
! SHEAR PARAMETERS
!-----------------------------------------------------------------------
! HERE ARE ENOUGH POSSIBILITIES FOR IMPROVEMENT
! THE INITIATION OF MOTION STARTS WITH THE SMALLEST GRAIN
! LEEDS TO BIGGER ACTIVE LAYER THICKNESSES
      tauc = ac(1)*((xmvs-xmve)*g*d50)
      taub = tob%R(j)
!
!-----------------------------------------------------------------------
! NEW ACTIVE LAYER THICKNESS
!-----------------------------------------------------------------------
!
      SELECT CASE (formula)
!-----------------------------------------------------------------------
! HUNZIKER & GUENTHER
!-----------------------------------------------------------------------
!
      CASE (1)
        cvsp_alt = 5.d0 * dmax
!
!-----------------------------------------------------------------------
! FREDSOE & DEIGAARD 1992
!-----------------------------------------------------------------------
!
      CASE (2)
        cvsp_alt = 2.d0 * taub / (g*(rho_s-rho))
     &     / tan(phised/180.0d0*pi) / (1.d0-pon)
!
!-----------------------------------------------------------------------
! VAN RIJN 1993
!-----------------------------------------------------------------------
!
      CASE (3)
        cvsp_alt = 0.d0
        IF(tauc.GT.0.d0) THEN
          IF(taub.GT.tauc) THEN
            dens = (rho_s - rho) / rho
            dstar = d50*(g*dens/vce**2)**(1.d0/3.d0)
            cvsp_alt = 0.3d0*(dstar**0.7d0)*
     &         ((taub-tauc)/tauc)**0.5*d50
          ENDIF
        ENDIF
!-----------------------------------------------------------------------
! WONG 2006
!-----------------------------------------------------------------------
      CASE (4)
        IF((taub/(rho_s-rho)/g/d50).LT.0.0549d0) THEN
          cvsp_alt = 0.d0
        ELSE
          cvsp_alt=5.0d0*d50*((taub/(rho_s-rho)/g/d50)
     &          -0.0549d0)**0.56d0
        ENDIF
!-----------------------------------------------------------------------
! MALCHEREK 2003
!-----------------------------------------------------------------------
      CASE (5)
        cvsp_alt = 0.d0
        IF (tauc.GT.0) THEN
          cvsp_alt = d90 / (1.d0-pon) * taub/tauc
        ENDIF
!-----------------------------------------------------------------------
! SISYPHE
!-----------------------------------------------------------------------
      CASE (6)
        cvsp_alt = 3.d0 * d50
!-----------------------------------------------------------------------
! CONSTANT FROM CAS FILE
!-----------------------------------------------------------------------
      CASE (0)
        cvsp_alt = elay0
!-----------------------------------------------------------------------
      CASE DEFAULT
        WRITE(lu,*)'NO VALID CHOICE FOR "C-VSM DYNAMIC ALT MODEL"'
        WRITE(lu,*)'MUST BE BETWEEN 0 AND 6'
        CALL plante(1)
!
      END SELECT
!-----------------------------------------------------------------------
! MINIMUM ALT VALUE = SMALLEST GRAIN SIZE
!-----------------------------------------------------------------------
      cvsp_alt = max(cvsp_alt, fdm(1))
!
!-----------------------------------------------------------------------
!
      RETURN
      END FUNCTION cvsp_alt