tob_sisyphe.f

Go to the documentation of this file.

!                   **********************
                    SUBROUTINE tob_sisyphe
!                   **********************
!
     & (tob, tobw, mu, ks,ksp, ksr,cf,fw,chestr,uetcar,
     &  cf_tel,ks_tel,code,
     &  kfrot,icr, kspratio, houle,grav,xmve,xmvs, vce, karman,
     &  zero,hmin,hn, acladm, unorm,uw, tw, npoin,kspred,iks)
!
!***********************************************************************
! SISYPHE   V7P2
!***********************************************************************
!
!brief    COMPUTES THE TOTAL STRESS AT THE BOTTOM DEPENDING
!+                ON WHETHER SISYPHE IS COUPLED OR NOT.
!
!history  CV
!+        **/04/05
!+
!+   CORRECTION WHEN SISYPHE IS RUN ALONE: DO NOT MODIFY
!
!history  C. VILLARET (LNHE)
!+        29/11/06
!+        V6P0
!+
!
!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), ON BEHALF OF CLEMENS DORFMAN
!+        12/02/2016
!+        V7P2
!+  Changing the constant 12.D0 into 11.036D0 because in the rest of the
!+  code the Karman constant has been changed from 0.41 to 0.40.
!
!history  RIADH ATA
!+        12/03/2018
!+        V8P0
!+  Change the call of OS by a call of OV since, when coupling, this
!+  subroutine can be called by vectors and not structures (arguments)
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| ACLADM         |-->| MEAN DIAMETER OF SEDIMENT
!| CF             |-->| QUADRATIC FRICTION COEFFICIENT
!| CF_TEL         |-->| QUADRATIC FRICTION COEFFICIENT (COUPLED T2D)
!| CHESTR         |-->| FRICTION COEFFICIENT (KEYWORD)
!| CODE           |-->| CALLING PROGRAM IN COUPLING
!| FW             |-->| QUADRATIC FRICTION COEFFICIENT (WAVE)
!| GRAV           |-->| ACCELERATION OF GRAVITY
!| HMIN           |-->| MINIMUM VALUE OF WATER DEPTH
!| HN             |-->| WATER DEPTH
!| HOULE          |-->| INCLUDE WAVES COMPUTATIONS
!| ICR            |-->| ICR=0: MU=1
!|                |   | ICR=1: SKIN FRICTION CORRECTION USE KSP
!|                |   | ICR=2: RIPPLE ROUGHNESS USE KSR, KSR
!| KARMAN         |-->| VON KARMAN CONSTANT
!| KFROT          |-->| FRICTION LAW
!| KS             |<--| RUGOSITE TOTALE
!| KSP            |<--| RUGOSITE DE PEAU
!| KSPRATIO       |-->| RATIO BETWEEN SKIN BED ROUGHNESS AND GRAIN DIAMETER
!| KSR            |<--| RUGOSITE DE RIDE
!| KS_TEL         |<--| RUGOSITE TOTALE
!| MU             |<->| CORRECTION FACTOR FOR BED ROUGHNESS
!| NPOIN          |-->| NUMBER OF POINTS
!| TOB            |<->| BED SHEAR STRESS (TOTAL FRICTION)
!| TOBW           |-->| WAVE INDUCED SHEAR STRESS
!| TW,UW          |-->| WAVE PERIOD AND ORBITAL VELOCITY
!| UETCAR         |-->| SQUARE OF THE FRICTION VELOCITY (COUPLED T3D)
!| UNORM          |-->| INTENSITE DU COURANT
!| VCE            |-->| WATER VISCOSITY
!| XMVE           |-->| FLUID DENSITY (MASS)
!| XMVS           |-->| SEDIMENT DENSITY (MASS)
!| ZERO           |-->| ZERO
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE interface_sisyphe, ex_tob_sisyphe=>tob_sisyphe
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER,            INTENT(IN)  :: NPOIN,KFROT,ICR, IKS
      LOGICAL,            INTENT(IN)  :: KSPRED
      LOGICAL,            INTENT(IN)  :: HOULE
      CHARACTER(LEN=24),  INTENT(IN)  :: CODE
      DOUBLE PRECISION,   INTENT(IN)  :: XMVE,XMVS, VCE,GRAV,KARMAN
      DOUBLE PRECISION,   INTENT(IN)  :: ZERO,HMIN,KSPRATIO
      TYPE(bief_obj), INTENT(IN)      :: UETCAR
      TYPE(bief_obj), INTENT(IN)      :: HN,UNORM
      TYPE(bief_obj), INTENT(IN)      :: TW,UW
      TYPE(bief_obj), INTENT(INOUT)   :: KS,KSP,KSR,KS_TEL
      TYPE(bief_obj), INTENT(INOUT)   :: CHESTR,MU
      TYPE(bief_obj), INTENT(IN)      :: ACLADM
      TYPE(bief_obj), INTENT(INOUT)   :: CF,TOB
      TYPE(bief_obj), INTENT(INOUT)   :: FW,TOBW
      TYPE(bief_obj), INTENT(IN)      :: CF_TEL
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER                     :: I
      DOUBLE PRECISION            :: A,B,C, HCLIP,KSMAX
!
!     12.D0 WAS EXP(8.5*0.41)/EXP(1.D0)
!     11.036D0 IS EXP(8.5*0.40)/EXP(1.D0)
!     CONSIDERING THAT EXP(8.5*0.40) IS 30 INSTEAD OF 29.9641...
!     0.40 IS THE KARMAN CONSTANT THAT SHOULD BE PARAMETERISED SOMEWHERE
      DOUBLE PRECISION, PARAMETER :: CSTE=11.036d0
!
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! ----------------------------------------------------------------------------------------------
!  QUADRATIC FRICTION COEFICIENT       :  ---> CF
!-----------------------------------------------------------------------
!
!     INTERNAL COUPLING WITH TELEMAC2D OR 3D
!               UETCAR IS CF IN TELEMAC-2D
!               UETCAR IS UETCAR IN TELEMAC3D ?
!  KSP : skin friction
!  KSR: ripple roughness
!  KS : total bed roughness
!  initialisation
!
      CALL os('X=CY    ', x=ksp, y=acladm, c=kspratio)
      CALL os('X=CY    ', x=ksr, y=acladm, c=kspratio)
!
      IF(kspred) THEN
!
!       bed roughness predictor
!
        CALL ks_sisyphe(iks,ks,ksp,ksr,kspratio,houle,
     &                  grav,xmve,xmvs,vce,
     &                  hn,acladm,unorm,uw,tw,npoin)
        CALL coefro_sisyphe(cf,hn,kfrot,ks,grav,npoin,hmin,karman)
        IF(code(1:7).EQ.'TELEMAC')
     &   CALL ov( 'X=Y     ', ks_tel%R, ks%R, ks%R, 0.d0, npoin)
!
      ELSE
!
! here the total bed roughness is calculated as a function of friction coefficient
! -- > issued from Telemac if coupling
! -- > from the steering file of Sisyphe
!
        IF(code(1:7).EQ.'TELEMAC') THEN
          CALL ov('X=Y     ', x=cf%R, y=cf_tel%R, dim1=cf%DIM1)
        ELSE
          CALL coefro_sisyphe(cf,hn,kfrot,chestr,grav,npoin,hmin,karman)
        ENDIF
        DO i =1,npoin
          a = -karman*sqrt(2.d0/max(cf%R(i),zero))
          ks%R(i)=cste*hn%R(i)*exp(a)
          ks%R(i)=max(ks%R(i),ksp%R(i))
        ENDDO
!
      ENDIF
!
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! Frottement total: loi quadratique sauf couplage 3D
!  --> TOB
!
!     INTERNAL COUPLING WITH TELEMAC3D
!     UETCAR CORRESPONDS TO THE FRICTION VELOCITY SQUARED
!
      IF(code(1:9).EQ.'TELEMAC3D') THEN
        CALL os( 'X=CY     ',x=tob,y=uetcar,c=xmve)
      ELSE
        DO i=1,npoin
          tob%R(i) = xmve*0.5d0*cf%R(i)*unorm%R(i)**2
        ENDDO
      ENDIF
!
! -----WAVE-INDUCED FRICTION -----------------------------
!  --> TOBW
!
      IF(houle) THEN
        CALL tobw_sisyphe(tobw%R,cf%R,fw%R,uw%R,tw%R,hn%R,npoin,xmve)
      ENDIF
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! SKIN FRICTION CORRECTOR
!                ---> MU = TOP/TOB
! ICR=0:    MU=1
! ICR=1     : SKIN FRICTION CORRECTION USE KSP
! ICR= 2    : RIPPLE ROUGHNESS USE KSR, KSR
! COUPLED WITH TELEMAC: MU>1 IS ACCEPTABLE
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
      IF(icr.EQ.0) THEN
        CALL os('X=C     ', x=mu, c=1.d0)
      ELSEIF(icr.EQ.1) THEN
        DO i= 1, npoin
          IF(cf%R(i).GT.zero.AND.hn%R(i).GT.ksp%R(i)) THEN
            hclip=max(hn%R(i),ksp%R(i))
            a = 2.5d0*log(cste*hclip/ksp%R(i))
            c =2.d0/a**2
            mu%R(i) = c/cf%R(i)
          ELSE
            mu%R(i) = 0.d0
          ENDIF
        ENDDO
      ELSEIF(icr.EQ.2) THEN
        DO i= 1, npoin
          ksmax=max(ksr%R(i),ksp%R(i))
          IF(hn%R(i).GT.ksmax.AND.cf%R(i).GT.zero)THEN
            hclip=max(hn%R(i),ksmax)
            a = log(cste*hclip/ksp%R(i))
            b = log(cste*hclip/ksr%R(i))
            c = 0.32d0/cf%R(i)
            mu%R(i) = c/sqrt(b*a**3)
          ELSE
            mu%R(i) = 0.d0
          ENDIF
        ENDDO
      ENDIF
!
!------------------------------------------------------------
!
      RETURN
      END