suspension_dispersion.f

Go to the documentation of this file.

!                   ******************************************
                    SUBROUTINE suspension_dispersion ! (_IMP_)
!                   ******************************************
!
     &  (tob, xmve,hn,  optdif, npoin, xkx, xky,
     &   t1, t2, t3, kx, ky, kz, disp,u2d,v2d,visc_tel,code)
!
!***********************************************************************
! SISYPHE   V6P1                                   21/07/2011
!***********************************************************************
!
!brief    COMPUTES THE DISPERSION PARAMETERS.
!
!history  C. MOULIN (LNH)
!+        13/12/2000
!+        V5P1
!+
!
!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
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CODE           |-->| HYDRODYNAMIC CODE IN CASE OF COUPLING
!| DISP           |-->| VISCOSITY COEFFICIENTS ALONG X,Y AND Z
!|                |   | IF P0 : PER ELEMENT
!|                |   | IF P1 : PER POINT
!| HN             |-->| WATER DEPTH
!| KX             |<->| COEFFICIENTS OF THE DISPERSION TENSOR (DIM. NPOIN)
!| KY             |<->| COEFFICIENTS OF THE DISPERSION TENSOR (DIM. NPOIN)
!| KZ             |<->| COEFFICIENTS OF THE DISPERSION TENSOR (DIM. NPOIN)
!| NPOIN          |-->| NUMBER OF POINTS
!| OPTDIF         |-->| OPTION FOR THE DISPERSION
!| T1             |<->| WORK BIEF_OBJ STRUCTURE
!| T2             |<->| WORK BIEF_OBJ STRUCTURE
!| T3             |<->| WORK BIEF_OBJ STRUCTURE
!| TOB            |-->| BED SHEAR STRESS (TOTAL FRICTION)
!| U2D            |-->| MEAN FLOW VELOCITY X-DIRECTION
!| V2D            |-->| MEAN FLOW VELOCITY Y-DIRECTION
!| VISC_TEL       |-->| VELOCITY DIFFUSIVITY (TELEMAC)
!| XKX            |-->| COEFFICIENT USED FOR COMPUTING THE DISPERSION
!|                |   | DEPENDS OF OPTIONS
!| XKY            |-->| COEFFICIENT USED FOR COMPUTING THE DISPERSION
!|                |   | DEPENDS OF OPTIONS
!| XMVE           |-->| FLUID DENSITY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe,
     &    ex_suspension_dispersion => suspension_dispersion
      USE bief
      USE declarations_special
      IMPLICIT NONE
      ! 2/ GLOBAL VARIABLES
      ! -------------------
      type(bief_obj),  INTENT(IN)    :: tob,hn,visc_tel
      INTEGER,          INTENT(IN)    :: OPTDIF, NPOIN
      DOUBLE PRECISION, INTENT(IN)    :: XMVE, XKX, XKY
      type(bief_obj),  INTENT(INOUT) :: t1, t2, t3
      type(bief_obj),  INTENT(INOUT) :: kx, ky, kz, disp
      type(bief_obj),  INTENT(IN)    :: u2d,v2d
      CHARACTER(LEN=24), INTENT(IN)   :: CODE
      ! 3/ LOCAL VARIABLES
      ! ------------------
      INTEGER                     :: K,DIMVISC
      DOUBLE PRECISION            :: UETH, COST, SINT
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
      ! ****************************************************** !
      ! IA - CONSTANT DISPERSION OR DISPERSION(I) = ALPHA(I)*H !
      ! ****************************************************** !
      IF (optdif == 2.OR.optdif == 1) THEN
!
        CALL cpstvc(u2d,t1)
        CALL cpstvc(u2d,t2)
!
        IF(optdif == 2) THEN
          DO k = 1, npoin
            ueth = sqrt(tob%R(k)/xmve)
            t1%R(k) = xkx * ueth * hn%R(k)
            t2%R(k) = xky * ueth * hn%R(k)
          ENDDO
        ELSE
          CALL os('X=C     ', x=t1, c=xkx)
          CALL os('X=C     ', x=t2, c=xky)
        ENDIF
!
        CALL os('X=N(Y,Z)', x=t3, y=u2d, z=v2d)
!
        DO k=1,npoin
!
          IF(t3%R(k).GE.1.d-6) THEN
            cost = u2d%R(k)/t3%R(k)
            sint = v2d%R(k)/t3%R(k)
          ELSE
            cost = 0.d0
            sint = 0.d0
          ENDIF
!
          kx%R(k) = (t1%R(k) - t2%R(k))*(cost**2) + t2%R(k)
          ky%R(k) = (t2%R(k) - t1%R(k))*(cost**2) + t1%R(k)
          kz%R(k) = (t1%R(k) - t2%R(k))*cost*sint
!
        ENDDO
!
!
      ! *********************************** !
      ! IB - DISPERSION GIVEN BY TELEMAC-2D ! (_IMP_)
      ! *********************************** !
      ELSEIF(optdif == 3) THEN
!
        IF(code(1:9).EQ.'TELEMAC2D') THEN
          IF(visc_tel%DIM2.EQ.1) THEN
            CALL os('X=Y     ', x=kx,y=visc_tel)
            CALL os('X=Y     ', x=ky,y=kx)
            CALL os('X=0     ', x=kz)
          ELSEIF(visc_tel%DIM2.EQ.3) THEN
            dimvisc=visc_tel%MAXDIM1
            DO k=1,npoin
              kx%R(k)=visc_tel%R(k)
              ky%R(k)=visc_tel%R(k+  dimvisc)
              kz%R(k)=visc_tel%R(k+2*dimvisc)
            ENDDO
          ELSE
            WRITE(lu,*) 'SUSPENSION_DISPERSION:'
            WRITE(lu,*) ' '
            WRITE(lu,*) 'UNEXPECTED DIMENSION OF VISC_TEL:',
     &                   visc_tel%DIM2
            CALL plante(1)
            stop
          ENDIF
        ELSE
          WRITE(lu,*) ' '
          WRITE(lu,*) 'SUSPENSION_DISPERSION:'
          WRITE(lu,*) ' '
          WRITE(lu,*) 'OPTION 3: DIFFUSIVITY GIVEN BY TELEMAC'
          WRITE(lu,*) 'NOT IMPLEMENTED OR IMPOSSIBLE WITH ',code
          CALL plante(1)
          stop
        ENDIF
!
      ! ***************************************** !
      ! IC - OPTION FOR DISPERSION NOT CODED UP   ! (_IMP_)
      ! ***************************************** !
      ELSE
        WRITE(lu,31) optdif
        CALL plante(1)
        stop
      ENDIF
!
      CALL ov_2('X=Y     ', disp%R, 1, kx%R, 1, kx%R, 1, 0.d0,
     &          disp%MAXDIM1, disp%DIM1)
      CALL ov_2('X=Y     ', disp%R, 2, ky%R, 1, ky%R, 1, 0.d0,
     &          disp%MAXDIM1, disp%DIM1)
      CALL ov_2('X=Y     ', disp%R, 3, kz%R, 1, kz%R, 1, 0.d0,
     &          disp%MAXDIM1, disp%DIM1)
      !----------------------------------------------------------------!
31    FORMAT('DISPERSION: OPTION FOR THE DISPERSION NOT AVAILABLE:',1i6)
      !----------------------------------------------------------------!
!======================================================================!
!======================================================================!
      RETURN
      END SUBROUTINE suspension_dispersion