bedload_main.f

Go to the documentation of this file.

!                   ***********************
                    SUBROUTINE bedload_main
!                   ***********************
!
     &(acladm,ksp,ksr, v2dpar,unsv2d,cf,ebor,fw,hn,liqbor,
     & mask, maskel, maskpt, qbor, u2d,
     & v2d, s,unladm,uw,thetaw,mu,tob,tobw,tw,zf,
     & debug, hidfac, icf, ielmt, kddl, kdir,
     & kent, klog, kneu, ksort,
     & npoin, nptfr, nsicla, optban, beta, fd90, fdm,
     & grav, hidi, hmin, vce, csf_sable, xmve, xmvs, xwc,
     & pi, karman, zero, karim_holly_yang,msk, susp, vf,
     & mesh,
     & elay, liebor, limtec, masktr,
     & it1, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11,
     & t12,t13,unorm,ac, dts,
     & avail, breach, calfa_cl, coefpn,
     & hiding, qscl_c, qscl_s, qs_c,
     & qsclxc, qsxc, qsclyc, qsyc, salfa_cl, zf_c, zfcl_c,
     & entets, seccurrent, slopeff,
     & phised, devia, beta2, bijk,sedco,houle,
     & u3d,v3d,code,flbcla,maxadv)
!
!***********************************************************************
! SISYPHE   V8P0                                   12/09/2018
!***********************************************************************
!
!brief    MAIN SUBROUTINE FOR THE BEDLOAD TRANSPORT.
!
!history  F. HUVELIN
!+        14/09/2004
!+
!+
!
!history  JMH
!+        21/12/2006
!+        V6P0
!+   BEDLOAD_TIMESTEP NO LONGER EXISTS
!
!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  J-M HERVOUET (EDF LAB, LNHE)
!+        28/03/2014
!+        V7P0
!+  Call to bedload_diffin changed.
!history  R KOPMANN (BAW)
!+        10/05/2016
!+        V7P2
!+ CALFA,SALFA dependent of grain classes
!
!history  F.CORDIER & P. TASSI(EDF-LNHE)
!+        12/09/2018
!+        V8P0
!+  Calculation of SANFRA for the formula of WILCOCK and CROWE (2003)
!
!history  R.KOPMANN (BAW)
!+        15/02/2019
!+        V7P2
!+  Calculating mass changes from Nestor
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AC             |<->| CRITICAL SHIELDS PARAMETER
!| ACLADM         |-->| MEAN DIAMETER OF SEDIMENT
!| AVAIL          |<->| VOLUME PERCENT OF EACH CLASS
!| BETA           |-->| COEFFICIENT FOR SLOPING BED EFFECT ( KOCH AND FLOKSTRA)
!| BETA2          |-->| COEFFICIENT FOR THE DEVIATION  (TALMON ET AL.)
!| BIJK           |-->| COEFFICIENT OF THE BIJKER FORMULA
!| BREACH         |<->| INDICATOR FOR NON ERODIBLE BED (FINITE VOLUMES SHEMES)
!| CALFA          |<->| COSINUS OF THE ANGLE BETWEEN MEAN FLOW AND TRANSPORT
!| CF             |-->| QUADRATIC FRICTION COEFFICIENT
!| CODE           |-->| HYDRODYNAMIC CODE IN CASE OF COUPLING
!| COEFPN         |<->| CORRECTION OF TRANSORT FOR SLOPING BED EFFECT
!| DEBUG          |-->| FLAG FOR DEBUGGING
!| DEVIA          |-->| SLOPE EFFECT FORMULA FOR DEVIATION
!| DTS            |<->| TIME STEP FOR SUSPENSION
!| EBOR           |<->| IMPOSED BOUNDARY CONDITION FOR BED EVOLUTION (DIRICHLET)
!| ELAY           |<->| THICKNESS OF SURFACE LAYER
!| ENTETS         |<->| LOGICAL, IF YES INFORMATION IS GIVEN ON MASS CONSERVATION FOR SUSPENSION
!| FD90           |-->| DIAMETER D90
!| FDM            |-->| DIAMETER DM FOR EACH CLASS
!| FLBCLA         |-->| BLOCK OF FLUXES AT BOUNDARY FOR EACH CLASS
!| FW             |-->| WAVE FRICTION FACTOR
!| GRAV           |-->| ACCELERATION OF GRAVITY
!| HIDFAC         |-->| HIDING FACTOR FORMULAS
!| HIDI           |-->| HIDING FACTOR FOR PARTICULAR SIZE CLASS (HIDFAC =0)
!| HIDING         |-->| HIDING FACTOR CORRECTION
!| HMIN           |-->| MINIMUM VALUE OF WATER DEPTH
!| HN             |-->| WATER DEPTH
!| HOULE          |-->| LOGICAL, FOR WAVE EFFECTS
!| ICF            |-->| BED-LOAD OR TOTAL LOAD TRANSPORT FORMULAS
!| IELMT          |-->| NUMBER OF ELEMENTS
!| IT1            |<->| INTEGER WORK ARRAY IN A BIEF_OBJ STRUCTURE
!| KARMAN         |-->| VON KARMAN CONSTANT
!| KDDL           |-->| CONVENTION FOR DEGREE OF FREEDOM
!| KDIR           |-->| CONVENTION FOR DIRICHLET POINT
!| KENT           |-->| CONVENTION FOR LIQUID INPUT WITH PRESCRIBED VALUE
!| KLOG           |-->| CONVENTION FOR SOLID BOUNDARY
!| KNEU           |-->| CONVENTION FOR NEUMANN CONDITION
!| KSORT          |-->| CONVENTION FOR FREE OUTPUT
!| KSP            |-->| BED SKIN ROUGHNESS
!| KSR            |-->| RIPPLE BED ROUGHNESS
!| LIEBOR         |<->| TYPE OF BOUNDARY CONDITIONS FOR BED EVOLUTION
!| LIMTEC         |<->| TECHNICAL BOUNDARY CONDITION (NEUMAN...)
!| LIQBOR         |-->| TYPE OF BOUNDARY CONDITION FOR QS
!| MASK           |-->| BLOCK OF MASKS, EVERY ONE FOR A TYPE OF BOUNDARY
!|                |   | SEE DIFFIN.F IN LIBRARY BIEF.
!| MASKEL         |-->| MASKING OF ELEMENTS
!| MASKPT         |-->| MASKING PER POINT
!| MASKTR         |<->| MASKING FOR TRACERS, PER POINT
!| MESH           |<->| MESH STRUCTURE
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS
!| MU             |<->| CORRECTION FACTOR FOR BED ROUGHNESS
!| NPOIN          |-->| NUMBER OF POINTS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| NSICLA         |-->| NUMBER OF SIZE CLASSES FOR BED MATERIALS
!| OPTBAN         |-->| OPTION FOR TIDAL FLATS
!| PHISED         |-->| ANGLE OF REPOSE OF THE SEDIMENT
!| PI             |-->| PI
!| QBOR           |-->| BOUNDARY CONDITION FOR TRANSPORT RATE
!| QSCLXC         |<->| TRANSPORT RATE FOR EACH CLASS X-DIRECTION
!| QSCLYC         |<->| TRANSPORT RATE FOR EACH CLASS Y-DIRECTION
!| QSCL_C         |<->| BEDLOAD TRANSPORT RATE
!| QSCL_S         |<->| SUSPENDED LOAD TRANSPORT RATE
!| QSXC           |<->| BEDLOAD TRANSPORT RATE X-DIRECTION
!| QSYC           |<->| BEDLOAD TRANSPORT RATE Y-DIRECTION
!| QS_C           |<->| BEDLOAD TRANSPORT RATE
!| S              |-->| VOID STRUCTURE
!| SALFA          |<->| SINUS OF THE ANGLE BETWEEN TRANSPORT RATE AND CURRENT
!| SANFRA         |-->| SAND FRACTION (WILCOCK AND CROWE, 2003)
!| SECCURRENT     |-->| LOGICAL, PARAMETRISATION FOR SECONDARY CURRENTS
!| SEDCO          |-->| LOGICAL, SEDIMENT COHESIVE OR NOT
!| SLOPEFF        |-->| LOGICAL, SLOPING BED EFFECT OR NOT
!| SUSP           |-->| LOGICAL, SUSPENSION
!| T1             |<->| WORK BIEF_OBJ STRUCTURE
!| T10            |<->| WORK BIEF_OBJ STRUCTURE
!| T11            |<->| WORK BIEF_OBJ STRUCTURE
!| T12            |<->| WORK BIEF_OBJ STRUCTURE
!| T13            |<->| WORK BIEF_OBJ STRUCTURE
!| T2             |<->| WORK BIEF_OBJ STRUCTURE
!| T3             |<->| WORK BIEF_OBJ STRUCTURE
!| T4             |<->| WORK BIEF_OBJ STRUCTURE
!| T5             |<->| WORK BIEF_OBJ STRUCTURE
!| T6             |<->| WORK BIEF_OBJ STRUCTURE
!| T7             |<->| WORK BIEF_OBJ STRUCTURE
!| T8             |<->| WORK BIEF_OBJ STRUCTURE
!| T9             |<->| WORK BIEF_OBJ STRUCTURE
!| THETAW         |-->| ANGLE BETWEEN WAVE AND CURRENT
!| TOB            |<->| BED SHEAR STRESS (TOTAL FRICTION)
!| TOBW           |-->| WAVE INDUCED SHEAR STRESS
!| TW             |-->| WAVE PERIOD
!| U2D            |<->| MEAN FLOW VELOCITY X-DIRECTION
!| U3D            |-->| THREE-DIMENSIONAL VELOCITY X-DIRECTION
!| UNLADM         |-->| MEAN DIAMETER OF ACTIVE STRATUM LAYER
!| UNORM          |<->| NORM OF THE MEAN FLOW VELOCITY
!| UNSV2D         |-->| INVERSE OF INTEGRALS OF TEST FUNCTIONS
!| UW             |-->| ORBITAL WAVE VELOCITY
!| V2D            |<->| MEAN FLOW VELOCITY Y-DIRECTION
!| V2DPAR         |-->| INTEGRAL OF TEST FUNCTIONS, ASSEMBLED IN PARALLEL
!| V3D            |-->| THREE-DIMENSIONAL VELOCITY Y-DIRECTION
!| VCE            |-->| WATER VISCOSITY
!| VF             |-->| LOGICAL, FINITE VOLUMES OR NOT
!| CSF_SABLE      |-->| BED VOLUME CONCENTRATION CSF = (1-POROSITY)
!| XMVE           |-->| FLUID DENSITY
!| XMVS           |-->| SEDIMENT DENSITY
!| XWC            |-->| SETTLING VELOCITY
!| ZERO           |-->| ZERO
!| ZF             |-->| ELEVATION OF BOTTOM
!| ZFCL_C         |<->| BEDLOAD EVOLUTION FOR EACH SEDIMENT CLASS
!| ZF_C           |<->| BEDLOAD EVOLUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE interface_sisyphe, ex_bedload_main => bedload_main
      USE declarations_sisyphe, ONLY : nestor,nomblay,sanfra
      USE declarations_sisyphe, ONLY : volu2d,volnestorcl,volnestorcla
      USE declarations_special
      USE interface_parallel, ONLY : p_dsum
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_obj),   INTENT(IN)    :: ACLADM, KSR,V2DPAR,UNSV2D
      TYPE(bief_obj),   INTENT(IN)    :: CF,FW,KSP,HN,LIQBOR
      TYPE(bief_obj),   INTENT(IN)    :: MASK, MASKEL, MASKPT
      TYPE(bief_obj),   INTENT(IN)    :: QBOR
      TYPE(bief_obj),   INTENT(INOUT) :: U2D,V2D,TOB,MU,UNORM,EBOR
      TYPE(bief_obj),   INTENT(IN)    :: S,UNLADM
      TYPE(bief_obj),   INTENT(IN)    :: UW, THETAW,  TOBW, TW
      TYPE(bief_obj),   INTENT(IN)    :: ZF
      INTEGER,          INTENT(IN)    :: DEBUG, HIDFAC, ICF,MAXADV
      INTEGER,          INTENT(IN)    :: IELMT, KDDL, KDIR, KENT
      INTEGER,          INTENT(IN)    :: KLOG, KNEU, KSORT
      INTEGER,          INTENT(IN)    :: NPOIN, NPTFR
      INTEGER,          INTENT(IN)    :: NSICLA, OPTBAN
      DOUBLE PRECISION, INTENT(IN)    :: BETA, FD90(nsicla)
      DOUBLE PRECISION, INTENT(IN)    :: FDM(nsicla),GRAV
      DOUBLE PRECISION, INTENT(IN)    :: HIDI(nsicla),HMIN,VCE
      DOUBLE PRECISION, INTENT(IN)    :: CSF_SABLE,XMVE,XMVS,XWC(nsicla)
      DOUBLE PRECISION, INTENT(IN)    :: PI,KARMAN,ZERO
      DOUBLE PRECISION, INTENT(IN)    :: KARIM_HOLLY_YANG
      LOGICAL,          INTENT(IN)    :: MSK, SUSP, VF
      LOGICAL,          INTENT(IN)    :: SECCURRENT
      LOGICAL,          INTENT(IN)    :: SEDCO(nsicla),HOULE
      TYPE(bief_mesh),  INTENT(INOUT) :: MESH
      TYPE(bief_obj),   INTENT(INOUT) :: ELAY,FLBCLA
      TYPE(bief_obj),   INTENT(INOUT) :: LIEBOR, LIMTEC, MASKTR
      TYPE(bief_obj),   INTENT(INOUT) :: IT1,T1,T2,T3,T4,T5,T6,T7
      TYPE(bief_obj),   INTENT(INOUT) :: T8,T9,T10,T11,T12,T13
      DOUBLE PRECISION, INTENT(INOUT) :: AC(nsicla), DTS
      DOUBLE PRECISION, INTENT(INOUT) :: AVAIL(npoin,nomblay,nsicla)
      TYPE(bief_obj),   INTENT(INOUT) :: BREACH, CALFA_CL, COEFPN
      TYPE(bief_obj),   INTENT(INOUT) :: HIDING
      TYPE(bief_obj),   INTENT(INOUT) :: QSCL_C,QSCL_S
      TYPE(bief_obj),   INTENT(INOUT) :: QS_C, QSCLXC, QSXC, QSCLYC
      TYPE(bief_obj),   INTENT(INOUT) :: QSYC, SALFA_CL, ZF_C, ZFCL_C
      LOGICAL,          INTENT(INOUT) :: ENTETS
      DOUBLE PRECISION,   INTENT(IN)  :: BETA2, PHISED
      INTEGER, INTENT (IN)            :: SLOPEFF, DEVIA
      DOUBLE PRECISION, INTENT(IN)    :: BIJK
      TYPE(bief_obj),    INTENT(IN)    :: U3D,V3D
      CHARACTER(LEN=24), INTENT(IN)    :: CODE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I,K,J
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
!
! FC
!     Calculation of sand fraction content at each node (Wilcock and Crowe, 2003)
      IF (icf == 10) THEN
        DO k = 1, npoin
          sanfra(k) = 0.0d0
          DO i = 1, nsicla
            IF (fdm(i).LT.2.d-3) THEN
              sanfra(k) = sanfra(k) + avail(k,1,i)
            ENDIF
          ENDDO
        ENDDO
      ENDIF
! ENDFC

!     INITIALISES TECHNICAL BOUNDARY CONDITIONS
!
      IF (debug > 0) WRITE(lu,*) 'BEDLOAD_DIFFIN'
      CALL bedload_diffin
     &        (u2d, v2d, mesh%NBOR, mesh%XNEBOR, mesh%YNEBOR,
     &         maskel, mesh%NELBOR, nptfr, kent, ksort, klog,
     &         kdir, kddl, kneu, msk, it1, liebor, masktr,limtec,
     &         mesh%IKLBOR%I,mesh%NELEB,mesh%NELEBX)
      IF (debug > 0) WRITE(lu,*) 'END_BEDLOAD_DIFFIN'
!
      DO i = 1, nsicla
!
!       FOR SAND
        IF(.NOT.sedco(i)) THEN
          IF (debug > 0) WRITE(lu,*)
     &      'BEDLOAD_SOLIDISCHARGE : ',i,'/',nsicla
          CALL bedload_solidischarge
     &       (mesh, u2d, v2d, unorm,hn, tw, uw, mu,tob,cf,
     &         tobw,fw,thetaw,avail(1:npoin,1,i),
     &         maskpt, maskel, acladm,
     &         unladm,ksp,ksr, liqbor, debug, npoin,
     &         nptfr, ielmt, icf, kent, optban, hidfac, grav,
     &         fdm(i), fd90(i), xwc(i), xmve, xmvs, vce, hmin,
     &         hidi(i),karman,zero,pi,
     &         karim_holly_yang,susp,msk,t1,t2,
     &         t3, t4, t5, t6, t7, t8, t9, t10, t11,t12, ac(i),
     &         hiding,qscl_c%ADR(i)%P,qscl_s%ADR(i)%P,
     &         slopeff,coefpn,phised,
     &         calfa_cl%ADR(i)%P,salfa_cl%ADR(i)%P,
     &         beta,zf,s,
     &         devia, beta2 , seccurrent, bijk,houle,unsv2d,
     &         u3d,v3d,code,sanfra)
          IF(debug > 0) WRITE(lu,*) 'END_BEDLOAD_SOLIDISCHARGE'
        ELSE
!         FOR COHESIVE SEDIMENT: ZERO BEDLOAD TRANSPORT RATE
!         TODO: IS THIS USEFUL ???
          CALL os('X=0     ',x=qscl_c%ADR(i)%P)
          CALL os('X=0     ',x=qsclxc%ADR(i)%P)
          CALL os('X=0     ',x=qsclyc%ADR(i)%P)
        ENDIF
!
      ENDDO
!
!     COMPUTES THE EVOLUTION FOR EACH CLASS
!
      DO i = 1, nsicla
!
        IF(.NOT.sedco(i)) THEN
!
          IF (debug > 0) WRITE(lu,*) 'BEDLOAD_EVOL : ',i,'/',nsicla
          CALL bedload_evol(s,elay,avail(1:npoin,1,i),
     &                      coefpn,calfa_cl%ADR(i)%P,salfa_cl%ADR(i)%P,
     &                      limtec,
     &                      ebor%ADR(i)%P,maskel,mask,
     &                      v2dpar,unsv2d,debug,npoin,nptfr,ielmt,
     &                      kent,kdir,kddl,
     &                      dts,
     &                      vf,entets,msk,
     &                      mesh,qscl_c%ADR(i)%P,
     &                      t1,t2,t3,t4,t8,t11,t12,
     &                      t13,csf_sable,breach,qsclxc%ADR(i)%P,
     &                      qsclyc%ADR(i)%P,zfcl_c%ADR(i)%P,slopeff,
     &                      i,flbcla,liqbor,qbor%ADR(i)%P,maxadv)
          IF(debug.GT.0) WRITE(lu,*) 'END_BEDLOAD_EVOL'
!
!         NOW DIVIDING BY CSF_SABLE TO GET THE EVOLUTION OF BED
!         INCLUDING VOIDS
!         TODO: NOTE JMH: IN BEDLOAD_EVOL THERE IS A PRELIMINARY MULTIPLICATION BY
!                   CSF_SABLE, SO THIS COULD BE SIMPLIFIED, BUT FOR THE
!                   FINITE ELEMENT OPTION ONLY, THE FINITE VOLUME IMPLEMENTATION
!                   SEEMS MORE COMPLICATED TO SORT OUT.
!
          CALL os('X=CX    ',x= zfcl_c%ADR(i)%P,c=1.d0/csf_sable)
!
        ELSE
!
!         NO EVOLUTION FOR COHESIVE SEDIMENT
          CALL os('X=0     ',x=zfcl_c%ADR(i)%P)
!
        ENDIF
!
      ENDDO
!
      IF(nestor) THEN
        DO i = 1,nsicla
            t8%R(i) = 0.d0
          DO j=1,npoin
            t8%R(i) = t8%R(i) + zfcl_c%ADR(i)%P%R(j)*volu2d%R(j)
          END DO
          IF(ncsize>1) t8%R(i) = p_dsum(t8%R(i))
          volnestorcla(i) = t8%R(i)
        END DO
        CALL nestor_interface(2)
        DO i = 1,nsicla
          t9%R(i) = 0.d0
          DO j=1,npoin
            t9%R(i) = t9%R(i) + zfcl_c%ADR(i)%P%R(j)*volu2d%R(j)
          END DO
          IF(ncsize>1) t9%R(i) = p_dsum(t9%R(i))
!         CUMULATED NESTOR CHANGES PER CLASS
          volnestorcl(i) = volnestorcl(i)+t9%R(i)-volnestorcla(i)
!          NESTOR CHANGES PER TIME STEP AND CLASS
          volnestorcla(i) = t9%R(i)-volnestorcla(i)
        END DO
      ENDIF
!
      ! *********************************************** !
      ! II - EVOLUTIONS AND QS FOR EACH CLASS ARE ADDED !
      ! *********************************************** !
      ! II.1 - INITIALISES
      ! ---------------------
      CALL os('X=0     ', x=qs_c)
      CALL os('X=0     ', x=zf_c)
      CALL os('X=0     ',x=qsxc)
      CALL os('X=0     ',x=qsyc)
      ! II.2 - ADDS THE CLASSES
      ! ----------------------
      !
      DO i=1,nsicla
        IF(.NOT.sedco(i)) THEN
          CALL os('X=X+Y   ', x=qs_c, y=qscl_c%ADR(i)%P)
          CALL os('X=X+Y   ', x=zf_c, y=zfcl_c%ADR(i)%P)
          CALL os('X=X+YZ  ', x=qsxc, y=qscl_c%ADR(i)%P,
     &                              z=calfa_cl%ADR(i)%P)
          CALL os('X=X+YZ  ', x=qsyc, y=qscl_c%ADR(i)%P,
     &                              z=salfa_cl%ADR(i)%P)

        ENDIF
      ENDDO
!
!     TIDAL FLATS WITH MASKING
!
      IF(optban.EQ.2) CALL os('X=XY    ',x=zf_c,y=maskpt)
!
!======================================================================!
!======================================================================!
!
      RETURN
      END