kepsil.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\kepsil.f
00002 !
00061                      SUBROUTINE KEPSIL
00062 !                    *****************
00063 !
00064      &(AK,EP,AKTILD,EPTILD,AKN,EPN,VISC,CF,U,V,HN,UCONV,VCONV,
00065      & KBOR,EBOR,LIMKEP,IELMK,IELME,
00066      & SMK,SME,TM1,MAK,MAE,CM2,TE1,TE2,NPTFR,DT,
00067      & MESH,T1,T2,T3,TB,CMU,C1,C2,SIGMAK,SIGMAE,ESTAR,SCHMIT,
00068      & KMIN,KMAX,EMIN,EMAX,
00069      & INFOKE,KDIR,MSK,MASKEL,MASKPT,S,SLVK,SLVEP,ICONV,OPTSUP)
00070 !
00071 !***********************************************************************
00072 ! TELEMAC2D   V6P3                                   21/08/2010
00073 !***********************************************************************
00074 !
00075 !
00076 !
00077 !
00078 !
00079 !
00080 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00081 !| AK             |<--| TURBULENT KINETIC ENERGY K AT TIME T(N+1)
00082 !| AKN            |-->| TURBULENT KINETIC ENERGY K AT TIME T(N)
00083 !| AKTILD         |-->| TURBULENT KINETIC ENERGY AFTER ADVECTION
00084 !| C1             |-->| CONSTANT OF K-EPSILON MODEL
00085 !| C2             |-->| CONSTANT OF K-EPSILON MODEL
00086 !| CF             |-->| ADIMENSIONAL FRICTION COEFFICIENT
00087 !| CM2            |<->| MATRIX
00088 !| CMU            |-->| CONSTANT OF K-EPSILON MODEL
00089 !| DT             |-->| TIME STEP
00090 !| EBOR           |<--| TURBULENT ENERGY DISSIPATION AT BOUNDARY
00091 !| EMIN           |-->| MINIMUM EPSILON IF CLIPPING
00092 !| EMAX           |-->| MAXIMUM EPSILON IF CLIPPING
00093 !| EP             |<--| TURBULENT ENERGY DISSIPATION AT TIME T(N+1)
00094 !| EPN            |-->| TURBULENT ENERGY DISSIPATION AT TIME T(N)
00095 !| EPTILD         |-->| TURBULENT ENERGY DISSIPATION AFTER ADVECTION
00096 !| ESTAR          |-->| CONSTANT OF K-EPSILON MODEL
00097 !| HN             |-->| WATER DEPTH AT TIME T(N)
00098 !| ICONV          |-->| TYPE OF ADVECTION ON K AND EPSILON
00099 !|                |   | 1 : CHARACTERISTICS
00100 !|                |   | 2 : SUPG, ...
00101 !| IELME          |-->| TYPE OF ELEMENT FOR K
00102 !| IELMK          |-->| TYPE OF ELEMENT FOR EPSILON
00103 !| INFOKE         |-->| IF YES, INFORMATION ON LINEAR SYSTEMS
00104 !| KBOR           |<--| TURBULENTE KINETIC ENERGY ON BOUNDARIES
00105 !| KDIR           |-->| CONVENTION FOR DIRICHLET POINT
00106 !| KMIN           |-->| MINIMUM K IF CLIPPING
00107 !| KMAX           |-->| K MINIMUM ET MAXIMUM EN CAS DE CLIPPING
00108 !| LIMKEP         |-->| BOUNDARY CONDITIONS ON K AND EPSILON
00109 !| MAE            |<->| MATRIX FOR EPSILON EQUATION
00110 !| MAK            |<->| MATRIX FOR K EQUATION
00111 !| MASKEL         |-->| MASKING OF ELEMENTS
00112 !|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
00113 !| MASKPT         |-->| MASKING PER POINT.
00114 !|                |   | =1. : NORMAL   =0. : MASKED
00115 !| MESH           |-->| MESH STRUCTURE
00116 !| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
00117 !| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
00118 !| OPTSUP         |-->| SUPG OPTION
00119 !| S              |-->| VOID STRUCTURE
00120 !| SCHMIT         |-->| CONSTANT OF K-EPSILON MODEL
00121 !| SIGMAE         |-->| CONSTANT OF K-EPSILON MODEL
00122 !| SIGMAK         |-->| CONSTANT OF K-EPSILON MODEL
00123 !| SLVEP          |-->| STRUCTURE WITH SOLVER OPTIONS FOR E
00124 !| SLVK           |-->| STRUCTURE WITH SOLVER OPTIONS FOR K
00125 !| SME            |<--| RIGHT-HAND SIDE OF EPSILON EQUATION
00126 !| SMK            |<--| RIGHT-HAND SIDE OF K EQUATION
00127 !| T1             |<->| WORK BIEF_OBJ STRUCTURE
00128 !| T2             |<->| WORK BIEF_OBJ STRUCTURE
00129 !| T3             |<->| WORK BIEF_OBJ STRUCTURE
00130 !| TB             |<->| BLOCK OF WORK ARRAYS
00131 !| TE1            |<->| WORK BIEF_OBJ STRUCTURE FOR ELEMENTS
00132 !| TE2            |<->| WORK BIEF_OBJ STRUCTURE FOR ELEMENTS
00133 !| TM1            |<->| DIFFUSION MATRIX
00134 !| UCONV          |-->| X-COMPONENT OF ADVECTION VELOCITY FIELD
00135 !| VCONV          |-->| Y-COMPONENT OF ADVECTION VELOCITY FIELD
00136 !| VISC           |-->| TURBULENT DIFFUSION
00137 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00138 !
00139       USE BIEF
00140 !
00141       IMPLICIT NONE
00142       INTEGER LNG,LU
00143       COMMON/INFO/LNG,LU
00144 !
00145 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00146 !
00147       TYPE(SLVCFG), INTENT(INOUT)  :: SLVK,SLVEP
00148       INTEGER, INTENT(IN)          :: ICONV,NPTFR,KDIR,LIMKEP(NPTFR,2)
00149       INTEGER, INTENT(IN)          :: OPTSUP,IELMK,IELME
00150       LOGICAL, INTENT(IN)          :: INFOKE,MSK
00151       DOUBLE PRECISION, INTENT(IN) :: KMIN,KMAX,EMIN,EMAX,SCHMIT
00152       DOUBLE PRECISION, INTENT(IN) :: CMU,C1,C2,SIGMAK,SIGMAE,ESTAR
00153       DOUBLE PRECISION, INTENT(IN) :: DT
00154 !     MATRIX STRUCTURES
00155       TYPE(BIEF_OBJ), INTENT(INOUT) :: TM1,MAK,MAE,CM2
00156 !     VECTOR STRUCTURES
00157       TYPE(BIEF_OBJ), INTENT(IN)    :: UCONV,VCONV,AKN,EPN,AKTILD,EPTILD
00158       TYPE(BIEF_OBJ), INTENT(IN)    :: HN,VISC,U,V,MASKEL,S,MASKPT,CF
00159       TYPE(BIEF_OBJ), INTENT(IN)    :: KBOR,EBOR
00160       TYPE(BIEF_OBJ), INTENT(INOUT) :: T1,T2,T3,AK,EP,SMK,SME,TE1,TE2
00161 !     MESH STRUCTURE
00162       TYPE(BIEF_MESH) :: MESH
00163 !     BLOCK STRUCTURE
00164       TYPE(BIEF_OBJ) :: TB
00165 !
00166 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00167 !
00168       DOUBLE PRECISION C,SL1,CEPS,USTAR,AGGLOK,AGGLOE,TETAK
00169 !
00170       INTEGER N
00171 !
00172 !-----------------------------------------------------------------------
00173 !
00174       INTRINSIC SQRT,MAX
00175 !
00176 !-----------------------------------------------------------------------
00177 !
00178 !     COMPUTES THE MASS MATRIX AND THE DIFFUSION MATRIX
00179 !
00180       SL1 = 1.D0/DT
00181 !
00182 !     -----------------------------
00183 !     COMPUTES THE MASS MATRIX
00184 !     -----------------------------
00185 !
00186       CALL MATRIX(MAK,'M=N     ','MATMAS          ',IELMK,IELMK,
00187      &            SL1,S,S,S,S,S,S,MESH,MSK,MASKEL)
00188       CALL MATRIX(MAE,'M=N     ','MATMAS          ',IELME,IELME,
00189      &            SL1,S,S,S,S,S,S,MESH,MSK,MASKEL)
00190 !
00191 !     MASS-LUMPING TEST
00192 !
00193       AGGLOK = 1.D0
00194       AGGLOE = 1.D0
00195       IF(AGGLOK.GT.0.001D0) THEN
00196         CALL LUMP(T1,MAK,MESH,AGGLOK)
00197         CALL OM( 'M=CN    ' , MAK , MAK , S  , 1.D0-AGGLOK , MESH )
00198         CALL OM( 'M=M+D   ' , MAK , MAK , T1 , C           , MESH )
00199       ENDIF
00200       IF(AGGLOE.GT.0.001D0) THEN
00201         CALL LUMP(T1,MAE,MESH,AGGLOE)
00202         CALL OM( 'M=CN    ' , MAE , MAE , S  , 1.D0-AGGLOE , MESH )
00203         CALL OM( 'M=M+D   ' , MAE , MAE , T1 , C           , MESH )
00204       ENDIF
00205 !
00206 !     --------------------------------------------------
00207 !     CONCATENATES THE MASS MATRIX: IN T3
00208 !     --------------------------------------------------
00209 !
00210       CALL LUMP(T3,MAK,MESH,DT)
00211 !
00212 !     ---------------------
00213 !     DIFFUSION MATRIX
00214 !     ---------------------
00215 !
00216       CALL MATRIX(TM1,'M=N     ','MATDIF          ',IELMK,IELMK,
00217      &            1.D0,S,S,S,VISC,VISC,VISC,MESH,MSK,MASKEL)
00218 !
00219 !***********************************************************************
00220 !
00221 !     EXPLICIT SOURCE TERMS: T1 FOR K, T2 FOR EPSILON                  *
00222 !                                                                      *
00223 !     EXPLICIT TERM FOR K :                                            *
00224 !                                            3                         *
00225 !                               N           U                          *
00226 !                              K             *
00227 !                              --   +  C  * --  +  PROD
00228 !                              DT       K   H
00229 !
00230 !
00231 !     EXPLICIT TERM FOR EPSILON:
00232 !
00233 !                                            4
00234 !                                N          U              N
00235 !                              EP            *           EP
00236 !                              --   +  C  * --  +  C   * -- * PROD
00237 !                              DT       E    2      E1    N
00238 !                                           H            K
00239 !
00240 !
00241 !                     2        2           2
00242 !                  DU       DV     DU   DV           N
00243 !      PROD = ( 2*(--) + 2*(--) + (-- + --)  ) * VISC
00244 !                  DX       DY     DY   DX
00245 !
00246 !
00247 !                           N
00248 !                         EP
00249 !      THE TERM  +  C1  * -- * PROD   IS WRITTEN AS :
00250 !                          N
00251 !                         K
00252 !
00253 !                               N
00254 !                   C1 * CMU * K * PROD / VISC
00255 !
00256 !***********************************************************************
00257 !
00258 !     --------------------------------
00259 !     TAKES ADVECTION INTO ACCOUNT
00260 !     --------------------------------
00261 !
00262       IF(ICONV.EQ.1) THEN
00263 !
00264         CALL MATVEC('X=AY    ',SMK,MAK,AKTILD,C,MESH)
00265         CALL MATVEC('X=AY    ',SME,MAE,EPTILD,C,MESH)
00266 !
00267       ELSEIF(ICONV.EQ.2) THEN
00268 !
00269         CALL MATVEC('X=AY    ',SMK,MAK,AKN,C,MESH)
00270         CALL MATVEC('X=AY    ',SME,MAE,EPN,C,MESH)
00271 !       CENTERED SEMI-IMPLICIT ADVECTION TERM : MATRIX
00272         CALL MATRIX(CM2,'M=N     ','MATVGR          ',IELMK,IELMK,
00273      &              1.D0,S,S,S,UCONV,VCONV,VCONV,MESH,MSK,MASKEL)
00274 !       SUPG CONTRIBUTION
00275         IF(OPTSUP.EQ.1) THEN
00276 !         CLASSICAL SUPG
00277           CALL KSUPG(TE1,TE2,1.D0,UCONV,VCONV,MESH)
00278           CALL MATRIX(CM2,'M=M+N   ','MASUPG          ',IELMK,IELMK,
00279      &                1.D0,TE1,TE2,S,UCONV,VCONV,VCONV,
00280      &                MESH,MSK,MASKEL)
00281         ELSEIF(OPTSUP.EQ.2) THEN
00282 !         MODIFIED SUPG
00283           CALL MATRIX(CM2,'M=M+N   ','MAUGUG          ',IELMK,IELMK,
00284      &                0.5D0*DT,S,S,S,UCONV,VCONV,VCONV,
00285      &                MESH,MSK,MASKEL)
00286         ENDIF
00287 !       END OF SUPG CONTRIBUTION
00288 !       EXPLICIT RIGHT HAND SIDES
00289         TETAK=0.6
00290         CALL MATVEC( 'X=X+CAY ',SMK,CM2,AKN,TETAK-1.D0,MESH)
00291         CALL MATVEC( 'X=X+CAY ',SME,CM2,EPN,TETAK-1.D0,MESH)
00292 !       ADDS SUPG MATRIX TO MAK AND MAE
00293         CALL OM( 'M=X(M)  ' , MAK , MAK , S , C , MESH )
00294         CALL OM( 'M=M+CN  ' , MAK , CM2 , S , TETAK , MESH )
00295         CALL OM( 'M=X(M)  ' , MAE , MAE , S , C , MESH )
00296         CALL OM( 'M=M+CN  ' , MAE , CM2 , S , TETAK , MESH )
00297 !
00298       ELSE
00299 !
00300         IF(LNG.EQ.1) WRITE(LU,100) ICONV
00301         IF(LNG.EQ.2) WRITE(LU,101) ICONV
00302 100     FORMAT(1X,'KEPSIL : FORME DE LA CONVECTION INCONNUE : ',1I4)
00303 101     FORMAT(1X,'KEPSIL: UNKNOWN TYPE OF ADVECTION:',1I4)
00304         CALL PLANTE(1)
00305         STOP
00306 !
00307       ENDIF
00308 !
00309 !     ------------------------------------------------
00310 !     CREATION TERM - BOTTOM FRICTION (P)
00311 !     ------------------------------------------------
00312 !
00313 !     BEWARE : MISSES 1.D0/(0.5D0*CF)**0.75 TO GET TRUE CEPS
00314 !     (TAKEN INTO ACCOUNT AFTERWARD)
00315 !
00316       CEPS = C2 * SQRT(CMU) / SQRT( ESTAR*SCHMIT )
00317 !
00318       CALL CPSTVC(SMK,T1)
00319       CALL CPSTVC(SMK,T2)
00320       DO N=1,T1%DIM1
00321         USTAR = SQRT(0.5D0*CF%R(N)*(U%R(N)**2+V%R(N)**2))
00322 !       T1 : PKV TERM OF THE EQUATION FOR K
00323         T1%R(N)=USTAR**3/MAX(SQRT(0.5D0*CF%R(N))*HN%R(N),1.D-6)
00324 !       LIMITS THE GROWTH OF K DUE TO BOTTOM FRICTION
00325 !       TO 50% PER TIMESTEP
00326         T1%R(N) = MIN (T1%R(N) , EP%R(N) + 0.5D0*AK%R(N)/DT )
00327 !       T2 : PEV TERM OF THE EQUATION FOR EPSILON
00328         T2%R(N) = CEPS * USTAR**4 /
00329      &          MAX(((0.5D0*CF%R(N))**0.75D0)*HN%R(N)**2,1.D-6)
00330 !       LIMITS THE GROWTH OF E TO 50% PER TIMESTEP
00331         T2%R(N) = MIN(T2%R(N),C2*EP%R(N)**2/MAX(AK%R(N),KMIN)
00332      &          +0.5D0*EP%R(N)/DT )
00333       ENDDO
00334 !
00335       CALL OS( 'X=XY    ' , T1  , T3 , T3 , C )
00336       CALL OS( 'X=X+Y   ' , SMK , T1 , T1 , C )
00337       CALL OS( 'X=XY    ' , T2  , T3 , T3 , C )
00338       CALL OS( 'X=X+Y   ' , SME , T2 , T2 , C )
00339 !
00340 !     -----------------------------------
00341 !     CREATION TERM - SHEAR
00342 !     -----------------------------------
00343 !
00344       CALL VECTOR(T1,'=','PRODF           ',IELMK,
00345      &            1.D0,S,S,S,U,V,S,MESH,MSK,MASKEL)
00346       CALL OS( 'X=XY    ' , T1  , VISC , VISC , C )
00347 !
00348 ! TEST JMH : LIMITS TURBULENCE FOR SHALLOW DEPTHS ( < 2CM )
00349 !
00350       DO N=1,SMK%DIM1
00351         IF(HN%R(N).LT.0.02D0) T1%R(N)=0.D0
00352       ENDDO
00353 !
00354 ! END OF TEST
00355 !
00356       CALL OS( 'X=X+Y   ' , SMK , T1   , T1   , C )
00357 !
00358       CALL VECTOR(T2,'=','PRODF           ',IELMK,
00359      &            CMU*C1,S,S,S,U,V,S,MESH,MSK,MASKEL)
00360       CALL OS( 'X=XY    ' , T2  , AK , AK , C )
00361       CALL OS( 'X=X+Y   ' , SME , T2 , T2 , C )
00362 !
00363 ! TEST JMH : LIMITS TURBULENCE FOR SHALLOW DEPTHS ( < 2CM )
00364 !
00365       DO N=1,SMK%DIM1
00366         SMK%R(N) = SMK%R(N) * (MIN(HN%R(N),0.02D0)/0.02D0)**2
00367       ENDDO
00368 !
00369 ! END OF TEST
00370 !
00371 !***********************************************************************
00372 !     IMPLICIT SOURCE TERMS : T1 FOR K , T2 FOR EPSILON                *
00373 !                                                                      *
00374 !     IMPLICIT TERM FOR K :           +      EP(N)/K(N) * K (N+1)      *
00375 !     IMPLICIT TERM FOR EPSILON:      + C2 * EP(N)/K(N) * EP(N+1)      *
00376 !***********************************************************************
00377 !
00378       CALL OS( 'X=Y/Z   ',T1,EP,AK,C  ,IOPT=2,INFINI=0.D0,ZERO=KMIN)
00379       CALL OS( 'X=CY    ',T2,T1,T1,C2 )
00380 !
00381 !     ---------------------------------------------------
00382 !     INTEGRATES THESE SOURCE TERMS IN THE MATRICES
00383 !     ---------------------------------------------------
00384 !
00385       CALL OS( 'X=XY    ' , T1 , T3 , T3 , C )
00386       CALL OS( 'X=XY    ' , T2 , T3 , T3 , C )
00387 !
00388 !     -------------------------------------------
00389 !     ADDS TO THE DIAGONAL OF THE MASS MATRIX
00390 !     -------------------------------------------
00391 !
00392       CALL OM( 'M=M+D   ' , MAK , MAK , T1 , C , MESH )
00393       CALL OM( 'M=M+D   ' , MAE , MAE , T2 , C , MESH )
00394 !
00395 !***********************************************************************
00396 !
00397 !     COMBINES THE MASS AND DIFFUSION MATRICES
00398 !
00399 !     MAK = MAK + TM1/SIGMAK
00400 !     MAE = MAE + TM1/SIGMAE
00401 !
00402 !***********************************************************************
00403 !
00404       CALL OM( 'M=M+CN  ' , MAK , TM1 , S , 1.D0/SIGMAK , MESH )
00405       CALL OM( 'M=M+CN  ' , MAE , TM1 , S , 1.D0/SIGMAE , MESH )
00406 !
00407 !***********************************************************************
00408 !     DIRICHLET TYPE BOUNDARY CONDITIONS
00409 !***********************************************************************
00410 !
00411       IF(NPTFR.GT.0) THEN
00412         CALL DIRICH(AK,MAK,SMK,KBOR,LIMKEP(1,1),TB,MESH,KDIR,MSK,MASKPT)
00413         CALL DIRICH(EP,MAE,SME,EBOR,LIMKEP(1,2),TB,MESH,KDIR,MSK,MASKPT)
00414       ENDIF
00415 !
00416 !***********************************************************************
00417 !     SOLVES THE TWO OBTAINED SYSTEMS
00418 !***********************************************************************
00419 !
00420       CALL SOLVE(AK,MAK,SMK,TB,SLVK ,INFOKE,MESH,TM1)
00421       CALL SOLVE(EP,MAE,SME,TB,SLVEP,INFOKE,MESH,TM1)
00422 !
00423 !***********************************************************************
00424 !     CLIPS SMALL VALUES                                               *
00425 !***********************************************************************
00426 !
00427       CALL CLIP(AK,0.D0,.TRUE.,KMAX,.FALSE.,0)
00428       CALL CLIP(EP,EMIN,.TRUE.,EMAX,.FALSE.,0)
00429 !
00430 !***********************************************************************
00431 !
00432       RETURN
00433       END

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