bedload_dibwat.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\sisyphe\bedload_dibwat.f
00002 !
00069                      SUBROUTINE BEDLOAD_DIBWAT !
00070 !                    ***************************
00071 !
00072      &(U2D,V2D,UCMOY, CF, TOB, TOBW, UW, TW, FW, THETAW, NPOIN,
00073      & XMVE, DENS, GRAV, DM, XWC, PI, ALPHAW, T2, T3, UCW, UCN,
00074      & UW1, UW2, TW1, TW2, THETAC, FCW, QSC,HOULE)
00075 !
00076 !***********************************************************************
00077 ! SISYPHE   V6P1                                   21/07/2011
00078 !***********************************************************************
00079 !
00080 !
00081 !
00082 !
00083 !
00084 !
00085 !
00086 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00087 !| ALPHAW         |<->| ANGLE BETWEEN WAVE AND CURRENT
00088 !| CF             |-->| QUADRATIC FRICTION COEFFICIENT
00089 !| DENS           |-->| RELATIVE DENSITY
00090 !| DM             |-->| SEDIMENT GRAIN DIAMETER
00091 !| FCW            |-->| WAVE-CURRENT FRICTION FACTOR
00092 !| FW             |-->| WAVE FRICTION FACTOR
00093 !| GRAV           |-->| ACCELERATION OF GRAVITY
00094 !| HOULE          |-->| LOGICAL, FOR WAVE EFFECTS
00095 !| NPOIN          |-->| NUMBER OF POINTS
00096 !| PI             |-->| PI
00097 !| QSC            |<->| BED LOAD TRANSPORT
00098 !| T2             |<->| WORK BIEF_OBJ STRUCTURE
00099 !| T3             |<->| WORK BIEF_OBJ STRUCTURE
00100 !| THETAC         |<->| CURRENT ANGLE TO THE X AXIS
00101 !| THETAW         |-->| ANGLE BETWEEN WAVE AND CURRENT
00102 !| TOB            |-->| BED SHEAR STRESS (TOTAL FRICTION)
00103 !| TOBW           |-->| WAVE INDUCED SHEAR STRESS
00104 !| TW             |-->| WAVE PERIOD
00105 !| TW1            |<->| MID PERIOD (U(T)>0)
00106 !| TW2            |<->| MID PERIOD (U(T)<0)
00107 !| U2D            |<->| MEAN FLOW VELOCITY X-DIRECTION
00108 !| UCMOY          |-->| MEAN CURRENT
00109 !| UCN            |<->| MEAN CURRENT AT AN ANGLE TO THE WAVE
00110 !| UCW            |<->| MEAN CURRENT PROJECTED IN THE WAVE DIRECTION
00111 !| UW             |-->| ORBITAL WAVE VELOCITY
00112 !| UW1            |<->| MEAN CURRENT IN THE WAVE DIRECTION
00113 !| UW2            |<->| MEAN CURRENT IN THE OPPOSITE DIRECTION
00114 !| V2D            |<->| MEAN FLOW VELOCITY Y-DIRECTION
00115 !| XMVE           |-->| FLUID DENSITY
00116 !| XWC            |-->| SETTLING VELOCITY
00117 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00118 !
00119       USE INTERFACE_SISYPHE,EX_BEDLOAD_DIBWAT => BEDLOAD_DIBWAT
00120       USE BIEF
00121       IMPLICIT NONE
00122       INTEGER LNG,LU
00123       COMMON/INFO/LNG,LU
00124 !
00125       ! 2/ GLOBAL VARIABLES
00126       ! -------------------
00127       TYPE(BIEF_OBJ),   INTENT(IN)    :: U2D,V2D,UCMOY, CF, TOB, TOBW
00128       TYPE(BIEF_OBJ),   INTENT(IN)    :: UW, TW, FW
00129       TYPE(BIEF_OBJ),   INTENT(IN)    :: THETAW
00130       INTEGER,          INTENT(IN)    :: NPOIN
00131       LOGICAL,          INTENT(IN)    :: HOULE
00132       DOUBLE PRECISION, INTENT(IN)    :: XMVE, DENS, GRAV, DM, XWC, PI
00133       TYPE(BIEF_OBJ),   INTENT(INOUT) :: ALPHAW          ! WORK ARRAY T1
00134       TYPE(BIEF_OBJ),   INTENT(INOUT) :: T2, T3          !
00135       TYPE(BIEF_OBJ),   INTENT(INOUT) ::  UCW, UCN ! WORK ARRAY T4, T5, T6
00136       TYPE(BIEF_OBJ),   INTENT(INOUT) :: UW1, UW2, TW1   ! WORK ARRAY T7, T8, T9
00137       TYPE(BIEF_OBJ),   INTENT(INOUT) :: TW2, THETAC     ! WORK ARRAY T10, T11
00138       TYPE(BIEF_OBJ),   INTENT(INOUT) :: FCW, QSC        ! WORK ARRAY T12
00139 !
00140       ! 3/ LOCAL VARIABLES
00141       ! ------------------
00142       INTEGER                     :: I
00143       DOUBLE PRECISION            :: UCW2, SHIELDS, OMEGACR, OMEGA1
00144       DOUBLE PRECISION            :: OMEGA2, W1, WP1, W2, WP2, GAMMAW
00145       DOUBLE PRECISION            :: GAMMAN, GAMMA, QSW, QSN
00146       DOUBLE PRECISION, PARAMETER :: ZERO  = 1.D-6
00147       DOUBLE PRECISION, PARAMETER :: ALPHA = 0.0001D0
00148       DOUBLE PRECISION, PARAMETER :: BETA  = 0.55D0
00149 !
00150 !======================================================================!
00151 !======================================================================!
00152 !                               PROGRAM                                !
00153 !======================================================================!
00154 !======================================================================!
00155 !
00156 !     ANGLE OF VELOCITY WITH OX
00157 !
00158       CALL BEDLOAD_DIRECTION(U2D,V2D,NPOIN,PI,THETAC)
00159 !
00160 !     PROJECTIONS OF UCMOY IN THE WAVE DIRECTION
00161 !
00162       CALL OS('X=CY    ', X=ALPHAW, Y=THETAW, C=-PI/180.D0)
00163       CALL OS('X=X+C   ', X=ALPHAW, C=0.5D0*PI)
00164       CALL OS('X=Y-Z   ', X=ALPHAW, Y=ALPHAW, Z=THETAC)
00165       CALL OS('X=COS(Y)', X=T2    , Y=ALPHAW)
00166       CALL OS('X=SIN(Y)', X=T3    , Y=ALPHAW)
00167       CALL OS('X=YZ    ', X=UCW   , Y=UCMOY , Z=T2)
00168       CALL OS('X=CYZ   ', X=UCN   , Y=UCMOY , Z=T3, C= -1.D0)
00169 !
00170 !     QUADRATIC VELOCITIES (UW1=U(T)> 0 AND UW2=U(T)
00171 !
00172       CALL BEDLOAD_CALCDW(UCW,UW,TW,NPOIN,PI,UW1,UW2,TW1,TW2)
00173 !
00174 !     FRICTION COEFFICIENT
00175 !
00176       CALL BEDLOAD_INTERACT
00177      &     (UCMOY,TOBW,TOB,ALPHAW,FW,CF,UW,NPOIN,XMVE,FCW)
00178 !
00179       DO I = 1, NPOIN
00180 !
00181 !       SHIELDS PARAMETER
00182         UCW2    = UCMOY%R(I)**2 + 0.5D0 * UW%R(I)**2
00183         SHIELDS = FCW%R(I)*UCW2/DENS/GRAV/DM
00184 !
00185 !       CRITICAL OMEGA (TEMPERVILLE AND GUIZA,2000)  !
00186         IF (SHIELDS <= 0.2D0) THEN
00187           OMEGACR = 0.03D0
00188         ELSEIF (SHIELDS <= 0.4D0) THEN
00189           OMEGACR = 1.D0 - SQRT(1.D0-((SHIELDS-0.2D0)/0.58D0)**2)
00190         ELSEIF (SHIELDS <= 1.5D0) THEN
00191           OMEGACR = 0.8104D0 * SQRT(SHIELDS) - 0.4225D0
00192         ELSE
00193           OMEGACR = 0.7236D0 * SQRT(SHIELDS) - 0.3162D0
00194         ENDIF
00195 !
00196 !       OMEGA1, OMEGA2  !
00197         IF(TW1%R(I) > ZERO) THEN
00198           OMEGA1 = UW1%R(I)**2 / (2.D0*DENS*GRAV*XWC*TW1%R(I))
00199         ELSE
00200           OMEGA1 = ZERO
00201         ENDIF
00202 !
00203         IF(TW2%R(I) > ZERO) THEN
00204           OMEGA2 = UW2%R(I)**2 / (2.D0*DENS*GRAV*XWC*TW2%R(I))
00205         ELSE
00206           OMEGA2 = ZERO
00207         ENDIF
00208 !
00209 !       QUANTITIES OF SAND DEPOSITED AND
00210 !       IN SUSPENSION AT EACH PHASE OF THE CYCLE
00211 !
00212 !       PHASE 1
00213         IF (OMEGA1 <= OMEGACR) THEN
00214           W1  = 2.D0 * OMEGA1 * XWC * TW1%R(I) / DM
00215           WP1 = 0.D0
00216         ELSE
00217           W1  = 2.D0 * OMEGACR * XWC * TW1%R(I) / DM
00218           WP1 = 2.D0 * (OMEGA1-OMEGACR) * XWC * TW1%R(I) / DM
00219         ENDIF
00220 !
00221 !       PHASE 2
00222         IF (OMEGA2 <= OMEGACR) THEN
00223           W2  = 2.D0 * OMEGA2 * XWC * TW2%R(I) / DM
00224           WP2 = 0.D0
00225         ELSE
00226           W2  = 2.D0 * OMEGACR * XWC * TW2%R(I) / DM
00227           WP2 = 2.D0 * (OMEGA2-OMEGACR) * XWC * TW2%R(I) / DM
00228         ENDIF
00229 !
00230 !       GAMMAW, GAMMAN,GAMMA
00231         IF ((UW2%R(I)*TW2%R(I) + UW1%R(I)*TW1%R(I)) > 0.D0 ) THEN
00232           GAMMAW = (  UW1%R(I) * TW1%R(I) * (W1**3+WP2**3)
00233      &              - UW2%R(I) * TW2%R(I) * (W2**3+WP1**3))
00234      &           / ( UW1%R(I)*TW1%R(I)+UW2%R(I)*TW2%R(I))
00235         ELSE
00236           GAMMAW = (2.D0 * UCW%R(I)**2 / DENS / GRAV / DM)**3
00237         ENDIF
00238         GAMMAN = (2.D0 * UCN%R(I)**2 / DENS / GRAV / DM)**3
00239         GAMMA  = MAX(SQRT(GAMMAW**2 + GAMMAN**2),1.D-10)
00240 !
00241 !       SOLID TRANSPORT IN THE WAVE DIRECTION : QSW
00242 !       AND IN THE NORMAL DIRECTION : QSN
00243 !
00244         QSW      = ALPHA * GAMMAW * XWC * DM / GAMMA**(1.D0-BETA)
00245         QSN      = ALPHA * GAMMAN * XWC * DM / GAMMA**(1.D0-BETA)
00246         QSC%R(I) = SQRT(QSW**2 + QSN**2)
00247 !
00248       ENDDO
00249 !
00250 !----------------------------------------------------------------------
00251 !
00252       RETURN
00253       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0