diff3d.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac3d\diff3d.f
00002 !
00123                      SUBROUTINE DIFF3D
00124 !                    *****************
00125 !
00126      &(FD,FC,FN,VISCF,SIGMAF,S0F,YAS0F,S1F,YAS1F,
00127      & FBORL,FBORF,FBORS,AFBORL,AFBORF,AFBORS,
00128      & BFBORL,BFBORF,BFBORS,LIFBOF,LIFBOL,LIFBOS,
00129      & FMIN,CLIMIN,FMAX,CLIMAX,SCHCF,SCHDF,SLVDIF,TRBAF,INFO,
00130      & NEWDIF,DT,T2_01,T2_02,T2_03,T3_01,T3_02,T3_03,T3_04,
00131      & NPOIN2,NPOIN3,INCHYD,SEM3D,YASEM3D,IT1,NPTFR3,NBOR3,MASKPT,
00132      & TRAV3,MESH2D,MESH3D,MTRA1,MTRA2,IELM3,MSUPG,IELM2H,IELM2V,
00133      & MDIFF,MATR2H,MASKBR,SVIDE,MSK,MASKEL,H,
00134      & NPLAN,OPTBAN,OPTDIF,TETADI,YAWCC,WCC,AGGLOD,VOLU,
00135      & YASCE,NSCE,FSCE,SOURCES,TETASUPG,VELOCITY,RAIN,PLUIE,TRAIN,
00136      & SIGMAG,IPBOT,SETDEP)
00137 !
00138 !***********************************************************************
00139 ! TELEMAC3D   V7P0                                  21/08/2010
00140 !***********************************************************************
00141 !
00142 !
00143 !
00144 !
00145 !
00146 !
00147 !
00148 !
00149 !
00150 !
00151 !
00152 !
00153 !
00154 !
00155 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00156 !| AFBORF         |-->| LOGARITHMIC LAW FOR COMPONENT ON THE BOTTOM:
00157 !|                |   |  NU*DF/DN = AFBORF*U + BFBORF
00158 !| AFBORL         |-->| LOGARITHMIC LAW FOR COMPONENT ON THE
00159 !|                |   | LATERAL BOUNDARIES:
00160 !|                |   | NU*DF/DN = AFBORL*U + BFBORL
00161 !| AFBORS         |-->| LOGARITHMIC LAW FOR COMPONENT AT THE SURFACE:
00162 !|                |   | NU*DF/DN = AFBORS*U + BFBORS
00163 !| AGGLOD         |-->| MASS-LUMPING IN DIFFUSION
00164 !| BFBORF         |-->| LOGARITHMIC LAW FOR COMPONENT ON THE BOTTOM:
00165 !|                |   |  NU*DF/DN = AFBORF*U + BFBORF
00166 !| BFBORL         |-->| LOGARITHMIC LAW FOR COMPONENT ON THE
00167 !|                |   | LATERAL BOUNDARIES:
00168 !|                |   | NU*DF/DN = AFBORL*U + BFBORL
00169 !| BFBORS         |-->| LOGARITHMIC LAW FOR COMPONENT AT THE SURFACE:
00170 !|                |   | NU*DF/DN = AFBORS*U + BFBORS
00171 !| CLIMAX         |-->| LOGICAL FOR CLIPPING (MAX VALUE)
00172 !| DT             |-->| TIME STEP
00173 !| FBORF          |<->| DIRICHLET CONDITIONS ON F AT THE BOTTOM
00174 !| FBORL          |<->| DIRICHLET CONDITIONS ON F ON LATERAL BOUNDARIES
00175 !| FBORS          |<->| DIRICHLET CONDITIONS ON F AT THE SURFACE
00176 !| FC             |<->| VARIABLE AFTER CONVECTION
00177 !| FD             |<->| VARIABLE AFTER DIFFUSION
00178 !| FMAX           |-->| MAX CLIPPING VALUE
00179 !| FMIN           |-->| MIN CLIPPING VALUE
00180 !| FN             |<->| VARIABLE F AT TIME N
00181 !| FSCE           |-->| SOURCE TERM OF F
00182 !| H              |-->| WATER DEPTH
00183 !| IELM2H         |-->| DISCRETISATION TYPE FOR 2D HORIZONTAL MESH
00184 !| IELM2V         |-->| DISCRETISATION TYPE FOR 2D VERTICAL MESH
00185 !| IELM3          |-->| DISCRETISATION TYPE FOR 3D
00186 !| INCHYD         |-->| IF YES, HYDROSTATIC INCONSISTENCY FILTER
00187 !| INFO           |-->| INFORMATIONS FOR SOLVERS
00188 !| IPBOT          |-->| PLANE NUMBER OF LAST CRUSHED PLANE (0 IF NONE)
00189 !| IT1            |<->| BIEF_OBJ STRUCTURES FOR INTEGER ARRAYS
00190 !| LIFBOF         |<->| TYPE OF BOUNDARY CONDITIONS AT THE BOTTOM
00191 !| LIFBOL         |<->| TYPE OF BOUNDARY CONDITIONS ON LATERAL BOUNDARIES
00192 !| LIFBOS         |<->| TYPE OF BOUNDARY CONDITIONS AT THE SURFACE
00193 !| MASKBR         |-->| 3D MASK ON LATERAL BOUNDARIES
00194 !| MASKEL         |-->| MASKING OF ELEMENTS
00195 !|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
00196 !| MASKPT         |-->| MASKING PER POINT.
00197 !|                |   | =1. : NORMAL   =0. : MASKED
00198 !| MATR2H         |<->| WORK MATRIX 2DH
00199 !| MDIFF          |<->| DIFFUSION MATRIX
00200 !| MESH2D         |<->| 2D MESH
00201 !| MESH3D         |<->| 3D MESH
00202 !| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
00203 !| MSUPG          |<->| NON SYMMETRIC SUPG MATRIX
00204 !| MTRA1          |<->| 3D WORK MATRIX
00205 !| MTRA2          |<->| 3D WORK MATRIX
00206 !| NBOR3          |-->| GLOBAL NUMBER OF 3D BOUNDARY POINTS
00207 !| NEWDIF         |-->| RECALCULATE OR NOT DIFFUSION MATRIX
00208 !| NPLAN          |-->| NUMBER OF PLANES IN THE 3D MESH
00209 !| NPOIN2         |-->| NUMBER OF 2D POINTS
00210 !| NPOIN3         |-->| NUMBER OF 3D POINTS
00211 !| NPTFR3         |-->| NUMBER OF LATERAL BOUNDARY POINTS IN 3D
00212 !| NSCE           |-->| NUMBER OF GIVEN POINTS FOR SOURCES
00213 !| OPTBAN         |-->| OPTION FOR TIDAL FLATS, IF 1, FREE SURFACE
00214 !|                |   | MODIFIED AND PIECE-WISE LINEAR
00215 !| OPTDIF         |-->| OPTION FOR THE DIFFUSION OF F
00216 !| PLUIE          |-->| RAIN IN M/S MULTIPLIED BY VOLU2D
00217 !| RAIN           |-->| IF YES, THERE IS RAIN OR EVAPORATION
00218 !| S0F            |<->| EXPLICIT SOURCE TERM (DIM=F/T)
00219 !| S1F            |<->| IMPLICIT SOURCE TERM (DIM=1/T)
00220 !| SCHCF          |-->| ADVECTION SCHEME OF F
00221 !| SCHDF          |-->| DIFFUSION SCHEME OF F
00222 !| SEM3D          |<->| SECOND MEMBERS (RIGHT HAND SIDE)
00223 !|                |   | FOR THE LINEAR EQUATIONS 3D
00224 !| SETDEP         |-->| ADVECTION-DIFFUSION SCHEME WITH SETTLING VELOCITY
00225 !|                |   | 0 : IMPLICIT TREATED IN DIFFUSION
00226 !|                |   | 1 : EXPLICIT TREATED IN ADVECTION
00227 !| SIGMAF         |-->| COEFFICIENT OF VISCOSITY REDUCTION
00228 !|                |   | ONLY USED FOR K AND EPSILON
00229 !| SIGMAG         |-->| LOGICAL FOR GENERALISED SIGMA TRANSFORMATION
00230 !| SLVDIF         |-->| SOLVER FOR DIFFUSION OF VELOCITIES
00231 !| SOURCES        |-->| RIGHT HAND SIDE OF CONTINUITY EQUATION WHEN SOURCES
00232 !| SVIDE          |-->| VOID STRUCTURE
00233 !| T2_01          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00234 !| T2_02          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00235 !| T2_03          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00236 !| T3_02          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00237 !| T3_03          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00238 !| T3_04          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00239 !| TETADI         |<->| IMPLICITATION COEFFICIENT OF THE DIAGONAL
00240 !|                |   | OF DIFFUSION IF OPTDIF = 2
00241 !|                |   | IMPLICITATION COEFFICIENT OF DIFFUSION
00242 !|                |   | IF OPTDIF = 1
00243 !| TETASUPG       |-->| IMPLICITATION COEFFICIENT FOR SUPG
00244 !| TRAIN          |-->| VALUE OF TRACER IN RAIN
00245 !| TRAV3          |<->| 3D WORK ARRAYS
00246 !| TRBAF          |-->| TREATMENT ON TIDAL FLATS FOR F
00247 !| VELOCITY       |-->| IF TRUE, COMPONENT IS VELOCITY: NOT USED
00248 !| VISCF          |<->| VISCOSITY COEFFICIENTS
00249 !|                |   | VISCF(*,1 OU 2) HORIZONTAL VISCOSITY
00250 !|                |   | VISCF(*,3)      VERTICAL VISCOSITY
00251 !| VOLU           |-->| VOLUME AROUND POINTS AT TIME N+1
00252 !| WCC            |-->| VELOCITY (POSITIVE IN 6.3 IF SEDIMENT SETTLING VELOCITY)
00253 !| YAS0F          |-->| LOGICAL TO TAKE INTO ACCOUNT S0F TERM IN DIFF3D
00254 !| YAS1F          |-->| LOGICAL TO TAKE INTO ACCOUNT S1F TERM IN DIFF3D
00255 !| YASCE          |-->| IF TRUE, THERE IS SOURCE
00256 !| YASEM3D        |-->| IF TRUE, RIGHT HAND SIDE HAS BEEN PARTLY
00257 !|                |   | COMPUTED BEFORE CALLING DIFF3D
00258 !| YAWCC          |-->| LOGICAL TO TAKE INTO ACCOUNT WCC FOR SEDIMENT
00259 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00260 !
00261       USE BIEF
00262       USE DECLARATIONS_TELEMAC
00263 !
00264 !     PROVISOIRE
00265 !
00266       USE DECLARATIONS_TELEMAC3D, ONLY:UCONV,VCONV,WSCONV,DM1,ZCONV
00267 !
00268       IMPLICIT NONE
00269       INTEGER LNG,LU
00270       COMMON/INFO/LNG,LU
00271 !
00272 !-----------------------------------------------------------------------
00273 !
00274       INTEGER, INTENT(IN)             :: SCHCF,SCHDF,TRBAF,OPTDIF,NSCE
00275       INTEGER, INTENT(IN)             :: NPOIN2,SETDEP
00276       TYPE(BIEF_OBJ), INTENT(INOUT)   :: FD,FC,FN,S0F,S1F,VISCF
00277       TYPE(BIEF_OBJ), INTENT(IN)      :: LIFBOL,LIFBOF,LIFBOS
00278       TYPE(BIEF_OBJ), INTENT(IN)      :: FBORL,FBORF,FBORS
00279       TYPE(BIEF_OBJ), INTENT(IN)      :: AFBORL,AFBORF,AFBORS
00280       TYPE(BIEF_OBJ), INTENT(IN)      :: BFBORL,BFBORF,BFBORS,H,WCC
00281       DOUBLE PRECISION, INTENT(IN)    :: SIGMAF,FMIN,FMAX,DT,AGGLOD
00282       DOUBLE PRECISION, INTENT(IN)    :: FSCE(NSCE),TETASUPG
00283       DOUBLE PRECISION, INTENT(INOUT) :: TETADI
00284       LOGICAL, INTENT(IN)             :: CLIMIN,CLIMAX,YASCE,YAS0F
00285       LOGICAL, INTENT(IN)             :: INFO,NEWDIF,YASEM3D,YAS1F
00286       LOGICAL, INTENT(IN)             :: SIGMAG
00287       TYPE(SLVCFG)                    :: SLVDIF
00288 !
00289       LOGICAL, INTENT(IN)             :: INCHYD,MSK,YAWCC,VELOCITY,RAIN
00290       TYPE(BIEF_OBJ), INTENT(IN)      :: MASKPT,MASKBR,MASKEL,VOLU
00291       TYPE(BIEF_OBJ), INTENT(IN)      :: NBOR3,SVIDE,SOURCES
00292       TYPE(BIEF_OBJ), INTENT(INOUT)   :: T3_01,T3_02,T3_03,T3_04
00293       TYPE(BIEF_OBJ), INTENT(INOUT)   :: T2_01,T2_02,T2_03
00294       TYPE(BIEF_MESH), INTENT(INOUT)  :: MESH3D,MESH2D
00295       INTEGER, INTENT(IN)             :: NPOIN3,IELM3,IELM2H,IELM2V
00296       INTEGER, INTENT(IN)             :: NPTFR3,NPLAN,OPTBAN
00297       INTEGER, INTENT(IN)             :: IPBOT(NPOIN2)
00298       TYPE(BIEF_OBJ), INTENT(INOUT)   :: SEM3D,IT1,TRAV3,MTRA1,MTRA2
00299       TYPE(BIEF_OBJ), INTENT(INOUT)   :: MSUPG,MDIFF
00300       TYPE(BIEF_OBJ), INTENT(INOUT)   :: MATR2H
00301       DOUBLE PRECISION, INTENT(IN)    :: PLUIE(NPOIN2),TRAIN
00302 !
00303 !-----------------------------------------------------------------------
00304 !
00305       INTEGER IPOIN3,IPOIN2,IPTFR3,ITERD,NITERD,IS,IIS,I,I2D,NP
00306       DOUBLE PRECISION C,AGGLO,QTOT,VTOT,FTOT
00307       CHARACTER(LEN=16) FORMUL
00308       DOUBLE PRECISION, POINTER :: VOLUME(:)
00309 !
00310 !-----------------------------------------------------------------------
00311 !
00312       NITERD = SLVDIF%NITMAX
00313       IF(OPTDIF.EQ.2) TETADI = 0.D0
00314 !
00315 !=======================================================================
00316 !   MASS MATRIX
00317 !=======================================================================
00318 !
00319 !     COMPUTES MTRA2 : MASS MATRIX DIVIDED BY DT
00320 !
00321       FORMUL='MATMAS          '
00322       CALL MATRIX(MTRA2, 'M=N     ', FORMUL, IELM3, IELM3, 1.D0/DT,
00323      &         SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,SVIDE, MESH3D, MSK, MASKEL)
00324 !
00325 ! MASS-LUMPING OF THE MASS MATRIX (FOR THE MOMENT ONLY IF
00326 ! EXPLICIT DIFFUSION OR 2 X IMPLICIT DIAGONAL)
00327 ! OF IF COMPATIBILITY WITH PSI SCHEME NEEDED (USE OF VOLU WHICH
00328 ! IS DIFFERENT FROM INTEGRAL OF TEST FUNCTIONS IF MASS-LUMPING
00329 ! IS DONE IN TELEMAC-2D)
00330 !
00331       AGGLO=AGGLOD
00332       IF(TETADI.LT.0.001D0.OR.SCHCF.EQ.ADV_NSC.OR.SCHCF.EQ.ADV_PSI
00333      &                    .OR.SCHCF.EQ.ADV_LPO.OR.SCHCF.EQ.ADV_LPO_TF
00334      &                    .OR.SCHCF.EQ.ADV_NSC_TF) AGGLO=1.D0
00335       IF(AGGLO.GT.0.001D0) THEN
00336 !       VOLU/DT REPLACES MTRA2 LUMPED
00337         CALL OS('X=CY    ',X=T3_01,Y=VOLU,C=AGGLO/DT)
00338 !       IF(MSK) THERE IS A RISK THAT VOLU=0 SOMEWHERE
00339 !               THIS MAY CAUSE PROBLEMS IN LINEAR SYSTEMS
00340         IF(MSK) CALL OS('X=+(Y,C)',X=T3_01,Y=T3_01,C=1.D-3)
00341 !       CALL LUMP(T3_01,MTRA2,MESH3D,AGGLO)
00342         CALL OM( 'M=CN    ',MTRA2,MTRA2 , SVIDE , 1.D0-AGGLO , MESH3D )
00343         CALL OM( 'M=M+D   ',MTRA2,MTRA2 , T3_01 , 0.D0       , MESH3D )
00344       ENDIF
00345 !
00346 !=======================================================================
00347 !   CASES OF ADVECTION OF SETTLING VELOCITY WITH EXPLICIT SCHEME
00348 !   FC IS AT THE BEGINNING THE RESULT OF ADVECTION WITHOUT SETTLING
00349 !   VELOCITY, THEN THE EFFECT OF SETTLING VELOCITY IS ADDED
00350 !=======================================================================
00351 !
00352       IF(SETDEP.EQ.2) THEN
00353         CALL OS('X=Y     ',X=T3_02,Y=FC)
00354         CALL SED_FALL(FC,T3_02,WCC,MESH2D,MESH3D,DT,VOLU,
00355      &                IPBOT,NPOIN2,NPOIN3,NPLAN,T3_01)
00356       ENDIF
00357 !
00358 !=======================================================================
00359 !   SECOND MEMBER OF ADVECTION SCHEME + EXPLICIT SOURCE TERM
00360 !=======================================================================
00361 !
00362 !   COMPUTES SEM3D = DT * S0F + FC
00363 !
00364       IF(S0F%TYPR.NE.'0'.AND.YAS0F) THEN
00365         CALL OS( 'X=Y+CZ  ' , T3_04 , FC , S0F , DT )
00366       ELSE
00367         CALL OS( 'X=Y     ' , X=T3_04 , Y=FC )
00368       ENDIF
00369 !
00370       IF(YASEM3D) THEN
00371 !       HAS STARTED COMPUTING SEM3D IN TRISOU OR SOURCE
00372         CALL MATVEC ('X=X+CAY  ',SEM3D, MTRA2,T3_04,DT, MESH3D)
00373 !                                                   DT BECAUSE MTRA2
00374 !                                                   HAS A FACTOR 1/DT
00375       ELSE
00376         CALL MATVEC ('X=AY    ',SEM3D, MTRA2,T3_04, C, MESH3D)
00377       ENDIF
00378 !
00379 !=======================================================================
00380 !
00381 !     SUPG SCHEME: TEST OF UPWINDING OF THE TIME DERIVATIVE
00382 !
00383       IF(SCHCF.EQ.ADV_SUP) THEN
00384 !       HERE OPTSUP=2 ONLY
00385         CALL MATRIX(MTRA1,'M=TN    ','MATVGR 2        ',
00386      &              IELM3,IELM3,0.5D0,DM1,ZCONV,SVIDE,
00387      &              UCONV,VCONV,WSCONV,MESH3D,MSK,MASKEL)
00388 !       COEFFICIENT 1/DT HAS ALREADY BEEN PUT IN MTRA1 (0.5=0.5*DT/DT)
00389         CALL OM('M=X(M)  ',MTRA2,MTRA2,SVIDE,C,MESH3D)
00390         CALL OM('M=M+N   ',MTRA2,MTRA1,SVIDE,C,MESH3D)
00391 !       RIGHT HAND SIDE PART
00392 !                                            =FN IN THIS OPTION
00393         CALL MATVEC ('X=X+CAY  ',SEM3D,MTRA1,FC,1.D0,MESH3D)
00394       ENDIF
00395 !
00396 !=======================================================================
00397 !
00398 !     SOURCES INSIDE THE DOMAIN
00399 !
00400       IF(YASCE.AND.NSCE.GT.0) THEN
00401         DO IS=1,NSCE
00402 !         IF INTAKE FSCE=F, SO NO EXTRA TERM
00403           IIS=IS
00404 !         IN PARALLEL MODE SOURCES WITHOUT PARCOM
00405           IF(NCSIZE.GT.1) IIS=IIS+NSCE
00406           DO I=1,NPOIN3
00407 !           EXPLICIT SOURCE TERM
00408             SEM3D%R(I) = SEM3D%R(I)
00409      &                 + MAX(SOURCES%ADR(IIS)%P%R(I),0.D0)*
00410      &            (FSCE(IS)-(1.D0-TETASUPG)*FN%R(I))
00411 !           IMPLICIT SOURCE TERM : SEE BELOW
00412           ENDDO
00413         ENDDO
00414       ENDIF
00415 !
00416 !=======================================================================
00417 !
00418 !     RAIN (ALL TRACERS) - EXPLICIT PART
00419 !     VALUE OF TRACER IN RAIN TAKEN INTO ACCOUNT ONLY IF RAIN POSITIVE
00420 !     NOT IN CASE OF EVAPORATION, HENCE THE MAX(PLUIE,0)
00421 !     THERE IS NO COEFFICIENT (1-TETASUPG) ON THIS TERM BECAUSE IT HAS
00422 !     NO IMPLICIT PART
00423 !
00424       IF(RAIN) THEN
00425         DO I=1,NPOIN2
00426           IPOIN3=NPOIN3-NPOIN2+I
00427           SEM3D%R(IPOIN3)=SEM3D%R(IPOIN3)
00428      &                    -PLUIE(I)*(1.D0-TETASUPG)*FN%R(IPOIN3)
00429      &                    +MAX(PLUIE(I),0.D0)*TRAIN
00430         ENDDO
00431       ENDIF
00432 !
00433 !=======================================================================
00434 !
00435 !     EXPLICIT ADVECTION TERM:
00436 !
00437       IF(SCHCF.EQ.ADV_SUP) THEN
00438         CALL MATVEC('X=X+CAY ',SEM3D,MSUPG,FN,TETASUPG-1.D0,MESH3D)
00439       ENDIF
00440 !
00441 !=======================================================================
00442 !   IMPLICIT SOURCE TERM (MASS-LUMPED AND ADDED TO THE DIAGONAL)
00443 !=======================================================================
00444 !
00445 ! JMH PROPOSITION (LUMPED IMPLICIT SOURCE TERM)
00446 !                  BEWARE THE + SIGN: THE TERM S1F HAS A + SIGN WHEN
00447 !                  TO THE LEFT OF THE = SIGN
00448 !                 (CAREFUL: OPPOSITE IN CVDFTR)
00449 !                  VOLU USED FOR MASS-CONSERVATION
00450 !                 (COMPATIBILITY WITH CONTINUITY EQUATION)
00451       IF(S1F%TYPR.NE.'0'.AND.YAS1F) THEN
00452         CALL OS( 'X=YZ    ' , X=T3_01 , Y=S1F , Z=VOLU )
00453         CALL OM( 'M=M+D   ' , MTRA2 , MTRA2 , T3_01 , C , MESH3D )
00454       ENDIF
00455 !
00456 !=======================================================================
00457 !
00458 !     SOURCES INSIDE THE DOMAIN
00459 !
00460       IF(YASCE.AND.NSCE.GT.0) THEN
00461         DO IS=1,NSCE
00462 !         IF INTAKE FSCE=T, SO NO EXTRA TERM
00463           IIS=IS
00464 !         IN PARALLEL MODE SOURCES WITHOUT PARCOM
00465           IF(NCSIZE.GT.1) IIS=IIS+NSCE
00466           DO I=1,NPOIN3
00467 !           IMPLICIT SOURCE TERM
00468             MTRA2%D%R(I)=MTRA2%D%R(I)+
00469      &                   MAX(SOURCES%ADR(IIS)%P%R(I),0.D0)*TETASUPG
00470           ENDDO
00471         ENDDO
00472       ENDIF
00473 !
00474 !     RAIN (ALL TRACERS) - IMPLICIT PART
00475 !
00476       IF(RAIN) THEN
00477         DO I=1,NPOIN2
00478           IPOIN3=NPOIN3-NPOIN2+I
00479           MTRA2%D%R(IPOIN3)=MTRA2%D%R(IPOIN3)+PLUIE(I)*TETASUPG
00480         ENDDO
00481       ENDIF
00482 !
00483 !=======================================================================
00484 !   DIFFUSION MATRIX + BOUNDARY TERMS
00485 !=======================================================================
00486 !
00487       IF(SCHDF.NE.0) THEN
00488 !
00489         IF(INFO) THEN
00490           IF(SCHCF.EQ.ADV_SUP) THEN
00491             IF(LNG.EQ.1) WRITE(LU,*) 'DIFFUSION DE ',FN%NAME,
00492      &                               ' AVEC CONVECTION PAR SUPG'
00493             IF(LNG.EQ.2) WRITE(LU,*) 'DIFFUSION OF ',FN%NAME,
00494      &                               ' WITH SUPG ADVECTION'
00495           ELSE
00496             IF(LNG.EQ.1) WRITE(LU,*) 'DIFFUSION DE ',FN%NAME
00497             IF(LNG.EQ.2) WRITE(LU,*) 'DIFFUSION OF ',FN%NAME
00498           ENDIF
00499         ENDIF
00500 !
00501         IF(NEWDIF) THEN
00502 !         RELEASE 5.7
00503 !         FORMUL='MATDIF          '
00504 !         RELEASE 5.8 : MONOTONICITY ENSURED
00505           FORMUL='MATDIF       MON'
00506           IF(INCHYD) FORMUL(7:7)='2'
00507           CALL MATRIX(MDIFF,'M=N     ',FORMUL,IELM3,IELM3,1.D0,
00508      &                VISCF%ADR(1)%P,VISCF%ADR(2)%P,VISCF%ADR(3)%P,
00509      &                SVIDE,SVIDE,SVIDE,MESH3D,MSK,MASKEL)
00510         ENDIF
00511 !
00512 !=======================================================================
00513 !
00514         IF(OPTDIF.EQ.1) THEN
00515 !
00516 !         SEMI-IMPLICITATION OF THE DIFFUSION
00517 !
00518           IF(TETADI.LT.0.999D0) THEN
00519 !           IF(TETADI.GT.0.001) THEN
00520             CALL MATVEC ('X=X+CAY  ',SEM3D,MDIFF,FN,
00521      &                   (TETADI-1.D0)/SIGMAF,MESH3D)
00522           ENDIF
00523 !         IF TETADI=0, NO NEED TO ADD MDIFF
00524           IF(TETADI.GT.0.001D0) THEN
00525             CALL OM('M=M+CN  ',MTRA2,MDIFF,SVIDE,TETADI/SIGMAF,MESH3D)
00526           ENDIF
00527 !
00528         ENDIF
00529 !
00530 ! TAKES THE IMPLICIT BOUNDARY TERMS INTO ACCOUNT (FRICTION FOR EXAMPLE)
00531 ! ---------------------------------------------
00532 !
00533 !   LATERAL FACES : (MASKBR SET ACCORDING TO MASKEL IN MASK3D)
00534 !  (MASS-LUMPED FORM)
00535 !
00536         IF(AFBORL%TYPR.NE.'0') THEN
00537           CALL VECTOR(T3_02,'=','MASBAS          ',IELM2V,-1.D0,SVIDE,
00538      &                SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH3D,MSK,MASKEL)
00539           CALL OS('X=XY    ',T3_02,AFBORL,AFBORL,0.D0)
00540           CALL OSDB( 'X=X+Y   ' , MTRA2%D , T3_02 , T3_02, C , MESH3D )
00541         ENDIF
00542 !
00543 !  (NON MASS-LUMPED FORM, BUT MATR2V NOW SUPPRESSED)
00544 !        FORMUL='FMATMA          '
00545 !        CALL MATRIX
00546 !    *   (MATR2V, 'M=N     ', FORMUL, IELM2V, IELM2V, -1.D0, AFBORL,
00547 !    *    SVIDE, SVIDE, SVIDE, SVIDE, SVIDE, MESH3D, MSK, MASKBR)
00548 !        CALL OM('M=M+N   ',MTRA2,MATR2V,SVIDE,C,MESH3D)
00549 !
00550 !   BOTTOM (MASS-LUMPED FORM AS IN 2D):
00551 !
00552         IF(AFBORF%TYPR.NE.'0') THEN
00553           CALL VECTOR(T2_03, '=','MASBAS          ',IELM2H,-1.D0,SVIDE,
00554      &                SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH2D,MSK,MASKEL)
00555           CALL OV('X=XY    ',T2_03%R,AFBORF%R,AFBORF%R,0.D0,NPOIN2)
00556 !         DRY ZONES OR CRUSHED ELEMENTS
00557 !         SEE EQUIVALENT TREATMENT IN WAVE_EQUATION
00558 !         JMH CORRECTION TIDAL FLATS 05/11/2013
00559           IF(SIGMAG.OR.OPTBAN.EQ.1) THEN
00560             DO IPOIN2 = 1,NPOIN2
00561 !             FLUXES ON THE BOTTOM MUST NOT BE TAKEN INTO ACCOUNT ON TIDAL FLATS
00562 !             E.G. NO SETTLING VELOCITY IF NO WATER!
00563 !             IF(IPBOT(IPOIN2).NE.NPLAN-1) THEN
00564 !               IF NOT A TIDAL FLAT... WE CAN HOWEVER HAVE CRUSHED POINTS
00565 !               SAME TREATMENT FROM PLANE 1 UP TO ACTUAL BOTTOM PLANE
00566 !               THE CRUSHED POINTS ARE TREATED AS THE POINT ON ACTUAL BOTTOM PLANE ABOVE THEM
00567 !               WARNING 1 (JMH ON 04/06/2014)
00568 !               THIS MAY NOT CORRECT FOR SEDIMENT, AS DEPOSITION WOULD BE TAKEN INTO ACCOUNT
00569 !               SEVERAL TIMES. TO BE INVESTIGATED ON MORE TESTS.
00570 !               THIS IS RATHER MEANT FOR VELOCITIES, TO STOP THEM ALL.
00571                 DO NP=0,IPBOT(IPOIN2)
00572                   I=NP*NPOIN2+IPOIN2
00573                   MTRA2%D%R(I)=MTRA2%D%R(I)+T2_03%R(IPOIN2)
00574                 ENDDO
00575 !             ENDIF
00576 !             PROPOSED ALTERNATIVE IMPLEMENTATION, ONLY FIRST FREE POINT TREATED
00577 !             AND NOTHING DONE ON TIDAL FLATS.
00578 !             IF(IPBOT(IPOIN2).NE.NPLAN-1) THEN
00579 !               I=IPBOT(IPOIN2)*NPOIN2+IPOIN2
00580 !               MTRA2%D%R(I)=MTRA2%D%R(I)+T2_03%R(IPOIN2)
00581 !             ENDIF
00582             ENDDO
00583           ELSE
00584             DO IPOIN2 = 1,NPOIN2
00585               MTRA2%D%R(IPOIN2)=MTRA2%D%R(IPOIN2)+T2_03%R(IPOIN2)
00586             ENDDO
00587           ENDIF
00588         ENDIF
00589 !
00590 !       SURFACE (MASS-LUMPED FORM):
00591 !
00592         IF(AFBORS%TYPR.NE.'0') THEN
00593           CALL VECTOR(T2_03, '=','MASBAS          ',IELM2H,-1.D0,SVIDE,
00594      &                SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH2D,MSK,MASKEL)
00595           CALL OV('X=XY    ',T2_03%R,AFBORS%R,AFBORS%R,0.D0,NPOIN2)
00596           DO I=1,NPOIN2
00597             MTRA2%D%R(NPOIN3-NPOIN2+I)=
00598      &      MTRA2%D%R(NPOIN3-NPOIN2+I)+T2_03%R(I)
00599           ENDDO
00600 !         CALL OV('X=X+Y   ',MTRA2%D%R(NPOIN3-NPOIN2+1:NPOIN3),
00601 !    &                       T2_03%R,T2_03%R,0.D0,NPOIN2)
00602         ENDIF
00603 !
00604 !       TAKES THE EXPLICIT BOUNDARY TERMS INTO ACCOUNT
00605 !       ---------------------------------------------
00606 !
00607         CALL T3D_STRESS(SEM3D,'X=X+Y   ',T2_01,T3_02,
00608      &                  BFBORL,BFBORF,BFBORS,NPOIN2,NPOIN3,MESH2D,
00609      &                  MESH3D,IELM3,IELM2H,IELM2V,SVIDE,
00610      &                  MSK,MASKBR,MASKEL,IPBOT,SIGMAG,OPTBAN,NPLAN)
00611 !
00612 !=======================================================================
00613 !   SEDIMENT-SPECIFIC ++++ This is for WC > 0 downwards
00614 !                                D
00615 !   THE MATRIX + PSI1(J) * WCC * --( PSI2(I) ) IS ADDED
00616 !                                DZ
00617 !                                                     D          N+1
00618 !   IT IS AN INTEGRATION BY PART OF TERM :  PSI2(I) * --( WCC * C    )
00619 !                                                     DZ
00620 !
00621 !   CORRESPONDING BOUNDARY TERMS ARE THE DEPOSITION
00622 !   THIS TERM IS INCLUDED IN ATABOF !!!
00623 !
00624 !   NOTE: IT IS DONE IF AND ONLY IF SED. DIFFUSION IS REQUIRED !
00625 !=======================================================================
00626 !
00627         IF(YAWCC.AND.SETDEP.EQ.0) THEN
00628 !         FOR BOUNDARY TERMS, SEE SUBROUTINE FLUSED
00629 !         NOTE: MATWC IS NOT PROGRAMMED WITH TETRAHEDRA...
00630           CALL MATRIX(MTRA1,'M=N     ','MATWC           ',
00631      &                IELM3,IELM3,1.D0,WCC,
00632      &                SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH3D,MSK,MASKEL)
00633           CALL OM('M=X(M)  ',MTRA2,MTRA2,SVIDE,C,MESH3D)
00634           CALL OM('M=M+N   ',MTRA2,MTRA1,SVIDE,C,MESH3D)
00635         ENDIF
00636 !
00637 !=======================================================================
00638 !   END OF DIFFUSION MATRIX + BOUNDARY TERMS
00639 !=======================================================================
00640 !
00641       ENDIF
00642 !
00643 !=======================================================================
00644 !
00645 !   ADDS SUPG MATRIX IF SCHCF=ADV_SUP
00646 !
00647 !=======================================================================
00648 !
00649       IF(SCHCF.EQ.ADV_SUP.AND.OPTDIF.EQ.1) THEN
00650         IF(MTRA2%TYPEXT.EQ.'S') THEN
00651           CALL OM('M=X(M)  ',MTRA2,MTRA2,SVIDE,C,MESH3D)
00652         ENDIF
00653         CALL OM ('M=M+CN  ',MTRA2, MSUPG, SVIDE, TETASUPG, MESH3D )
00654       ENDIF
00655 !
00656 !=======================================================================
00657 !
00658 !   BOUNDARY CONDITIONS + PRECONDITIONING + SOLVER
00659 !
00660 !=======================================================================
00661 !
00662 !   BOUNDARY CONDITIONS FOR BOUNDARY POINTS (POINTS OF TYPE DIRICHLET)
00663 !
00664       CALL CPSTVC(MTRA2%D,T3_03)
00665 !
00666       DO IPOIN3 = 1,NPOIN3
00667         IT1%I(IPOIN3) = KDDL
00668         T3_03%R(IPOIN3) = 0.D0
00669       ENDDO
00670 !
00671 !   LATERAL BOUNDARY CONDITIONS
00672 !
00673       DO IPTFR3 = 1,NPTFR3
00674         IF(LIFBOL%I(IPTFR3).EQ.KENT.OR.
00675      &    LIFBOL%I(IPTFR3).EQ.KENTU.OR.LIFBOL%I(IPTFR3).EQ.KADH) THEN
00676           IT1%I(NBOR3%I(IPTFR3)) = KDIR
00677           T3_03%R(NBOR3%I(IPTFR3)) = FBORL%R(IPTFR3)
00678         ENDIF
00679       ENDDO
00680 !
00681 !   BOTTOM AND SURFACE BOUNDARY CONDITIONS
00682 !
00683       DO IPOIN2 = 1,NPOIN2
00684         IF(LIFBOF%I(IPOIN2).EQ.KENT.OR.LIFBOF%I(IPOIN2).EQ.KADH) THEN
00685           IT1%I(IPOIN2) = KDIR
00686           T3_03%R(IPOIN2) = FBORF%R(IPOIN2)
00687         ENDIF
00688         IF(LIFBOS%I(IPOIN2).EQ.KENT.OR.LIFBOS%I(IPOIN2).EQ.KADH) THEN
00689           IT1%I(NPOIN3-NPOIN2+IPOIN2) = KDIR
00690           T3_03%R(NPOIN3-NPOIN2+IPOIN2) = FBORS%R(IPOIN2)
00691         ENDIF
00692       ENDDO
00693 !
00694 !   CRUSHED POINTS AND TIDAL FLATS: FIRST TREATED AS DIRICHLET
00695 !
00696       IF(SIGMAG.OR.OPTBAN.EQ.1) THEN
00697         IF(NCSIZE.GT.1) THEN
00698 !         ONLY DIFFERENCE : VALUE OF MTRA2 EQUAL TO FAC INSTEAD OF 1
00699           CALL OS('X=Y     ',X=T3_02,Y=VOLU)
00700 !         T3_02 WILL BE THE ASSEMBLED VOLUME
00701           CALL PARCOM(T3_02,2,MESH3D)
00702           DO I2D=1,NPOIN2
00703             IF(IPBOT(I2D).GT.0) THEN
00704               IF(LIFBOF%I(I2D).EQ.KENT.OR.LIFBOF%I(I2D).EQ.KADH) THEN
00705 !               DIRICHLET POINT ON THE BOTTOM: ALL POINTS UP TO THE
00706 !               FIRST FREE ARE TREATED AS DIRICHLET WITH FBORF VALUE
00707                 DO I=0,IPBOT(I2D)
00708                   IPOIN3 = I2D+I*NPOIN2
00709                   MTRA2%D%R(IPOIN3)=MESH3D%FAC%R(IPOIN3)
00710                   IT1%I(IPOIN3)   = KDIR
00711                   T3_03%R(IPOIN3) = FBORF%R(I2D)
00712                 ENDDO
00713               ELSE
00714 !               POINTS WITH NO VOLUME PROVISIONALLY SET TO FN
00715 !               TESTED UP TO FIRST FREE POINT (IPBOT+1),
00716 !               WHICH SHOULD NOT HAVE VOLUME=0 EXCEPT ON TIDAL
00717 !               FLATS BECAUSE IN THIS CASE IPBOT=NPLAN-1
00718 !               DO I=0,IPBOT(I2D)-1
00719                 DO I=0,IPBOT(I2D)
00720                   IPOIN3 = I2D+I*NPOIN2
00721                   IF(T3_02%R(IPOIN3).LT.1.D-10) THEN
00722                     MTRA2%D%R(IPOIN3)=MESH3D%FAC%R(IPOIN3)
00723                     IT1%I(IPOIN3)   = KDIR
00724                     T3_03%R(IPOIN3) = FN%R(IPOIN3)
00725                   ENDIF
00726                 ENDDO
00727               ENDIF
00728             ENDIF
00729           ENDDO
00730         ELSE
00731           DO I2D=1,NPOIN2
00732             IF(IPBOT(I2D).GT.0) THEN
00733               IF(LIFBOF%I(I2D).EQ.KENT.OR.LIFBOF%I(I2D).EQ.KADH) THEN
00734 !               DIRICHLET POINT ON THE BOTTOM: ALL POINTS UP TO THE
00735 !               FIRST FREE ARE TREATED AS DIRICHLET WITH FBORF VALUE
00736                 DO I=0,IPBOT(I2D)
00737                   IPOIN3 = I2D+I*NPOIN2
00738                   MTRA2%D%R(IPOIN3)=1.D0
00739                   IT1%I(IPOIN3) = KDIR
00740                   T3_03%R(IPOIN3) = FBORF%R(I2D)
00741                 ENDDO
00742               ELSE
00743 !               POINTS WITH NO VOLUME PROVISIONALLY SET TO FN
00744 !               TESTED UP TO FIRST FREE POINT (IPBOT+1),
00745 !               WHICH SHOULD NOT HAVE VOLUME=0 EXCEPT ON TIDAL
00746 !               FLATS BECAUSE IN THIS CASE IPBOT=NPLAN-1
00747 !               DO I=0,IPBOT(I2D)-1
00748                 DO I=0,IPBOT(I2D)
00749                   IPOIN3 = I2D+I*NPOIN2
00750                   IF(VOLU%R(IPOIN3).LT.1.D-10) THEN
00751                     MTRA2%D%R(IPOIN3)=1.D0
00752                     IT1%I(IPOIN3) = KDIR
00753                     T3_03%R(IPOIN3) = FN%R(IPOIN3)
00754                   ENDIF
00755                 ENDDO
00756               ENDIF
00757             ENDIF
00758           ENDDO
00759         ENDIF
00760       ENDIF
00761 !
00762 !-----------------------------------------------------------------------
00763 !
00764 !---> IMPLICIT DIFFUSION (STANDARD)
00765 !
00766       IF(OPTDIF.EQ.1) THEN
00767 !
00768         CALL DIRI01(FD,MTRA2,SEM3D,T3_03,IT1%I,T3_01,T3_02,MESH3D,
00769      &              KDIR,MSK,MASKPT)
00770 !
00771 !   SOLVES THE LINEAR SYSTEM:
00772 !
00773         IF(TETADI.GT.0.001D0.OR.SCHCF.EQ.ADV_SUP) THEN
00774           CALL SOLVE (FD, MTRA2, SEM3D, TRAV3, SLVDIF,
00775      &                INFO, MESH3D, MTRA1)
00776         ELSE
00777           IF(NCSIZE.GT.1) CALL PARCOM(MTRA2%D,2,MESH3D)
00778           CALL OS('X=Y/Z   ',FD,SEM3D,MTRA2%D,0.D0,2,0.D0,1.D-10)
00779         ENDIF
00780 !
00781 !---> EXPLICIT DIFFUSION / PARTIAL IMPLICITATION OF 2XDIAGONAL
00782 !
00783       ELSEIF(OPTDIF.EQ.2) THEN
00784 !---> BACKS UP IF NITERD> 1
00785         CALL OS( 'X=Y     ' , T3_04 , SEM3D , SVIDE , C )
00786 !---> NITER=1
00787         CALL OS('X=CY    ',T3_01,MDIFF%D,T3_01,2.D0/SIGMAF )
00788         CALL OM( 'M=M+D   ',MTRA2,MTRA2,T3_01,C,MESH3D)
00789 ! CALL OM( 'M=N     ',MTRA1,MTRA2,MTRA2,C,MESH3D)
00790 !
00791         CALL MATVEC
00792      &       ('X=X+CAY ',SEM3D,MDIFF,FN,-1.D0/SIGMAF,MESH3D)
00793         CALL OS('X=X+CYZ ',SEM3D,FN,MDIFF%D,2.D0/SIGMAF)
00794 !
00795         CALL DIRI01(FD,MTRA2,SEM3D,T3_03,IT1%I,T3_01,T3_02,MESH3D,
00796      &            KDIR,MSK,MASKPT)
00797 !
00798 !-->  SOLVES THE LINEAR SYSTEM:
00799 !
00800         CALL LUMP(T3_01, MTRA2, MESH3D, 1.D0)
00801         CALL OS( 'X=Y/Z   ' , FD , SEM3D , T3_01 , C )
00802 !
00803 !---> GAUSS-SEIDEL ITERATIONS
00804 !
00805         IF(NITERD.GE.2) THEN
00806         DO ITERD = 2,NITERD
00807           CALL OS( 'X=Y     ' , SEM3D , T3_04 , SVIDE , C )
00808           CALL MATVEC
00809      &       ('X=X+CAY ',SEM3D,MDIFF,FD,-1.D0/SIGMAF,MESH3D)
00810           CALL OS('X=X+CYZ ',SEM3D,FD,MDIFF%D,2.D0/SIGMAF)
00811 !
00812           CALL DIRI01(FD,MTRA2,SEM3D,T3_03,IT1%I,T3_01,T3_02,MESH3D,
00813      &              KDIR,MSK,MASKPT)
00814 !
00815 !---> SOLVES THE DIAGONAL LINEAR SYSTEM:
00816 !
00817           CALL LUMP(T3_01, MTRA2, MESH3D, 1.D0)
00818           CALL OS( 'X=Y/Z   ' , FD , SEM3D , T3_01 , C )
00819 !
00820         ENDDO
00821         ENDIF
00822 !
00823       ENDIF
00824 !
00825 !=======================================================================
00826 !
00827 !   TREATS THE POINTS WHICH ARE COMPLETELY DRY
00828 !    (BY DEFAULT FORCED TO 0)
00829 !
00830       IF(MSK.AND.TRBAF.EQ.1) THEN
00831         DO IPOIN3 = 1,NPOIN3
00832           IF(MASKPT%R(IPOIN3).LT.0.5D0) FD%R(IPOIN3) = FN%R(IPOIN3)
00833         ENDDO
00834       ENDIF
00835 !
00836 !=======================================================================
00837 !
00838 !   CRUSHED POINTS AND THEIR FREE POINT ABOVE MUST HAVE THE SAME VALUE
00839 !   IT IS DONE HERE WITH MASS CONSERVATION: THE VOLUMIC AVERAGE IS
00840 !   TAKEN. TO AVOID SPURIOUS GRADIENTS...
00841 !   HEAVY BUT SAVES TIME IN COMPLEX APPLICATIONS!
00842 !
00843 !=======================================================================
00844 !
00845 !     WARNING 2 (JMH 04/06/2014)
00846 !     THIS SECTION MAY CHANGE THE MASS-BALANCE, AS THE POINT TAKEN INTO
00847 !     ACCOUNT FOR ENTERING AFBORF WILL HAVE ITS VALUE CHANGED.
00848 !     TO SOLVE THE PROBLEM, EITHER WE REMOVE THIS OR WE HAVE A MORE
00849 !     SOPHISTICATED COMPUTATION OF THE MASS-BALANCE, WITH TERMS DEPENDING
00850 !     ON AFBORF TREATED BEFORE THIS SECTION.
00851 !
00852       IF(SIGMAG.OR.OPTBAN.EQ.1) THEN
00853         IF(NCSIZE.GT.1) THEN
00854           CALL OS('X=Y     ',X=T3_02,Y=VOLU)
00855 !         T3_02 IS THE ASSEMBLED VOLUME
00856           CALL PARCOM(T3_02,2,MESH3D)
00857           VOLUME => T3_02%R
00858         ELSE
00859           VOLUME => VOLU%R
00860         ENDIF
00861         DO I2D=1,NPOIN2
00862           IF(IPBOT(I2D).GT.0) THEN
00863             VTOT=0.D0
00864             QTOT=0.D0
00865             DO I=0,IPBOT(I2D)
00866               VTOT=VTOT+VOLUME(I2D+I*NPOIN2)
00867               QTOT=QTOT+VOLUME(I2D+I*NPOIN2)*FD%R(I2D+I*NPOIN2)
00868             ENDDO
00869             IF(VTOT.GT.1.D-10) THEN
00870               FTOT=QTOT/VTOT
00871             ELSE
00872               FTOT=FD%R(I2D+IPBOT(I2D)*NPOIN2)
00873             ENDIF
00874             DO I=0,IPBOT(I2D)
00875               FD%R(I2D+I*NPOIN2)=FTOT
00876             ENDDO
00877           ENDIF
00878         ENDDO
00879       ENDIF
00880 !
00881 !=======================================================================
00882 !
00883 !   CLIPS F (IF REQUESTED WITH CLIMIN AND CLIMAX)
00884 !
00885 !=======================================================================
00886 !
00887       CALL CLIP(FD,FMIN,CLIMIN,FMAX,CLIMAX,0)
00888 !
00889 !=======================================================================
00890 !
00891       RETURN
00892       END

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