calcfw.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\artemis\calcfw.f
00002 !
00058                      SUBROUTINE CALCFW
00059 !                    *****************
00060 !
00061      &(I,H,C,CG,K,HMU,
00062      & NPOIN,OMEGA,GRAV,
00063      & VISCO,DIAM90,DIAM50,MVSED,MVEAU,
00064      & FORMFR,REGIDO,RICOEF,
00065      & ENTREG,ENTRUG,FFW)
00066 !
00067 !***********************************************************************
00068 ! ARTEMIS   V6P1                                   21/08/2010
00069 !***********************************************************************
00070 !
00071 !
00072 !
00073 !
00074 !
00075 !
00076 !
00077 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00078 !| C              |-->| PHASE VELOCITY
00079 !| CG             |-->| GROUP VELOCITY
00080 !| DIAM50         |-->| MAXIMUM GRAIN DIAMETER WHICH DEFINES 50% OF THE
00081 !|                |   | TOTAL WEIGTH OF SEDIMENT
00082 !| DIAM90         |-->| MAXIMUM GRAIN DIAMETER WHICH DEFINES 90% OF THE
00083 !|                |   | TOTAL WEIGTH OF SEDIMENT
00084 !| ENTREG         |-->| LOGICAL USED TO IMPOSE THE HYDRAULIC REGIME
00085 !|                |   | (FOR FRICTION CALCULATION)
00086 !| ENTRUG         |-->| LOGICAL USED TO RESTRICT THE TOTAL ROUGHNESS TO
00087 !|                |   | THE SKIN ROUGHNESS (FOR FRICTION CALCULATION)
00088 !| FFW            |<--| FRICTION FACTOR
00089 !| FORMFR         |-->| FRICTION LAW
00090 !| GRAV           |-->| GRAVITY
00091 !| H              |-->| WATER HEIGHT
00092 !| HMU            |-->| WAVE HEIGHT
00093 !| I              |-->| CURRENT POINT NUMBER
00094 !| K              |-->| WAVE NUMBER
00095 !| MVEAU          |-->| FLUID SPECIFIC WEIGTH
00096 !| MVSED          |-->| SEDIMENT SPECIFIC WEIGTH
00097 !| NPOIN          |-->| NUMBER OF POINT
00098 !| OMEGA          |-->| WAVE  PULSATION
00099 !| REGIDO         |-->| TYPE OF HYDRAULIC REGIME
00100 !| RICOEF         |-->| RIPPLES COEFFICIENT
00101 !| VISCO          |-->| KINEMATIC VISCOSITY OF THE FLUID
00102 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00103 !
00104       USE BIEF
00105       USE INTERFACE_ARTEMIS, EX_CALCFW => CALCFW
00106 !
00107       IMPLICIT NONE
00108       INTEGER LNG,LU
00109       COMMON/INFO/LNG,LU
00110 !
00111       INTEGER NPOIN,I
00112       INTEGER REGIME, REGIDO
00113       INTEGER FORMFR
00114 !
00115       DOUBLE PRECISION C(NPOIN),CG(NPOIN)
00116       DOUBLE PRECISION K(NPOIN)
00117       DOUBLE PRECISION HMU(NPOIN)
00118 !
00119       DOUBLE PRECISION H(NPOIN)
00120       DOUBLE PRECISION GRAV,OMEGA
00121 !
00122       DOUBLE PRECISION VISCO, DIAM90, DIAM50
00123       DOUBLE PRECISION MVSED, MVEAU, RICOEF
00124       DOUBLE PRECISION KSRUGO, AEX, UEX, FFW
00125       DOUBLE PRECISION RAP1, RAP2
00126       DOUBLE PRECISION LL1
00127       DOUBLE PRECISION RMAX1, RMIN1, RMAX2, RMIN2
00128       DOUBLE PRECISION SP, PSI, KS1, RH, RS, KS2
00129       DOUBLE PRECISION DIAMAD, TETACR, TAUCR
00130 !
00131       LOGICAL ENTREG,ENTRUG
00132 !
00133 !
00134       INTRINSIC ABS,EXP
00135 !
00136 !
00137 !  COMPUTES FW OVER THE WHOLE DOMAIN
00138 !--------------------------------------------------------
00139       RMAX1 = 0.D0
00140       RMIN1 = 1.D7
00141       RMAX2 = 0.D0
00142       RMIN2 = 1.D7
00143 !
00144       KSRUGO = 3.D0 * DIAM90
00145 !
00146 ! COMPUTES AEX AND UEX AT EACH GRID NODE
00147 ! THUS THE HYDRAULIC REGIME WILL ALSO BE GIVEN AT EACH
00148 ! GRID NODE
00149 !
00150 !
00151       AEX = HMU(I)/(2.D0*SINH(K(I)*H(I)))
00152       UEX = OMEGA * AEX
00153       RAP1 = (UEX*AEX)/VISCO
00154       RAP2 = AEX/KSRUGO
00155 !
00156 !   COMPARES AGAINST THE MAXIMUM REYNOLDS NUMBER
00157 !   AND THE EXCURSION RATIO ON ROUGHNESS
00158 !
00159       IF (RAP1 .GT. RMAX1) THEN
00160         RMAX1 = RAP1
00161       ENDIF
00162 !
00163       IF (RAP1 .LT. RMIN1) THEN
00164         RMIN1 = RAP1
00165       ENDIF
00166 !
00167       IF (RAP2 .GT. RMAX2) THEN
00168         RMAX2 = RAP2
00169       ENDIF
00170 !
00171       IF (RAP2 .LT. RMIN2) THEN
00172         RMIN2 = RAP2
00173       ENDIF
00174 !
00175 !-----------------------------------------------------------------
00176 ! DETERMINES THE HYDRAULIC REGIME
00177 !-----------------------------------------------------------------
00178 !
00179       IF (ENTREG) THEN
00180         REGIME = REGIDO
00181       ELSE
00182 !
00183 !     INITIALIZES THE REGIME (SETS TO 0) BEFORE EACH NEW ITERATION
00184 !
00185         REGIME = 0
00186         LL1 = 0.D0
00187 !
00188         IF (RAP1 .LE. 10250.D0) THEN
00189           LL1 = 0.0322D0*RAP1+3.33D0
00190           IF (RAP2 .GE. LL1) THEN
00191             REGIME = 1
00192 !             WRITE(*,*) 'LE REGIME HYDRAULIQUE EST LAMINAIRE'
00193 !             WRITE(*,*) 'THE HYDRAULIC REGIME IS LAMINAR'
00194           ENDIF
00195         ENDIF
00196 !
00197         IF (RAP1 .GE. 3.D4) THEN
00198           LL1 = 0.009792D0*RAP1+208.33D0
00199           IF (RAP2 .GE. LL1) THEN
00200             REGIME = 2
00201 !             WRITE(*,*) 'LE REGIME HYDRAULIQUE EST TURBULENT LISSE'
00202 !             WRITE(*,*) 'THE HYDRAULIC REGIME IS SMOOTH TURBULENT'
00203           ENDIF
00204         ENDIF
00205 !
00206         IF (RAP1 .GE. 5.D3 .AND. RAP1 .LE. 2.D4) THEN
00207           LL1 = 0.026D0*RAP1-12.D0
00208           IF (RAP2 .LE. LL1) THEN
00209             REGIME = 3
00210 !             WRITE(*,*) 'LE REGIME HYDRAULIQUE EST TURBULENT RUGUEUX'
00211 !             WRITE(*,*) 'THE HYDRAULIC REGIME IS ROUGH TURBULENT'
00212           ENDIF
00213         ENDIF
00214 !
00215         IF (RAP1 .GT. 2.D4) THEN
00216           LL1 = 0.00099D0*RAP1+30.30D0
00217           IF (RAP2 .LE. LL1) THEN
00218             REGIME = 3
00219 !             WRITE(*,*) 'LE REGIME HYDRAULIQUE EST TURBULENT RUGUEUX'
00220 !             WRITE(*,*) 'THE HYDRAULIC REGIME IS ROUGH TURBULENT'
00221           ENDIF
00222         ENDIF
00223 !
00224         IF (REGIME .EQ. 0) THEN
00225           REGIME = 3
00226 !           WRITE(*,*) 'LE REGIME EST TURBULENT RUGEUX'
00227 !             WRITE(*,*) 'THE HYDRAULIC REGIME IS ROUGH TURBULENT'
00228         ENDIF
00229 !
00230       ENDIF
00231 !
00232 !
00233 !--------------------------------------------------------------
00234 ! COMPUTES FW: THE COEFFICIENT OF FRICTION
00235 !--------------------------------------------------------------
00236 !
00237       IF (REGIME .EQ. 1) THEN
00238         FFW = 2.D0*(((AEX*UEX)/VISCO)**(-0.5D0))
00239       ENDIF
00240 !
00241       IF (REGIME .EQ. 2) THEN
00242         FFW = 0.09*(((AEX*UEX)/VISCO)**(-0.2D0))
00243       ENDIF
00244 !
00245       IF (REGIME .EQ. 3) THEN
00246 !
00247         IF (ENTRUG) THEN
00248 !
00249           KSRUGO = 3.D0 * DIAM90
00250 !
00251         ELSE
00252 !
00253           SP = MVSED/MVEAU
00254           PSI = (UEX**2.D0)/((SP-1.D0)*GRAV*DIAM50)
00255 !
00256           IF (PSI .LT. 250.D0) THEN
00257             KS1 = 3.D0*DIAM90
00258 !
00259             IF (PSI .LT. 10.D0) THEN
00260               RH = 0.22D0*AEX
00261               RS = 0.18D0
00262             ELSE
00263               RH = AEX*(2.4D-13)*((250.D0-PSI)**5.D0)
00264               RS = (2.D-7)*((250.D0-PSI)**2.5D0)
00265             ENDIF
00266 !
00267             KS2 = 20.D0*RICOEF*RH*RS
00268 !
00269 !  RICOEF = 1     RIPPLES ONLY
00270 !  RICOEF = 0.7D0 RIPPLES AND SAND WAVES
00271 !
00272             ELSE
00273               KS1 = 3.D0*(0.04D0*PSI-9.D0)*DIAM90
00274               KS2 = 0.D0
00275               IF (KS1 .LT. 0.01D0) THEN
00276                 KS1 = 3.D0*DIAM90
00277               ENDIF
00278             ENDIF
00279 !
00280             KSRUGO = KS1+KS2
00281             KS1 = (AEX/KSRUGO)
00282 !
00283           ENDIF
00284 !
00285           FFW = EXP(-6.D0+5.2D0*((AEX/KSRUGO)**(-0.19D0)))
00286           IF (FFW .GT. 0.3D0) THEN
00287             FFW = 0.3D0
00288           ENDIF
00289         ENDIF
00290 !
00291         IF (REGIME .EQ. 4) THEN
00292 !
00293 !     TRANSIENT STATE NOT TAKEN INTO ACCOUNT
00294 !
00295           IF (ENTRUG) THEN
00296 !
00297             KSRUGO = 3.D0 * DIAM90
00298 !
00299           ELSE
00300 !
00301             DIAMAD = (((MVSED/MVEAU)-1.D0)*GRAV)/(VISCO**2.D0)
00302             DIAMAD = DIAMAD**(1/3)
00303             DIAMAD = DIAMAD * DIAM50
00304 !
00305             TETACR = 0.14D0*(DIAMAD**(-0.64D0))
00306 !
00307             TAUCR = TETACR*(MVSED-MVEAU)*GRAV*DIAM50
00308 !
00309             KSRUGO = (3.D0*DIAM90)+((3.3D0*VISCO)/
00310      &               ((TAUCR/MVEAU)**0.5D0))
00311 !
00312           ENDIF
00313 !
00314           FFW = EXP(-6.D0+5.2D0*((AEX/KSRUGO)**(-0.19D0)))
00315           IF (FFW .GT. 0.3D0) THEN
00316             FFW = 0.3D0
00317           ENDIF
00318 !
00319         ENDIF
00320 !
00321       RETURN
00322       END

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