calcmu.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\artemis\calcmu.f
00002 !
00097                      SUBROUTINE CALCMU
00098 !                    *****************
00099      &(ITERMU)
00100 !
00101 !
00102 !***********************************************************************
00103 ! ARTEMIS   V7P0                                   06/2014
00104 !***********************************************************************
00105 !
00106 !
00107 !
00108 !
00109 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00110 !| ITERMU             |-->| INDICE OF THE CURRENT CALCULATION
00111 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00112 !
00113       USE BIEF
00114       USE DECLARATIONS_TELEMAC
00115       USE DECLARATIONS_ARTEMIS
00116       IMPLICIT NONE
00117 
00118       INTEGER I , ITERMU
00119       DOUBLE PRECISION CBID,HEFF,Q1,Q2,Q3,HM,FFW,HMUE
00120       DOUBLE PRECISION PI
00121       PARAMETER(PI = 3.1415926535897932384626433D0)
00122 !-----------------------------------------------------------------------
00123 !
00124       INTRINSIC ABS,MIN,MAX,LOG
00125       DOUBLE PRECISION P_DMAX
00126       EXTERNAL P_DMAX
00127 !
00128 !----------------------------------------------------------------------
00129 
00130 
00131 !------------------------------------------------------------------------------
00132 ! K   = MEAN WAVE NUMBER
00133 ! C   = MEAN PHASE CELERITY
00134 ! CG  = MEAN GROUP CELERITY
00135 ! HALE= SIGNIFICANT WAVE HEIGHT
00136 !
00137 !
00138 !     ============================================================
00139 !
00140 !     COMPUTES THE TOTAL DISSIPATION COEFFICIENT MU_DEFERL + MU_FROTTE
00141 !                                                  (MU2)       (T1)
00142 !     IF BREAKING OR BOTTOM FRICTION TAKEN INTO ACCOUNT
00143 !     ============================================================
00144 !
00145 !
00146 !      ECRHMU = 0.D0
00147 !
00148 !     --------------------------------------------
00149 !     INITIALISES MU2 AND T3: SET TO 0
00150 !     MU2: NEW DISSIPATION COEFFICIENT
00151 !     T3: QB FOR THE CURRENT PERIOD
00152 !     --------------------------------------------
00153 !
00154       CALL OS( 'X=C     ' , MU2 , SBID , SBID , 0.D0 )
00155       CALL OS( 'X=C     ' , T3  , SBID , SBID , 0.D0 )
00156       CALL OS( 'X=N(Y,Z)' , T1  , PHIR , PHII , CBID )
00157 !
00158 !        ----------------------------------------------------
00159 !        COMPUTES THE WAVE HEIGHT HMU CORRESPONDING TO
00160 !        THE SOLUTION OF THE SYSTEM
00161 !
00162       IF ((.NOT. ALEMON).AND.(.NOT. ALEMUL)) THEN
00163 !       REGULAR WAVE HEIGHT GIVEN BY THE POTENTIAL
00164         IF (COURANT) THEN
00165 !         WE USE WR TO COMPUTE THE FREE SURFACE
00166           CALL OS( 'X=CY    ', X=T2   ,Y=WR, C=2.D0/GRAV)
00167           CALL OS( 'X=YZ    ', X=HHO  ,Y=T1, Z=T2 )
00168         ELSE
00169 !         WE USE OMEGA TO COMPUTE THE FREE SURFACE
00170           CALL OS( 'X=CY    ', HMU , T1   , SBID , 2.D0*OMEGA/GRAV)
00171         ENDIF
00172       ELSE
00173 !       IRREGULAR WAVE HEIGHT : HALE=SIGNIFICANT WAVE HEIGHT
00174         CALL OS( 'X=CY    ', HMU , HALE   , SBID , 1.D0)
00175       ENDIF
00176 !
00177 !     --------------
00178 !     IF BREAKING
00179 !     --------------
00180 !
00181       IF (DEFERL) THEN
00182 !
00183 !     ------------------------------------------------------
00184 !     TESTS IF HMU > HM (THERE IS BREAKING) OR NOT,
00185 !     AND CALCULATES MU2 ACCORDING TO DALLY OR BATTJES & JANSSEN
00186 !     (IF REGULAR WAVES)
00187 !     ------------------------------------------------------
00188 !
00189         IF ((.NOT. ALEMON).AND.(.NOT. ALEMUL)) THEN
00190           DO I = 1,NPOIN
00191             HM = 0.88D0/K%R(I)*TANH(GAMMAS*K%R(I)*H%R(I)/0.88D0)
00192 !
00193 !           HMUE = HMU/SQRT(2)
00194 !           (ENERGETIC WAVE HEIGHT)
00195             HMUE = HMU%R(I)/1.4142D0
00196             HEFF=MIN(HMUE,HM)
00197             HEFF=MAX(HEFF,1.D-5)
00198             Q1 = 1.D-10
00199             Q2 = (HEFF/HM)**2
00200 !           ADDED BY JMH BECAUSE OF THE LOG FUNCTION, LATER ON
00201             Q2 = MAX(Q2,1.D-9)
00202 !
00203 !           ------------
00204 !           COMPUTES QB
00205 !           ------------
00206 !
00207             CALL CALCQB(Q1,Q2,Q3)
00208 !
00209 !           ALGORITHM SPECIFIC TO REGULAR WAVES
00210 !           FOR THE COMPUTATION OF THE RATE OF BREAKING
00211 !
00212             IF (ITERMU.EQ.0) THEN
00213               IF (Q3.LT.0.19D0) THEN
00214                 T3%R(I) = 0.D0
00215               ELSE
00216                 T3%R(I) = 1.D0
00217               ENDIF
00218 !
00219 !             T3 COMPUTED AT ITERMU = 0
00220 !             IS TEMPORARILY STORED IN QB
00221 !
00222                      QB%R(I) = T3%R(I)
00223             ELSE
00224               IF (QB%R(I).EQ.1.D0) THEN
00225                 IF (Q3.LT.0.1D0) THEN
00226                   T3%R(I) = 0.D0
00227                 ELSE
00228                   T3%R(I) = 1.D0
00229                 ENDIF
00230               ENDIF
00231             ENDIF
00232           ENDDO
00233 !
00234 !
00235 !         --------------------------------
00236 !         DALLY AND AL 1985
00237 !         --------------------------------
00238 !
00239           IF (IBREAK.EQ.2) THEN
00240             DO I = 1,NPOIN
00241               HM = 0.88D0/K%R(I)*TANH(GAMMAS*K%R(I)*H%R(I)/0.88D0)
00242               HEFF=MIN(HMU%R(I),HM)
00243               HEFF=MAX(HEFF,1.D-5)
00244               MU2%R(I)=T3%R(I)*KDALLY*
00245      &              (1.D0-(GDALLY*H%R(I)/HEFF)**2)/H%R(I)
00246             ENDDO
00247           ENDIF
00248 !
00249 !     -------------------------------------
00250 !     BATTJES & JANSSEN 1978
00251 !     -------------------------------------
00252 !
00253           IF (IBREAK.EQ.1) THEN
00254             DO I = 1,NPOIN
00255               HM = 0.88D0/K%R(I)*TANH(GAMMAS*K%R(I)*H%R(I)/0.88D0)
00256               HEFF=MIN(HMU%R(I),HM)
00257               HEFF=MAX(HEFF,1.D-5)
00258               MU2%R(I) = T3%R(I)*2.D0*HEFF/(H%R(I)*CG%R(I)*PER)
00259             ENDDO
00260           ENDIF
00261 !
00262 !     -------------------------------------------------------------
00263 !     IRREGULAR WAVES :
00264 !     COMPUTES FIRST QB=T3, PROPORTION OF BREAKING OR BROKEN WAVES,
00265 !     THEN MU2 ACCORDING TO B&J 78 (RANDOM SEAS)
00266 !     -------------------------------------------------------------
00267 !
00268         ELSE
00269           DO I = 1,NPOIN
00270 
00271             HM = 0.88D0/K%R(I)*TANH(GAMMAS*K%R(I)*H%R(I)/0.88D0)
00272 !
00273 !           HMUE = HMU/SQRT (2)
00274 !           (ENERGETIC WAVE HEIGHT)
00275             HMUE = HMU%R(I)/1.4142D0
00276             HEFF=MIN(HMUE,HM)
00277             HEFF=MAX(HEFF,1.D-5)
00278             Q1 = 1.D-10
00279             Q2 = (HEFF/HM)**2
00280 !           ADDED BY JMH BECAUSE OF THE LOG FUNCTION, LATER ON
00281             Q2 = MAX(Q2,1.D-9)
00282 !
00283 !           ------------
00284 !           COMPUTES QB
00285 !           ------------
00286 !
00287             CALL CALCQB(Q1,Q2,Q3)
00288             T3%R(I) = Q3
00289 !
00290 !           -------------------------
00291 !           COMPUTES MU2
00292 !           -------------------------
00293 !
00294             HEFF = MIN((HMU%R(I)/1.4142D0),HM)
00295             HEFF=MAX(HEFF,1.D-5)
00296             MU2%R(I)=ALFABJ*OMEGAM%R(I)*T3%R(I)*((HM/HEFF)**2)/
00297      &             (PI*CG%R(I))
00298           ENDDO
00299 !
00300 !         --------------------------------
00301 !          STOCK QB FOR IRREGULAR WAVES
00302 !         --------------------------------
00303           CALL OS( 'X=Y     ', QB,T3,SBID,CBID)
00304 !
00305         ENDIF
00306 !
00307 !     ------------------
00308 !     END 'IF BREAKING'
00309 !     ------------------
00310 !
00311       ENDIF
00312 !
00313 !     --------------------------------
00314 !     RE-INITIALISES T1 = 0 BECAUSE
00315 !     T1 REPRESENTS MU_FROTTEMENT IN THE FOLLOWING
00316 !     --------------------------------
00317 !
00318       CALL OS( 'X=C     ' , T1 , C , CG , 0.D0 )
00319 !
00320 !     ---------------------
00321 !     IF BOTTOM FRICTION
00322 !     ---------------------
00323 !
00324       IF (FROTTE) THEN
00325 !
00326 !       ------------------------------------------------
00327 !       IF ENTFW=TRUE, THE FRICTION COEFFICIENT FW
00328 !       IS THE SAME EVERYWHERE IN THE DOMAIN
00329 !       ------------------------------------------------
00330 !
00331         IF (ENTFW) THEN
00332           CALL FWSPEC(FW%R,FWCOEF,MESH%X%R,MESH%Y%R,
00333      &                     NPOIN,PRIVE,ZF%R)
00334         ELSE
00335           IF ((.NOT. ALEMON).AND.(.NOT. ALEMUL)) THEN
00336             DO I = 1,NPOIN
00337               CALL CALCFW
00338      &                   (I,H%R,C%R,CG%R,K%R,HMU%R,
00339      &                    NPOIN,OMEGA,GRAV,
00340      &                    VISCO,DIAM90,DIAM50,MVSED,MVEAU,
00341      &                    FORMFR,REGIDO,RICOEF,
00342      &                    ENTREG,ENTRUG,FFW)
00343               FW%R(I) = FFW
00344             ENDDO
00345           ELSE
00346             DO I = 1,NPOIN
00347               CALL CALCFW
00348      &           (I,H%R,C%R,CG%R,K%R,HMU%R,
00349      &            NPOIN,OMEGAM%R(I),GRAV,
00350      &            VISCO,DIAM90,DIAM50,MVSED,MVEAU,
00351      &            FORMFR,REGIDO,RICOEF,
00352      &            ENTREG,ENTRUG,FFW)
00353               FW%R(I) = FFW
00354             ENDDO
00355           ENDIF
00356         ENDIF
00357 !
00358 !       -----------------------------------------
00359 !       COMPUTES THE DISSIPATION COEFFICIENT FOR
00360 !       BOTTOM FRICTION
00361 !       -----------------------------------------
00362 !
00363 !
00364         IF (FORMFR .EQ. 1) THEN
00365 !         ---------------------------------------------------
00366 !         COMPUTES AN EFFECTIVE SPEED
00367 !         UE = 1.2D0*(0.5*((DPHIR/DX)**2 + (DPHIR/DY)**2
00368 !                     +(DPHII/DX)**2 + (DPHII/DY)**2))**0.5
00369 !         UE IS STORED IN T4 HERE
00370 !         ---------------------------------------------------
00371 !
00372           IF ((.NOT. ALEMON).AND.(.NOT. ALEMUL)) THEN
00373             CALL CALCUE
00374 !                   WORK TABLE MODIFIED : T1,T2
00375 !                   RESULTS WORK TABLE  : T4
00376           ELSE
00377             CALL OS( 'X=Y     ', T4 , UEB , SBID  , CBID )
00378           ENDIF
00379 !
00380 !         ----------------------------------------
00381 !         THE DISSIPATION COEFFICIENT MU FOR
00382 !         FRICTION IS STORED IN T1
00383 !         ----------------------------------------
00384 !
00385           CALL OS( 'X=C     ' , T1 , SBID , SBID , 0.D0 )
00386 !
00387           DO I = 1,NPOIN
00388             T1%R(I) = (0.5D0*FW%R(I)*T4%R(I))/
00389      &                (H%R(I)*((COSH(K%R(I)*H%R(I)))**2))
00390             T1%R(I) = T1%R(I)/CG%R(I)
00391           ENDDO
00392         ENDIF
00393 !
00394         IF (FORMFR .EQ. 2) THEN
00395 !         ---------------------------------------------------
00396 !         COMPUTES BOTTOM SPEED FROM STOKES THEORY
00397 !         ---------------------------------------------------
00398           CALL OS( 'X=C     ' , T1 , SBID , SBID , 0.D0 )
00399 !
00400           IF ((.NOT. ALEMON).AND.(.NOT. ALEMUL)) THEN
00401             DO I = 1,NPOIN
00402               T1%R(I) = (2*FW%R(I)*HMU%R(I)*
00403      &                ((OMEGA/SINH(K%R(I)*H%R(I)))**3))
00404               T1%R(I) = T1%R(I)/(3.D0*3.14159D0*GRAV)
00405               T1%R(I) = T1%R(I)/CG%R(I)
00406             ENDDO
00407           ELSE
00408             CALL OS( 'X=Y     ', T4 , UEB , SBID  , CBID )
00409             DO I = 1,NPOIN
00410 !  FOR MEMORY : DIRECT CALCULATION OF MU BASED ON MEAN OMEGA
00411 !  CP suspects that th(kD)=kD approximation was used here
00412 !              T1%R(I) = (2*FW%R(I)*HMU%R(I)*
00413 !     &               ((OMEGAM%R(I)/SINH(K%R(I)*H%R(I)))**3))
00414 !              T1%R(I) = T1%R(I)/(3.D0*3.14159D0*GRAV)
00415 !              T1%R(I) = T1%R(I)/CG%R(I)
00416               T1%R(I) = (0.5D0*FW%R(I)*T4%R(I))/
00417      &                (H%R(I)*((COSH(K%R(I)*H%R(I)))**2))
00418               T1%R(I) = T1%R(I)/CG%R(I)
00419             ENDDO
00420           ENDIF
00421         ENDIF
00422 !
00423 !     -------------------------
00424 !     END 'IF BOTTOM FRICTION'
00425 !     -------------------------
00426 !
00427       ENDIF
00428 !
00429       DO I = 1,NPOIN
00430 !       --------------------------
00431 !       MU = MU_DEFERL + MU_FROTTE
00432 !       --------------------------
00433 !
00434         MU2%R(I) = MU2%R(I) + T1%R(I)
00435       ENDDO
00436 !
00437 !
00438 ! END OF THE CALCULATION OF THE DISSIPATION TERM MU
00439 !
00440 !-----------------------------------------------------------------------
00441 !
00442       RETURN
00443       END SUBROUTINE

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