init_sediment.f

Go to the documentation of this file.

!                   ************************
                    SUBROUTINE init_sediment
!                   ************************
!
     &(nsicla,elay,zf,zr,npoin,avail,fracsed_gf,ava0,
     & lgrafed,calwc,xmvs,xmve,grav,vce,xwc,fdm,
     & calac,ac,sedco,es,es_sable, es_vase ,nomblay,conc_vase,
     & ms_sable,ms_vase,acladm,unladm,toce_sable,
     & conc,debu,mixte)
!
!***********************************************************************
! SISYPHE   V7P2                                   27/06/2016
!***********************************************************************
!
!brief
!
!history  C. VILLARET (LNHE)
!+        30/12/2008
!+
!+
!
!history  JMH
!+        16/09/2009
!+        V6P0
!+   AVAIL(NPOIN,10,NSICLA)
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  C.VILLARET (EDF-LNHE), P.TASSI (EDF-LNHE)
!+        19/07/2011
!+        V6P1
!+  Name of variables
!+
!
!history  P.TASSI (EDF-LNHE)
!+        30/05/2012
!+        V6P2
!+  Case DSTAR > 150 AC(I) = 0.045D0
!+
!
!history  P.TASSI (EDF-LNHE)
!+        06/07/2012
!+        V6P2
!+  Line MIXTE=.FALSE. added.
!
!history  R.KOPMANN (BAW)
!+        27/06/2016
!+        V7P2
!+  CONTINUOUS APPROSIMATION CRITICAL SHILEDS CURVE DSTAR > 72 AC(I) = 0.045D0
!
!history  U.MERKEL AND R.KOPMANN (BAW)
!+        07/12/2017
!+        V7P2
!+  CALL TO CVSP_INIT
!
!history  R.KOPMANN (BAW)
!+        15/02/2019
!+        V7P2
!+  CALCULATION INITIAL VOLUME FOR SINGLE CLASS MODE
!+  ZEROING VARIABLES FOR MASS BALANCE
!
!history  R.KOPMANN (BAW)
!+        17/04/2019
!+        V8P2
!+  INITIALISING ES(:,1)=ELAY(:) IN CASE OF NSICLA=1
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AC             |<->| CRITICAL SHIELDS PARAMETER
!| ACLADM         |-->| MEAN DIAMETER OF SEDIMENT
!| AT0            |<->| TIME IN S
!| AVAIL          |<->| VOLUME PERCENT OF EACH CLASS
!| CALAC          |-->| IF YES, SHIELDS PARAMETER FOUND IN PARAMETER FILE
!| CALWC          |-->| IF YES, SETTLING VELOCITIES FOUND IN PARAMETER FILE
!| CONC_VASE      |<->| MUD CONCENTRATION FOR EACH LAYER
!| ELAY           |<->| THICKNESS OF SURFACE LAYER
!| ES             |<->| LAYER THICKNESSES AS DOUBLE PRECISION
!| ES_SABLE       |<->| LAYER THICKNESSES OF SAND AS DOUBLE PRECISION
!| ES_VASE        |<->| LAYER THICKNESSES OF MUD AS DOUBLE PRECISION
!| FDM            |-->| DIAMETER DM FOR EACH CLASS
!| FRACSED_GF     |-->|
!| GRAV           |-->| ACCELERATION OF GRAVITY
!| LGRAFED        |-->|
!| MS_SABLE       |<->| MASS OF SAND PER LAYER (KG/M2)
!| MS_VASE        |<->| MASS OF MUD PER LAYER (KG/M2)
!| ES_SABLE       |<->| THICKNESS OF SAND LAYER (M)
!| ES_VASE        |<->| THICKNESS OF MUD LAYER  (M)
!| MIXTE          |<->| SEDIMENT MIXTE  (SABLE /VASE)
!| NOMBLAY        |-->| NUMBER OF BED LAYERS
!| NPOIN          |-->| NUMBER OF POINTS
!| NSICLA         |-->| NUMBER OF SEDIMENT CLASSES
!| SEDCO          |-->| LOGICAL, SEDIMENT COHESIVE OR NOT
!| UNLADM         |-->| MEAN DIAMETER OF ACTIVE STRATUM LAYER
!| VCE            |-->| WATER VISCOSITY
!| XMVE           |-->| FLUID DENSITY
!| XMVS           |-->| WATER DENSITY
!| XWC            |-->| SETTLING VELOCITY
!| ZF             |-->| ELEVATION OF BOTTOM
!| ZR             |-->| NON ERODABLE BED
!| CONC           |<->| CONCENTRATION OF BED LAYER
!| DEBU           |-->| FLAG, RESTART ON SEDIMENTOLOGICAL FILE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE interface_sisyphe, ex_init_sediment => init_sediment
      USE declarations_sisyphe, ONLY : vsmtype
     & ,alt_model, pro_max_max, cvsmpperiod,nit,voltot,volu2d,volini
     & ,rmascl,vcumucl,volnestorcl

      USE interface_parallel, ONLY : p_dsum
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER,           INTENT(IN)     :: NSICLA,NPOIN,NOMBLAY
      TYPE(bief_obj),    INTENT(INOUT)  :: ELAY,ZF,ZR
      TYPE(bief_obj), INTENT(INOUT)     :: MS_SABLE, MS_VASE
      TYPE(bief_obj),    INTENT(INOUT)  :: ACLADM, UNLADM
      LOGICAL,           INTENT(IN)     :: LGRAFED,CALWC
      LOGICAL,           INTENT(IN)     :: CALAC
      DOUBLE PRECISION,  INTENT(IN)     :: XMVS,XMVE,GRAV,VCE
      DOUBLE PRECISION,  INTENT(INOUT)  :: AVA0(nsicla)
      DOUBLE PRECISION,  INTENT(INOUT)  :: AVAIL(npoin,nomblay,nsicla)
      DOUBLE PRECISION,  INTENT(INOUT)  :: FRACSED_GF(nsicla)
      DOUBLE PRECISION,  INTENT(INOUT)  :: FDM(nsicla),XWC(nsicla)
      DOUBLE PRECISION,  INTENT(INOUT)  :: AC(nsicla),TOCE_SABLE
      LOGICAL,           INTENT(IN)     :: SEDCO(nsicla), DEBU
      LOGICAL,           INTENT(IN)     :: MIXTE
      DOUBLE PRECISION, INTENT(IN)    :: CONC_VASE(nomblay)
      DOUBLE PRECISION, INTENT(INOUT) :: ES(npoin,nomblay)
      DOUBLE PRECISION, INTENT(INOUT) :: ES_SABLE(npoin,nomblay)
      DOUBLE PRECISION, INTENT(INOUT) :: ES_VASE(npoin,nomblay)
      DOUBLE PRECISION, INTENT(INOUT) :: CONC(npoin,nomblay)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER            :: I,J
      DOUBLE PRECISION   :: DENS,DSTAR
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
!  ------ BED COMPOSITION
!
      CALL os('X=Y-Z   ',x=elay,y=zf,z=zr)
!
!     ONLY ONE CLASS
!
      IF(nsicla.EQ.1) THEN
        DO i=1,npoin
          avail(i,1,1) = 1.d0
          acladm%R(i) = fdm(1)
        ENDDO
!
        IF(nomblay.EQ.1) THEN
          DO i=1,npoin
            es(i,1) = elay%R(i)
          ENDDO
        ENDIF
!
!     COMPUTES THE INITIAL VOLUME OF SEDIMENTS
        voltot(1)=0.d0
        DO j=1,npoin
            voltot(1) = voltot(1) + elay%R(j)*volu2d%R(j)
        ENDDO
        IF(ncsize.GT.1) THEN
          voltot(1) = p_dsum(voltot(1))
        ENDIF
        volini(1) = voltot(1)
!     INITIALISING VOLUMES
        rmascl = 0.d0
        vcumucl = 0.d0
        volnestorcl = 0.d0
!
!       PURE MUD ONLY
        IF(sedco(1)) CALL init_mixte(xmvs,npoin,avail,nsicla,es,
     &                               es_sable, es_vase,
     &                               elay%R,nomblay,conc_vase,
     &                                ms_sable%R,ms_vase%R,zf%R,
     &                               zr%R,ava0,conc,debu,.false.)
!
      ELSE
!
!     NON-COHESIVE, MULTI-CLASSES
!
        IF(.NOT.mixte) THEN
!
!
          CALL init_avai
!         CALL MEAN_GRAIN_SIZE
!         THIS PART CAN BE INTEGRATED INTO INIT_AVAI
          DO j=1,npoin
            acladm%R(j) = 0.d0
            unladm%R(j) = 0.d0
            DO i=1,nsicla
              IF(avail(j,1,i).GT.0.d0) THEN
                acladm%R(j) = acladm%R(j) + fdm(i)*avail(j,1,i)
                unladm%R(j) = unladm%R(j) + fdm(i)*avail(j,2,i)
              ENDIF
            ENDDO
            acladm%R(j)=max(acladm%R(j),0.d0)
            unladm%R(j)=max(unladm%R(j),0.d0)
          ENDDO
        ELSE
!
          CALL init_mixte(xmvs,npoin,avail,nsicla,es,
     &               es_sable, es_vase, elay%R,
     &               nomblay,conc_vase,ms_sable%R,
     &               ms_vase%R,zf%R,zr%R,ava0,conc,debu,mixte)
          DO i=1,npoin
            acladm%R(i) = fdm(1)
          ENDDO
        ENDIF
!
      ENDIF
!
      IF(lgrafed) THEN
        DO i=1, nsicla
          fracsed_gf(i)=ava0(i)
        ENDDO
      ENDIF
!
!     SETTLING VELOCITY
!
      IF(.NOT.calwc) THEN
        dens = (xmvs - xmve) / xmve
        DO i = 1, nsicla
          CALL vitchu_sisyphe(xwc(i),dens,fdm(i),grav,vce)
        ENDDO
      ENDIF
!
!     SHIELDS PARAMETER
!
      IF(.NOT.calac) THEN
        dens  = (xmvs - xmve )/ xmve
        DO i = 1, nsicla
          dstar = fdm(i)*(grav*dens/vce**2)**(1.d0/3.d0)
          IF (dstar <= 4.d0) THEN
            ac(i) = 0.24d0/dstar
          ELSEIF (dstar <= 10.d0) THEN
            ac(i) = 0.14d0*dstar**(-0.64d0)
          ELSEIF (dstar <= 20.d0) THEN
            ac(i) = 0.04d0*dstar**(-0.1d0)
!           CORRECTION 27/06/2016
!          ELSEIF (DSTAR <= 150.D0) THEN
          ELSEIF (dstar <= 72.d0) THEN
            ac(i) = 0.013d0*dstar**0.29d0
          ELSE
!           CORRECTION 30/05/2012
!           AC(I) = 0.055D0
            ac(i) = 0.045d0
          ENDIF
        ENDDO
      ENDIF
!
!     FOR MIXED SEDIMENTS
!
      IF(mixte) toce_sable=ac(1)*fdm(1)*grav*(xmvs - xmve)
!
!-----------------------------------------------------------------------
!
      IF(vsmtype.EQ.1) THEN
        WRITE(lu,*) ' '
        WRITE(lu,*) '--------------------------------------------------'
        WRITE(lu,*) 'C-VSM MODEL'
        WRITE(lu,*) 'CONTINUOUS VERTICAL GRAIN SORTING STRATIGRAPHY'
        WRITE(lu,*) ' '
        WRITE(lu,*) 'ACTIVE LAYER THICKNESS MODEL:', alt_model
!
        IF (pro_max_max .GT. 250) THEN
          WRITE(lu,*) 'HIGH NUMBER OF SECTIONS IS EXPENSIVE',pro_max_max
          WRITE(lu,*) 'BETTER < 250 and > 4 + 4 x NUMBER OF FRACTIONS'
        ENDIF
!
        IF ((cvsmpperiod / nit) .GT. 5) THEN
          WRITE(lu,*) 'HIGH NUMBER OF FULL CVSM PRINTOUTS'
          WRITE(lu,*) 'ATTENTION: DISK SPACE AND SIMULATION TIME'
          WRITE(lu,*) 'ADAPT C-VSM FULL PRINTOUT PERIOD'
        ENDIF
!
        WRITE(lu,*) ' '
        CALL cvsp_init()
        WRITE(lu,*)'CVSM INITIALISED!'
        WRITE(lu,*) '--------------------------------------------------'
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END