The TELEMAC-MASCARET system  trunk
bedload_dibwat_gaia.f
Go to the documentation of this file.
1 ! ******************************
2  SUBROUTINE bedload_dibwat_gaia
3 ! ******************************
4 !
5  &(u2d,v2d,unorm, cf, tob, tobw, uw, tw, fw, thetaw, npoin,
6  & xmve, dens, grav, dcla, xwc, pi, alphaw, t2, t3, ucw, ucn,
7  & uw1, uw2, tw1, tw2, thetac, fcw, qsc,xmvs)
8 !
9 !***********************************************************************
10 ! GAIA
11 !***********************************************************************
12 !
14 !
15 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
46 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
47 !
48  USE interface_gaia,ex_bedload_dibwat => bedload_dibwat_gaia
49  USE bief
51  IMPLICIT NONE
52 !
53 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
54 !
55  TYPE(bief_obj), INTENT(IN) :: U2D,V2D,UNORM, CF, TOB, TOBW
56  TYPE(bief_obj), INTENT(IN) :: UW, TW, FW
57  TYPE(bief_obj), INTENT(IN) :: THETAW
58  INTEGER, INTENT(IN) :: NPOIN
59  DOUBLE PRECISION, INTENT(IN) :: XMVE, DENS, GRAV, DCLA, XWC, PI
60  DOUBLE PRECISION, INTENT(IN) :: XMVS
61  TYPE(bief_obj), INTENT(INOUT) :: ALPHAW ! WORK ARRAY T1
62  TYPE(bief_obj), INTENT(INOUT) :: T2, T3 !
63  TYPE(bief_obj), INTENT(INOUT) :: UCW, UCN ! WORK ARRAY T4, T5, T6
64  TYPE(bief_obj), INTENT(INOUT) :: UW1, UW2, TW1 ! WORK ARRAY T7, T8, T9
65  TYPE(bief_obj), INTENT(INOUT) :: TW2, THETAC ! WORK ARRAY T10, T11
66  TYPE(bief_obj), INTENT(INOUT) :: FCW, QSC ! WORK ARRAY T12
67 !
68 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
69 !
70 ! LOCAL VARIABLES
71 ! ------------------
72  INTEGER :: I
73  DOUBLE PRECISION :: UCW2, SHIELDS, OMEGACR, OMEGA1
74  DOUBLE PRECISION :: OMEGA2, W1, WP1, W2, WP2, GAMMAW
75  DOUBLE PRECISION :: GAMMAN, GAMMA, QSW, QSN
76  DOUBLE PRECISION, PARAMETER :: ZERO = 1.d-6
77  DOUBLE PRECISION, PARAMETER :: ALPHA = 0.0001d0
78  DOUBLE PRECISION, PARAMETER :: BETA = 0.55d0
79 !
80 !======================================================================!
81 !======================================================================!
82 ! PROGRAM !
83 !======================================================================!
84 !======================================================================!
85 !
86 ! ANGLE OF VELOCITY WITH OX
87 !
88  CALL bedload_direction_gaia(u2d,v2d,npoin,pi,thetac)
89 !
90 ! PROJECTIONS OF UNORM IN THE WAVE DIRECTION
91 !
92  CALL os('X=CY ', x=alphaw, y=thetaw, c=-pi/180.d0)
93  CALL os('X=X+C ', x=alphaw, c=0.5d0*pi)
94  CALL os('X=Y-Z ', x=alphaw, y=alphaw, z=thetac)
95  CALL os('X=COS(Y)', x=t2 , y=alphaw)
96  CALL os('X=SIN(Y)', x=t3 , y=alphaw)
97  CALL os('X=YZ ', x=ucw , y=unorm , z=t2)
98  CALL os('X=CYZ ', x=ucn , y=unorm , z=t3, c= -1.d0)
99 !
100 ! QUADRATIC VELOCITIES (UW1=U(T)> 0 AND UW2=U(T)
101 !
102  CALL bedload_calcdw_gaia(ucw,uw,tw,npoin,pi,uw1,uw2,tw1,tw2)
103 !
104 ! FRICTION COEFFICIENT
105 !
107  & (unorm,tobw,tob,alphaw,fw,cf,uw,npoin,xmve,fcw)
108 !
109  DO i = 1, npoin
110 !
111 ! SHIELDS PARAMETER
112  ucw2 = unorm%R(i)**2 + 0.5d0 * uw%R(i)**2
113  shields = fcw%R(i)*ucw2/dens/grav/dcla
114 !
115 ! CRITICAL OMEGA (TEMPERVILLE AND GUIZA,2000) !
116  IF (shields <= 0.2d0) THEN
117  omegacr = 0.03d0
118  ELSEIF (shields <= 0.4d0) THEN
119  omegacr = 1.d0 - sqrt(1.d0-((shields-0.2d0)/0.58d0)**2)
120  ELSEIF (shields <= 1.5d0) THEN
121  omegacr = 0.8104d0 * sqrt(shields) - 0.4225d0
122  ELSE
123  omegacr = 0.7236d0 * sqrt(shields) - 0.3162d0
124  ENDIF
125 !
126 ! OMEGA1, OMEGA2 !
127  IF(tw1%R(i) > zero) THEN
128  omega1 = uw1%R(i)**2 / (2.d0*dens*grav*xwc*tw1%R(i))
129  ELSE
130  omega1 = zero
131  ENDIF
132 !
133  IF(tw2%R(i) > zero) THEN
134  omega2 = uw2%R(i)**2 / (2.d0*dens*grav*xwc*tw2%R(i))
135  ELSE
136  omega2 = zero
137  ENDIF
138 !
139 ! QUANTITIES OF SAND DEPOSITED AND
140 ! IN SUSPENSION AT EACH PHASE OF THE CYCLE
141 !
142 ! PHASE 1
143  IF (omega1 <= omegacr) THEN
144  w1 = 2.d0 * omega1 * xwc * tw1%R(i) / dcla
145  wp1 = 0.d0
146  ELSE
147  w1 = 2.d0 * omegacr * xwc * tw1%R(i) / dcla
148  wp1 = 2.d0 * (omega1-omegacr) * xwc * tw1%R(i) / dcla
149  ENDIF
150 !
151 ! PHASE 2
152  IF (omega2 <= omegacr) THEN
153  w2 = 2.d0 * omega2 * xwc * tw2%R(i) / dcla
154  wp2 = 0.d0
155  ELSE
156  w2 = 2.d0 * omegacr * xwc * tw2%R(i) / dcla
157  wp2 = 2.d0 * (omega2-omegacr) * xwc * tw2%R(i) / dcla
158  ENDIF
159 !
160 ! GAMMAW, GAMMAN,GAMMA
161  IF ((uw2%R(i)*tw2%R(i) + uw1%R(i)*tw1%R(i)) > 0.d0 ) THEN
162  gammaw = ( uw1%R(i) * tw1%R(i) * (w1**3+wp2**3)
163  & - uw2%R(i) * tw2%R(i) * (w2**3+wp1**3))
164  & / ( uw1%R(i)*tw1%R(i)+uw2%R(i)*tw2%R(i))
165  ELSE
166  gammaw = (2.d0 * ucw%R(i)**2 / dens / grav / dcla)**3
167  ENDIF
168  gamman = (2.d0 * ucn%R(i)**2 / dens / grav / dcla)**3
169  gamma = max(sqrt(gammaw**2 + gamman**2),1.d-10)
170 !
171 ! SOLID TRANSPORT IN THE WAVE DIRECTION : QSW
172 ! AND IN THE NORMAL DIRECTION : QSN
173 !
174  qsw = alpha * gammaw * xwc * dcla / gamma**(1.d0-beta)
175  qsn = alpha * gamman * xwc * dcla / gamma**(1.d0-beta)
176  qsc%R(i) = sqrt(qsw**2 + qsn**2)
177 !
178  ENDDO
179 !
180 ! SOLID DISCHARGE IS TRANSFORMED IN [kg/(m*s)]
181 !
182  CALL os('X=CX ', x=qsc, c=xmvs)
183 !----------------------------------------------------------------------
184 !
185  RETURN
186  END
subroutine bedload_interact_gaia(UCMOY, TOBW, TOB, ALPHAW, FW, CF, UW, NPOIN, XMVE, FCW)
subroutine bedload_dibwat_gaia(U2D, V2D, UNORM, CF, TOB, TOBW, UW, TW, FW, THETAW, NPOIN, XMVE, DENS, GRAV, DCLA, XWC, PI, ALPHAW, T2, T3, UCW, UCN, UW1, UW2, TW1, TW2, THETAC, FCW, QSC, XMVS)
subroutine bedload_direction_gaia(U2D, V2D, NPOIN, PI, THETAC)
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
subroutine bedload_calcdw_gaia(UCW, UW, TW, NPOIN, PI, UW1, UW2, TW1, TW2)
Definition: bief.f:3