cvdf3d.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac3d\cvdf3d.f
00002 !
00106                      SUBROUTINE CVDF3D
00107 !                    *****************
00108 !
00109      &(FD,FC,FN,VISCF,SIGMAF,S0F,YAS0F,S1F,YAS1F,
00110      & FBORL,FBORF,FBORS,AFBORL,AFBORF,AFBORS,
00111      & BFBORL,BFBORF,BFBORS,LIFBOL,LIFBOF,LIFBOS,
00112      & FLUXF,FLUEXT,FLUEXTPAR,FMIN,CLIMIN,FMAX,CLIMAX,
00113      & SCHCF,SCHDF,SLVDIF,TRBAF,INFOR,NEWDIF,CALFLU,
00114      & T2_01,T2_02,T2_03,
00115      & T3_01,T3_02,T3_03,T3_04,MESH3D,IKLE3,MASKEL,MTRA1,
00116      & W1,NPTFR3,MMURD,MURD_TF,VOLU,VOLUPAR,VOLUN,VOLUNPAR,
00117      & NBOR3,NPOIN3,NPOIN2,DT,MSK,NELEM2,NELEM3,
00118      & NPLAN,LV,IELM3,MSUPG,IELM2H,IELM2V,MDIFF,MTRA2,
00119      & INCHYD,MASKBR,MASKPT,SEM3D,YASEM3D,SVIDE,IT1,IT2,
00120      & TRAV3,MESH2D,MATR2H,H,OPTBAN,OPTDIF,TETADI,
00121      & YAWCC,WCC,AGGLOD,NSCE,SOURCES,FSCE,NUMLIQ,DIRFLU,NFRLIQ,
00122      & VOLUT,ZT,ZPROP,RAIN,PLUIE,PARAPLUIE,TRAIN,
00123      & FLODEL,FLOPAR,SIGMAG,IPBOT,MAXADV,FLUDPT,FLUDP,FLUER,
00124      & VOLU2D,V2DPAR,SETDEP)
00125 !
00126 !***********************************************************************
00127 ! TELEMAC3D   V7P0                                   21/08/2010
00128 !***********************************************************************
00129 !
00130 !
00131 !
00132 !
00133 !
00134 !
00135 !
00136 !
00137 !
00138 !
00139 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00140 !| AFBORF         |-->| LOGARITHMIC LAW FOR COMPONENT ON THE BOTTOM:
00141 !|                |   |  NU*DF/DN = AFBORF*U + BFBORF
00142 !| AFBORL         |-->| LOGARITHMIC LAW FOR COMPONENT ON THE
00143 !|                |   | LATERAL BOUNDARIES:
00144 !|                |   | NU*DF/DN = AFBORL*U + BFBORL
00145 !| AFBORS         |-->| LOGARITHMIC LAW FOR COMPONENT AT THE SURFACE:
00146 !|                |   | NU*DF/DN = AFBORS*U + BFBORS
00147 !| AGGLOD         |-->| MASS-LUMPING IN DIFFUSION
00148 !| BFBORF         |-->| LOGARITHMIC LAW FOR COMPONENT ON THE BOTTOM:
00149 !|                |   |  NU*DF/DN = AFBORF*U + BFBORF
00150 !| BFBORL         |-->| LOGARITHMIC LAW FOR COMPONENT ON THE
00151 !|                |   | LATERAL BOUNDARIES:
00152 !|                |   | NU*DF/DN = AFBORL*U + BFBORL
00153 !| BFBORS         |-->| LOGARITHMIC LAW FOR COMPONENT AT THE SURFACE:
00154 !|                |   | NU*DF/DN = AFBORS*U + BFBORS
00155 !| CALFLU         |-->| INDICATE IF FLUX IS CALCULATED FOR BALANCE
00156 !| CLIMAX         |-->| LOGICAL FOR CLIPPING (MAX VALUE)
00157 !| CLIMIN         |-->| LOGICAL FOR CLIPPING (MIN VALUE)
00158 !| DIRFLU         |-->| TREATMENT OF FLUXES AT THE BOUNDARIES
00159 !| DT             |-->| TIME STEP
00160 !|    SETDEP      |-->| INTEGER FOR SETTLING AND DEPOSITION CONVECTION SCHEME
00161 !| FBORF          |<->| DIRICHLET CONDITIONS ON F AT THE BOTTOM
00162 !| FBORL          |<->| DIRICHLET CONDITIONS ON F ON LATERAL BOUNDARIES
00163 !| FBORS          |<->| DIRICHLET CONDITIONS ON F AT THE SURFACE
00164 !| FC             |<->| VARIABLE AFTER CONVECTION
00165 !| FD             |<->| VARIABLE AFTER DIFFUSION
00166 !| FLODEL         |-->| FLUXES BY SEGMENT
00167 !| FLOPAR         |-->| FLUXES BY SEGMENT, ASSEMBLED IN PARALLEL
00168 !| FLUDPT         |-->| IMPLICIT DEPOSITION FLUX (SEDIMENT)
00169 !| FLUER          |-->| EROSION FLUX (SEDIMENT)
00170 !| FLUEXT         |-->| OUTPUT FLUX BY NODE
00171 !| FLUEXTPAR      |-->| OUTPUT FLUX BY NODE IN PARALLEL
00172 !| FLUXF          |<->| FLUX FOR F
00173 !| FMAX           |-->| MAX CLIPPING VALUE
00174 !| FMIN           |-->| MIN CLIPPING VALUE
00175 !| FN             |<->| VARIABLE F AT TIME N
00176 !| FSCE           |-->| SOURCE TERM OF F
00177 !| H              |-->| WATER DEPTH AT TIME
00178 !| IELM2H         |-->| DISCRETISATION TYPE FOR 2D HORIZONTAL MESH
00179 !| IELM2V         |-->| DISCRETISATION TYPE FOR 2D VERTICAL MESH
00180 !| IELM3          |-->| DISCRETISATION TYPE FOR 3D
00181 !| IKLE3          |-->| GLOBAL 3D CONNECTIVITY
00182 !| INCHYD         |-->| IF YES, HYDROSTATIC INCONSISTENCY FILTER
00183 !| INFOR          |-->| INFORMATIONS FOR SOLVERS
00184 !| IPBOT          |-->| PLANE NUMBER OF LAST CRUSHED PLANE (0 IF NONE)
00185 !| IT1            |<->| BIEF_OBJ STRUCTURES FOR INTEGER ARRAYS
00186 !| IT2            |<->| BIEF_OBJ STRUCTURES FOR INTEGER ARRAYS
00187 !| LIFBOF         |<->| TYPE OF BOUNDARY CONDITIONS AT THE BOTTOM
00188 !| LIFBOL         |<->| TYPE OF BOUNDARY CONDITIONS ON LATERAL BOUNDARIES
00189 !| LIFBOS         |<->| TYPE OF BOUNDARY CONDITIONS AT THE SURFACE
00190 !| LV             |-->| VECTOR LENGTH OF THE MACHINE
00191 !| MASKBR         |-->| 3D MASK ON LATERAL BOUNDARIES
00192 !| MASKEL         |-->| MASKING OF ELEMENTS
00193 !|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
00194 !| MASKPT         |-->| MASKING PER POINT.
00195 !|                |   | =1. : NORMAL   =0. : MASKED
00196 !| MATR2H         |<->| WORK MATRIX 2DH
00197 !| MAXADV         |-->| MAXIMUM NUMBER OF ITERATIONS FOR ADVECTION SCHEMES
00198 !| MDIFF          |<->| DIFFUSION MATRIX
00199 !| MESH2D         |<->| 2D MESH
00200 !| MESH3D         |<->| 3D MESH
00201 !| MMURD          |<->| NON SYMMETRIC MURD MATRIX
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 !| MURD_TF        |<->| MURD MATRIX FOR TIDAL FLAT
00207 !| NBOR3          |-->| GLOBAL NUMBER OF 3D BOUNDARY POINTS
00208 !| NELEM2         |-->| NUMBER OF ELEMENTS IN 2D
00209 !| NELEM3         |-->| NUMBER OF ELEMENTS IN 3D
00210 !| NEWDIF         |-->| RECALCULATE OR NOT DIFFUSION MATRIX
00211 !| NFRLIQ         |-->| NUMBER OF LIQUID BOUNDARIES
00212 !| NPLAN          |-->| NUMBER OF PLANES IN THE 3D MESH
00213 !| NPOIN2         |-->| NUMBER OF 2D POINTS
00214 !| NPOIN3         |-->| NUMBER OF 3D POINTS
00215 !| NPTFR3         |-->| NUMBER OF LATERAL BOUNDARY POINTS IN 3D
00216 !| NSCE           |-->| NUMBER OF GIVEN POINTS FOR SOURCES
00217 !| NUMLIQ         |-->| LIQUID BOUNDARY NUMBER OF BOUNDARY POINTS
00218 !| OPTBAN         |-->| OPTION FOR TIDAL FLATS, IF 1, FREE SURFACE
00219 !|                |   | MODIFIED AND PIECE-WISE LINEAR
00220 !| OPTDIF         |-->| OPTION FOR THE DIFFUSION
00221 !| PARAPLUIE      |-->| RAIN IN M/S MULTIPLIED BY VOLU2D
00222 !|                |   | (IN ASSEMBLED MODE IN PARALLEL)
00223 !| PLUIE          |-->| RAIN IN M/S MULTIPLIED BY VOLU2D
00224 !| RAIN           |-->| IF YES, THERE IS RAIN OR EVAPORATION
00225 !| S0F            |<->| EXPLICIT SOURCE TERM (DIM=F/T)
00226 !| S1F            |<->| IMPLICIT SOURCE TERM (DIM=1/T)
00227 !| SCHCF          |-->| ADVECTION SCHEME OF F
00228 !| SCHDF          |-->| DIFFUSION SCHEME OF F
00229 !| SEM3D          |<->| SECOND MEMBERS (RIGHT HAND SIDE)
00230 !|                |   | FOR THE LINEAR EQUATIONS 3D
00231 !| SIGMAF         |-->| COEFFICIENT OF VISCOSITY REDUCTION
00232 !| SIGMAG         |-->| LOGICAL FOR GENERALISED SIGMA TRANSFORMATION
00233 !| SLVDIF         |-->| SOLVER FOR DIFFUSION OF VELOCITIES
00234 !| SOURCES        |-->| RIGHT HAND SIDE OF CONTINUITY EQUATION WHEN SOURCES
00235 !| SVIDE          |-->| VOID STRUCTURE
00236 !| T2_01          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00237 !| T2_02          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00238 !| T2_03          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00239 !| T3_01          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00240 !| T3_02          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00241 !| T3_03          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00242 !| T3_04          |<->| BIEF_OBJ STRUCTURE FOR LOCAL WORK
00243 !| TETADI         |<->| IMPLICITATION RATE FOR DIFFUSION
00244 !| TRAIN          |-->| VALUE OF TRACER IN RAIN
00245 !| TRAV3          |<->| BLOCK OF 3D BIEF_OBJ STRUCTURES (AT LEAST 10)
00246 !| TRBAF          |-->| TREATMENT ON TIDAL FLATS FOR F
00247 !| VISCF          |<->| VISCOSITY COEFFICIENTS
00248 !|                |   | VISCF(*,1 OR 2) HORIZONTAL VISCOSITY
00249 !|                |   | VISCF(*,3)      VERTICAL VISCOSITY
00250 !| VOLU           |-->| VOLUME AROUND POINTS AT TIME N+1
00251 !| VOLUN          |-->| VOLUME AROUND POINTS AT TIME N
00252 !| VOLUNPAR       |-->| VOLUME AROUND POINTS AT TIME N, IN PARALLEL
00253 !| VOLUPAR        |-->| VOLUME AROUND POINTS AT TIME N+1, IN PARALLEL
00254 !| VOLUT          |<->| VOLUME AFTER SEMI-IMPLICITATION FOR TRACER
00255 !| W1             |<->| WORK ARRAY (MATRICES COMPUTATION...)
00256 !| WCC            |-->| VELOCITY (NEGATIVE IF SEDIMENT SETTLING VELOCITY)
00257 !| YASEM3D        |-->| IF TRUE, RIGHT HAND SIDE HAS BEEN PARTLY
00258 !|                |   | COMPUTED BEFORE CALLING DIFF3D
00259 !| YAS0F          |-->| LOGICAL TO TAKE INTO ACCOUNT S0F TERM IN DIFF3D
00260 !| YAS1F          |-->| LOGICAL TO TAKE INTO ACCOUNT S1F TERM IN DIFF3D
00261 !| YAWCC          |-->| LOGICAL TO TAKE INTO ACCOUNT WCC FOR SEDIMENT
00262 !| ZPROP          |-->| VERTICAL COORDINATES FOR PROPAGATION STEP
00263 !| ZT             |<->| Z: DISTRIBUTION
00264 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00265 !
00266       USE BIEF
00267       USE DECLARATIONS_TELEMAC
00268       USE INTERFACE_TELEMAC3D, EX_CVDF3D => CVDF3D
00269 !
00270       IMPLICIT NONE
00271       INTEGER LNG,LU
00272       COMMON/INFO/LNG,LU
00273 !
00274 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00275 !
00276       TYPE(BIEF_OBJ), INTENT(INOUT)   :: FD, FC, FN
00277       TYPE(BIEF_OBJ), INTENT(INOUT)   :: S1F, VISCF
00278       TYPE(BIEF_OBJ), TARGET, INTENT(INOUT)   :: S0F
00279       TYPE(BIEF_OBJ), INTENT(INOUT)   :: LIFBOL, LIFBOF, LIFBOS
00280       TYPE(BIEF_OBJ), INTENT(INOUT)   :: FBORL, FBORF, FBORS
00281       TYPE(BIEF_OBJ), INTENT(IN)      :: AFBORL, AFBORF, AFBORS
00282       TYPE(BIEF_OBJ), INTENT(IN)      :: BFBORL, BFBORF, BFBORS
00283       TYPE(BIEF_OBJ), INTENT(IN)      :: FLUEXT,PLUIE,PARAPLUIE
00284       TYPE(BIEF_OBJ), INTENT(IN)      :: FLUEXTPAR
00285       DOUBLE PRECISION, INTENT(IN)    :: SIGMAF,FMIN,FMAX,DT,TRAIN
00286       DOUBLE PRECISION, INTENT(IN)    :: AGGLOD
00287       DOUBLE PRECISION, INTENT(INOUT) :: FLUXF,TETADI
00288       INTEGER, INTENT(IN)             :: SCHCF,SCHDF,TRBAF,NPTFR3,NFRLIQ
00289       INTEGER, INTENT(IN)             :: NUMLIQ(*),DIRFLU(0:NFRLIQ)
00290       LOGICAL, INTENT(IN)             :: CLIMIN,CLIMAX,RAIN,YAS0F,YAS1F
00291       LOGICAL, INTENT(IN)             :: INFOR,NEWDIF,CALFLU,MSK,SIGMAG
00292       TYPE(SLVCFG)                    :: SLVDIF
00293       TYPE(BIEF_OBJ), INTENT(IN)      :: MASKEL,IKLE3,FLODEL,FLOPAR
00294       TYPE(BIEF_OBJ), INTENT(IN)      :: NBOR3,WCC,SOURCES,ZPROP
00295       TYPE(BIEF_OBJ), INTENT(INOUT)   :: T3_01,T3_02,T3_03,T3_04,W1
00296       TYPE(BIEF_OBJ), INTENT(INOUT)   :: T2_01,T2_02,T2_03,ZT
00297       TYPE(BIEF_OBJ), TARGET, INTENT(INOUT) :: VOLUT
00298       TYPE(BIEF_MESH), INTENT(INOUT)  :: MESH3D
00299       INTEGER, INTENT(IN)             :: NPOIN3,NPOIN2,MAXADV
00300       INTEGER, INTENT(IN)             :: IPBOT(NPOIN2)
00301       INTEGER, INTENT(IN)             :: NPLAN,NELEM2,NELEM3,LV
00302       INTEGER, INTENT(IN)             :: OPTBAN,OPTDIF
00303       TYPE(BIEF_OBJ), INTENT(INOUT)   :: MMURD,MURD_TF,MTRA1
00304       TYPE(BIEF_OBJ), INTENT(IN)      :: VOLUN,VOLUNPAR,VOLUPAR
00305       TYPE(BIEF_OBJ), TARGET, INTENT(IN) :: VOLU
00306       LOGICAL, INTENT(IN)             :: INCHYD,YASEM3D
00307       LOGICAL, INTENT(INOUT)          :: YAWCC
00308       TYPE(BIEF_OBJ), INTENT(IN)      :: MASKPT,MASKBR,H,SVIDE
00309       TYPE(BIEF_MESH), INTENT(INOUT)  :: MESH2D
00310       INTEGER, INTENT(IN)             :: IELM3,IELM2H,IELM2V,NSCE
00311       TYPE(BIEF_OBJ), INTENT(INOUT)   :: SEM3D,IT1,IT2,TRAV3,MTRA2
00312       TYPE(BIEF_OBJ), INTENT(INOUT)   :: MSUPG,MDIFF,MATR2H
00313       DOUBLE PRECISION, INTENT(IN)    :: FSCE(NSCE)
00314       TYPE(BIEF_OBJ), INTENT(IN)      :: FLUDPT, VOLU2D, V2DPAR
00315       TYPE(BIEF_OBJ), INTENT(INOUT)   :: FLUDP, FLUER
00316       INTEGER, INTENT(IN)             :: SETDEP
00317 !
00318 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00319 !
00320       INTEGER IP,K,IPTFR,IS,I,IIS,PARA,DIM1X
00321       DOUBLE PRECISION, POINTER, DIMENSION(:) :: SAVEZ
00322       DOUBLE PRECISION STOFD,TETASUPG
00323       CHARACTER(LEN=1) :: S0FTYPR,NUZTYPR
00324       TYPE(BIEF_OBJ), POINTER :: VOLUME,S0F2
00325 !
00326 !     FUNCTIONS
00327 !
00328       DOUBLE PRECISION P_DSUM,LAMBDA
00329       EXTERNAL         P_DSUM
00330 !
00331       LOGICAL YADIRFLU,YASCE,VELOCITY,YARAIN
00332 !
00333 !***********************************************************************
00334 !
00335 !     FOR DIMENSIONING XA AND XB IN MURD3D.
00336 !
00337       DIM1X=BIEF_DIM2_EXT(IELM3,IELM3,1,'Q',MESH3D)
00338 !
00339 !     SAVING S0F%TYPR AND NUZ%TYPR
00340 !
00341       S0FTYPR=S0F%TYPR
00342       NUZTYPR=VISCF%ADR(3)%P%TYPR
00343 !
00344 !     DEALING WITH A VELOCITY ?
00345 !
00346       VELOCITY=.FALSE.
00347       IF(FN%NAME(1:1).EQ.'U'.OR.
00348      &   FN%NAME(1:1).EQ.'V'.OR.
00349      &   FN%NAME(1:1).EQ.'W') VELOCITY=.TRUE.
00350 !
00351 !     EVEN IF.NOT.CALFLU
00352 !
00353       FLUXF = 0.D0
00354 !
00355 !     WITH DISTRIBUTIVE SCHEMES : COMPUTES PRESCRIBED VALUES THAT
00356 !     WILL ENSURE THE CORRECT FLUX (REAL PRESCRIBED VALUES DISCARDED)
00357 !     THESE CORRECTED PRESCRIBED VALUES ARE SET BEFORE ADVECTION
00358 !
00359 !     YADIRFLU=.TRUE. : THERE IS AT LEAST ONE BOUNDARY WITH
00360 !                       TREATMENT OF FLUXES AT BOUNDARIES = 2
00361       YADIRFLU=.FALSE.
00362 !     DIRFLU DISCARDED FOR VELOCITIES
00363       IF(NFRLIQ.GT.0.AND..NOT.VELOCITY) THEN
00364         DO K=1,NFRLIQ
00365           IF(DIRFLU(K).EQ.2) YADIRFLU=.TRUE.
00366         ENDDO
00367       ENDIF
00368 !
00369 !=======================================================================
00370 !
00371 !     FOR TRACERS (=NOT VELOCITY) : DIRICHLET VALUES ARE NOT RESPECTED IF EXIT
00372 !     THERE IS NO NEED TO TEST KENTU OR KADH FOR TRACERS
00373 !
00374       IF(NPTFR3.GT.0.AND.NFRLIQ.GT.0.AND..NOT.VELOCITY) THEN
00375         DO IPTFR=1,NPTFR3
00376           IF(LIFBOL%I(IPTFR).EQ.KENT) THEN
00377 !           EXITS ARE TREATED AS FREE BOUNDARIES
00378             IP=NBOR3%I(IPTFR)
00379             IF(FLUEXTPAR%R(IP).GE.0.D0) LIFBOL%I(IPTFR)=KSORT
00380           ENDIF
00381         ENDDO
00382       ENDIF
00383 !
00384 !=======================================================================
00385 !
00386 !     A PRIORI CORRECTION OF FN FOR REAL ENTRANCES
00387 !     I.E. LIFBOL STILL KENT DESPITE ABOVE CHANGE
00388 !
00389       IF((SCHCF.EQ.ADV_SUP   .OR.SCHCF.EQ.ADV_NSC    .OR.
00390      &    SCHCF.EQ.ADV_PSI   .OR.SCHCF.EQ.ADV_LPO    .OR.
00391      &    SCHCF.EQ.ADV_LPO_TF.OR.SCHCF.EQ.ADV_NSC_TF)
00392      &                                              .AND.YADIRFLU) THEN
00393 !
00394         IF(NPTFR3.GT.0) THEN
00395 !
00396         DO IP=1,NPTFR3
00397           IF(DIRFLU(NUMLIQ(IP)).EQ.2.AND.LIFBOL%I(IP).EQ.KENT) THEN
00398             I=NBOR3%I(IP)
00399             LAMBDA=-FLUEXTPAR%R(I)*DT/
00400      &      (MAX(VOLUNPAR%R(I),1.D-10)-FLUEXTPAR%R(I)*DT)
00401             FN%R(I)=FN%R(I)+LAMBDA*(FBORL%R(IP)-FN%R(I))
00402 !           CORRECTION OF FLUX
00403 !           IN THE PROOF OF MASS-CONSERVATION, FLUEXT IS MULTIPLIED
00404 !           BY FN INSTEAD OF FBOR, TO INTERPRET THE ADDED MASS AS
00405 !           A FLUX THIS CORRECTION IS NECESSARY
00406 !           HERE IT IS THE FN MODIFIED ABOVE
00407 !           EVEN IF NOT CALFLU (CHEAPER)
00408             FLUXF=FLUXF+(FBORL%R(IP)-FN%R(I))*FLUEXT%R(I)*DT
00409 !           AVOIDS A DIRICHLET TREATMENT HEREAFTER AND BY DIFF3D -
00410 !           WILL BE RESTORED AFTER DIFF3D
00411             LIFBOL%I(IP)=KSORT
00412           ENDIF
00413         ENDDO
00414 !
00415         ENDIF
00416 !
00417       ENDIF
00418 !
00419 !=======================================================================
00420 !
00421 !     PUTS DIRICHLET VALUES IN FN
00422 !     MAY HAVE NO EFFECT IF TREATMENT OF FLUXES AT THE BOUNDARIES=2
00423 !     BECAUSE LIFBOL CHANGED ABOVE
00424 !
00425       IF(NPTFR3.GT.0) THEN
00426         DO IPTFR=1,NPTFR3
00427           IF(LIFBOL%I(IPTFR).EQ.KENT .OR.
00428      &       LIFBOL%I(IPTFR).EQ.KENTU.OR.
00429      &       LIFBOL%I(IPTFR).EQ.KADH) THEN
00430              FN%R(NBOR3%I(IPTFR)) = FBORL%R(IPTFR)
00431           ENDIF
00432         ENDDO
00433       ENDIF
00434 !
00435 !     HERE BOTTOM AND FREE SURFACE SHOULD BE TREATED AS WELL
00436 !
00437 !=======================================================================
00438 !
00439 !     3D ADVECTION (OTHER THAN SUPG)
00440 !
00441 !=======================================================================
00442 !
00443 !     WITH DISTRIBUTIVE SCHEMES, RIGHT-HAND SIDE MUST BE
00444 !     IN INTEGRATED FORM
00445 !
00446       IF(SCHCF.EQ.ADV_NSC.OR.SCHCF.EQ.ADV_PSI.OR.SCHCF.EQ.ADV_LPO.OR.
00447      &   SCHCF.EQ.ADV_NSC_TF.OR.SCHCF.EQ.ADV_LPO_TF) THEN
00448 !
00449 !       S0F2 WILL BE MASS-MATRIX * S0F ASSEMBLED IN PARALLEL
00450 !
00451         S0F2=>TRAV3%ADR(10)%P
00452         S0F2%TYPR=S0F%TYPR
00453 !
00454         IF(S0F%TYPR.NE.'0') THEN
00455 !
00456           CALL VECTOR(S0F2,'=','MASVEC          ',IELM3,1.D0,
00457      &                S0F,S0F,S0F,S0F,S0F,S0F,MESH3D,MSK,MASKEL)
00458           IF(NCSIZE.GT.1) CALL PARCOM(S0F2,2,MESH3D)
00459 !
00460         ENDIF
00461 !
00462       ENDIF
00463 !
00464 !-----------------------------------------------------------------------
00465 !
00466 !     ADVECTION BY CHARACTERISTICS
00467 !
00468       IF(SCHCF.EQ.ADV_CAR) THEN
00469 !
00470 !       THIS IS NOW DONE IN CHARAC CALLED BY PRECON
00471 !
00472 !       CALL CARA3D(FC%R,FN%R,SHP%R,SHZ%R,ELT%I,ETA%I,IKLE2%I,
00473 !    &              NELEM2,NPOIN2,NPOIN3,DT,INFOR)
00474 !       IF(NCSIZE.GT.1) CALL PARCOM(FC,2,MESH3D)
00475 !
00476 !-----------------------------------------------------------------------
00477 !
00478 !     ADVECTION BY MURD DISTRIBUTIVE SCHEME, OPTION N OR PSI
00479 !
00480       ELSEIF(SCHCF.EQ.ADV_NSC.OR.SCHCF.EQ.ADV_PSI) THEN
00481 !
00482         CALL MURD3D(FC%R,FN%R,VOLU%R,VOLUN%R,T3_01%R,T3_01,
00483      &              MMURD%D%R,MMURD%X%R,DIM1X,
00484      &              T3_02%R,T3_03%R,T3_04%R,T3_02,T3_03,T3_04,
00485      &              IKLE3%I,MESH3D,
00486      &              NELEM3,NPOIN3,DT,SCHCF,LV,MSK,MASKEL%R,INFOR,
00487      &              CALFLU,FLUXF,FLUEXT%R,S0F2,NSCE,SOURCES,FSCE,
00488      &              RAIN,PARAPLUIE%R,TRAIN,NPOIN2,
00489      &              TRAV3%ADR(5)%P,TRAV3%ADR(6)%P,MASKPT%R,OPTBAN,
00490      &              FLODEL%R,FLOPAR%R,MESH3D%GLOSEG%I,
00491      &              MESH3D%GLOSEG%DIM1,MESH2D%NSEG,NPLAN,IELM3)
00492 !
00493 !       S0F CANCELLED TO AVOID A DUPLICATE TREATMENT
00494 !       IF DIFF3D IS CALLED AFTER
00495         S0F%TYPR='0'
00496 !
00497 !-----------------------------------------------------------------------
00498 !
00499 !     ADVECTION BY UPWIND EXPLICIT FINITE VOLUME SCHEME
00500 !
00501       ELSEIF(SCHCF.EQ.ADV_LPO) THEN
00502 !
00503         CALL MURD3D(FC%R,FN%R,VOLU%R,VOLUN%R,T3_01%R,T3_01,
00504      &              MESH3D%M%D%R,MESH3D%M%X%R,DIM1X,
00505      &              T3_02%R,T3_03%R,T3_04%R,T3_02,T3_03,T3_04,
00506      &              IKLE3%I,MESH3D,
00507      &              NELEM3,NPOIN3,DT,SCHCF,LV,MSK,MASKEL%R,INFOR,
00508      &              CALFLU,FLUXF,FLUEXT%R,S0F2,NSCE,SOURCES,FSCE,
00509      &              RAIN,PARAPLUIE%R,TRAIN,NPOIN2,
00510      &              TRAV3%ADR(5)%P,TRAV3%ADR(6)%P,MASKPT%R,OPTBAN,
00511      &              FLODEL%R,FLOPAR%R,MESH3D%GLOSEG%I,
00512      &              MESH3D%GLOSEG%DIM1,MESH2D%NSEG,NPLAN,IELM3)
00513 !
00514 !       S0F CANCELLED TO AVOID A DUPLICATE TREATMENT
00515 !       IF DIFF3D IS CALLED AFTER
00516         S0F%TYPR='0'
00517 !
00518 !-----------------------------------------------------------------------
00519 !
00520 !     ADVECTION BY UPWIND EXPLICIT FINITE VOLUME SCHEME
00521 !
00522       ELSEIF(SCHCF.EQ.ADV_LPO_TF) THEN
00523 !
00524         CALL MURD3D_POS(FC%R,FN%R,VOLU%R,VOLU,VOLUN%R,VOLUN,
00525      &                  T3_01%R,T3_01,MESH3D%M%X%R,
00526      &                  T3_02%R,T3_03%R,T3_04%R,T3_02,T3_03,T3_04,
00527      &                  MESH2D,MESH3D,
00528      &                  NELEM3,NPOIN3,DT,SCHCF,MSK,MASKEL%R,INFOR,
00529      &                  CALFLU,FLUXF,FLUEXT%R,S0F2,NSCE,SOURCES,FSCE,
00530      &                  RAIN,PARAPLUIE%R,TRAIN,NPOIN2,OPTBAN,
00531      &                  FLODEL%R,FLOPAR%R,MESH3D%GLOSEG%I,
00532      &                  MESH3D%GLOSEG%DIM1,MESH2D%NSEG,NPLAN,
00533      &                  TRAV3%ADR(6)%P,TRAV3%ADR(7)%P,
00534      &                  TRAV3%ADR(8)%P,
00535      &                  TRAV3%ADR(9)%P,2,IELM3,MAXADV)
00536 !
00537 !       S0F CANCELLED TO AVOID A DUPLICATE TREATMENT
00538 !       IF DIFF3D IS CALLED AFTER
00539         S0F%TYPR='0'
00540 !
00541 !-----------------------------------------------------------------------
00542 !
00543 !     ADVECTION BY UPWIND EXPLICIT FINITE VOLUME SCHEME
00544 !
00545       ELSEIF(SCHCF.EQ.ADV_NSC_TF) THEN
00546 !
00547         PARA=0
00548         IF(NCSIZE.GT.1) PARA=MESH3D%NSEG
00549         CALL MURD3D_POS(FC%R,FN%R,VOLU%R,VOLU,VOLUN%R,VOLUN,
00550      &                  T3_01%R,T3_01,MESH3D%M%X%R,
00551      &                  T3_02%R,T3_03%R,T3_04%R,T3_02,T3_03,T3_04,
00552      &                  MESH2D,MESH3D,
00553      &                  NELEM3,NPOIN3,DT,SCHCF,MSK,MASKEL%R,INFOR,
00554      &                  CALFLU,FLUXF,FLUEXT%R,S0F2,NSCE,SOURCES,FSCE,
00555      &                  RAIN,PARAPLUIE%R,TRAIN,NPOIN2,OPTBAN,
00556      &                  MURD_TF%X%R(1     :MESH3D%NSEG     ),
00557      &                  MURD_TF%X%R(1+PARA:MESH3D%NSEG+PARA),
00558      &                  MESH3D%GLOSEG%I,
00559      &                  MESH3D%GLOSEG%DIM1,MESH2D%NSEG,NPLAN,
00560      &                  TRAV3%ADR(6)%P,TRAV3%ADR(7)%P,
00561      &                  TRAV3%ADR(8)%P,
00562      &                  TRAV3%ADR(9)%P,2,IELM3,MAXADV)
00563 !
00564 !       S0F CANCELLED TO AVOID A DUPLICATE TREATMENT
00565 !       IF DIFF3D IS CALLED AFTER
00566         S0F%TYPR='0'
00567 !
00568 !-----------------------------------------------------------------------
00569 !
00570 !     OTHER CASES (SUPG OR NO ADVECTION)
00571 !
00572       ELSE
00573 !
00574         CALL OS ( 'X=Y     ' , X=FC , Y=FN )
00575 !
00576       ENDIF
00577 !
00578 !-----------------------------------------------------------------------
00579 !
00580 !     IMPLICITT VERTICAL DIFFUSION SCHEME FOR SEDIMENTS ADDED HERE
00581 !
00582       IF(YAWCC.AND.SETDEP.EQ.1) THEN
00583 
00584         CALL SET_DIF(FC%R,VOLU2D,MESH3D%Z%R,NPOIN2,NPOIN3,DT,FLUXF,
00585      &             NPLAN,WCC,FLUDPT,FLUDP,FLUER,IPBOT,VISCF%ADR(3)%P)
00586 !
00587 !       VERTICAL DIFFUSION HAS ALREADY BEEN TREATED THEN IT IS SET
00588 !       TO 0 (VIRTUALLY, WITH COMPONENT TYPR, SEE HOW MT02PP IS WRITTEN).
00589 !       IDEM WITH SETTLING VELOCITY
00590 !
00591         VISCF%ADR(3)%P%TYPR='0'
00592         YAWCC=.FALSE.
00593 !
00594       ENDIF
00595 !
00596 !-----------------------------------------------------------------------
00597 !
00598 !     RE-ENFORCES DIRICHLET POINTS (MAY CAUSE MASS ERRORS)
00599 !     IN FACT NOT DONE IF LIFBOL HAS BEEN CHANGED ABOVE INTO KSORT
00600 !     HENCE NO EFFECT WHEN YADIRFLU=.TRUE.
00601 !
00602       IF(NPTFR3.GT.0) THEN
00603       DO IP=1,NPTFR3
00604         IF(LIFBOL%I(IP).EQ.KENT .OR.
00605      &     LIFBOL%I(IP).EQ.KENTU.OR.
00606      &     LIFBOL%I(IP).EQ.KADH) THEN
00607            I=NBOR3%I(IP)
00608            FC%R(I) = FBORL%R(IP)
00609         ENDIF
00610       ENDDO
00611       ENDIF
00612 !
00613 !     BOTTOM AND FREE SURFACE
00614 !
00615       K = NPOIN3 - NPOIN2
00616       DO IP = 1,NPOIN2
00617         IF(LIFBOF%I(IP).EQ.KENT.OR.LIFBOF%I(IP).EQ.KADH) THEN
00618           FC%R(IP)   = FBORF%R(IP)
00619         ENDIF
00620         IF(LIFBOS%I(IP).EQ.KENT.OR.LIFBOS%I(IP).EQ.KADH) THEN
00621           FC%R(IP+K) = FBORS%R(IP)
00622         ENDIF
00623       ENDDO
00624 !
00625 !=======================================================================
00626 !
00627 !  SUPG ADVECTION AND/OR DIFFUSION
00628 ! (IN THIS CASE IT IS NECESSARY TO SOLVE A LINEAR SYSTEM)
00629 !
00630 !=======================================================================
00631 !
00632       IF(SCHCF.EQ.ADV_SUP.OR.SCHDF.NE.0) THEN
00633 !
00634         IF(SCHCF.EQ.ADV_SUP) THEN
00635 !         IT SEEMS THAT WITH SUBITERATIONS ONLY TETASUPG=1-TETAH WORKS
00636 !         FOR MASS-CONSERVATION. SEE ALSO DIFF3D WITH ANOTHER COMMENT
00637           TETASUPG=0.55D0
00638         ELSE
00639           TETASUPG=1.D0
00640         ENDIF
00641 !
00642         IF(SCHCF.EQ.ADV_SUP.AND..NOT.VELOCITY) THEN
00643           CALL OS('X=CY    ',X=VOLUT,Y=VOLUN    ,C=     TETASUPG)
00644           CALL OS('X=X+CY  ',X=VOLUT,Y=VOLU     ,C=1.D0-TETASUPG)
00645           CALL OS('X=CY    ',X=ZT   ,Y=ZPROP    ,C=     TETASUPG)
00646           CALL OS('X=X+CY  ',X=ZT   ,Y=MESH3D%Z ,C=1.D0-TETASUPG)
00647 !         ZT IS TEMPORARILY PUT IN MESH3D%Z
00648           SAVEZ=>MESH3D%Z%R
00649           MESH3D%Z%R=>ZT%R
00650           VOLUME=>VOLUT
00651         ELSE
00652           VOLUME=>VOLU
00653         ENDIF
00654 !
00655         IF(SCHCF.EQ.ADV_CAR.OR.SCHCF.EQ.ADV_SUP) THEN
00656 !         SOURCES HAVE TO BE TREATED
00657           YASCE=.TRUE.
00658           YARAIN=RAIN
00659         ELSE
00660 !         SOURCES HAVE ALREADY BEEN TREATED BY DISTRIBUTIVE SCHEMES
00661           YASCE=.FALSE.
00662 !         RAIN HAS ALREADY BEEN TREATED BY DISTRIBUTIVE SCHEMES
00663           YARAIN=.FALSE.
00664         ENDIF
00665 !
00666         CALL DIFF3D(FD,FC,FN,VISCF,SIGMAF,
00667      &              S0F,YAS0F,S1F,YAS1F,
00668      &              FBORL,FBORF,FBORS,AFBORL,AFBORF,AFBORS,
00669      &              BFBORL,BFBORF,BFBORS,LIFBOF,LIFBOL,LIFBOS,
00670      &              FMIN,CLIMIN,FMAX,CLIMAX,
00671      &              SCHCF,SCHDF,SLVDIF,TRBAF,INFOR,NEWDIF,
00672      &              DT,T2_01,T2_02,T2_03,T3_01,T3_02,T3_03,T3_04,
00673      &              NPOIN2,NPOIN3,INCHYD,SEM3D,YASEM3D,IT1,
00674      &              NPTFR3,NBOR3,MASKPT,TRAV3,MESH2D,
00675      &              MESH3D,MTRA1,MTRA2,IELM3,MSUPG,IELM2H,IELM2V,
00676      &              MDIFF,MATR2H,MASKBR,SVIDE,MSK,MASKEL,H,
00677      &              NPLAN,OPTBAN,OPTDIF,TETADI,YAWCC,WCC,AGGLOD,
00678      &              VOLUME,YASCE,NSCE,FSCE,SOURCES,TETASUPG,
00679      &              VELOCITY,YARAIN,PLUIE%R,TRAIN,SIGMAG,IPBOT,
00680      &              SETDEP)
00681 !
00682         IF(SCHCF.EQ.ADV_SUP.AND..NOT.VELOCITY) THEN
00683 !         MESH3D%Z RESTORED
00684           MESH3D%Z%R=>SAVEZ
00685         ENDIF
00686 !
00687       ELSE
00688         CALL OS ( 'X=Y     ', X=FD, Y=FC )
00689       ENDIF
00690 !
00691 !-----------------------------------------------------------------------
00692 !
00693 !     ADVECTIVE FLUXES AND SOURCES
00694 !
00695       IF(CALFLU) THEN
00696 !
00697         IF(SCHCF.EQ.ADV_CAR) THEN
00698           DO IP = 1,NPOIN3
00699             FLUXF = FLUXF + FN%R(IP)*FLUEXT%R(IP)*DT
00700           ENDDO
00701         ELSEIF(SCHCF.EQ.ADV_SUP) THEN
00702           DO IP = 1,NPOIN3
00703             FLUXF = FLUXF + FLUEXT%R(IP)*DT*
00704      &                      (TETASUPG*FD%R(IP)+(1.D0-TETASUPG)*FN%R(IP))
00705           ENDDO
00706         ENDIF
00707 !
00708 !       CHARACTERISTICS OR SUPG : FLUX DUE TO SOURCES
00709 !       (FOR DISTRIBUTIVE SCHEMES IT IS DONE IN MURD3D)
00710 !
00711         IF(NSCE.GT.0.AND.(SCHCF.EQ.ADV_CAR.OR.SCHCF.EQ.ADV_SUP)) THEN
00712           DO IS=1,NSCE
00713             IIS=IS
00714 !           HERE IN PARALLEL SOURCES WITHOUT PARCOM
00715             IF(NCSIZE.GT.1) IIS=IIS+NSCE
00716             DO IP=1,NPOIN3
00717               IF(SOURCES%ADR(IS)%P%R(IP).GT.0.D0) THEN
00718                 FLUXF=FLUXF-FSCE(IS)*SOURCES%ADR(IIS)%P%R(IP)*DT
00719               ELSE
00720 !                           FN FOR CHARACTERISTICS ?
00721                 FLUXF=FLUXF-FD%R(IP)*SOURCES%ADR(IIS)%P%R(IP)*DT
00722               ENDIF
00723             ENDDO
00724           ENDDO
00725         ENDIF
00726 !
00727       ENDIF
00728 !
00729 !-----------------------------------------------------------------------
00730 !
00731 !     A POSTERIORI CORRECTION OF SUPG RESULTS
00732 !
00733       IF(SCHCF.EQ.ADV_SUP.AND.YADIRFLU) THEN
00734 !
00735 !       CORRECTED VALUE AND CORRESPONDING FLUX CORRECTION
00736 !
00737         IF(NPTFR3.GT.0) THEN
00738         DO IP=1,NPTFR3
00739           IF(DIRFLU(NUMLIQ(IP)).EQ.2) THEN
00740             IF(LIFBOL%I(IP+NPTFR3).EQ.KENT .OR.
00741      &        LIFBOL%I(IP+NPTFR3).EQ.KENTU.OR.
00742      &        LIFBOL%I(IP+NPTFR3).EQ.KADH     ) THEN
00743 !             ONLY ENTRANCES
00744               I=NBOR3%I(IP)
00745               IF(FLUEXTPAR%R(I).LT.0.D0) THEN
00746                 STOFD=FD%R(I)
00747                 LAMBDA=-FLUEXTPAR%R(I)*TETASUPG*DT/
00748      &                  MAX(VOLUPAR%R(I),1.D-10)
00749                 FD%R(I)=STOFD+LAMBDA*(FN%R(I)-STOFD)
00750 !               CORRECTION OF FLUX
00751 !               A POSTERIORI ADDED MASS DUE TO CORRECTION
00752                 FLUXF=FLUXF-VOLU%R(I)*(FD%R(I)-STOFD)
00753               ENDIF
00754             ENDIF
00755           ENDIF
00756         ENDDO
00757         ENDIF
00758 !
00759       ENDIF
00760 !
00761 !-----------------------------------------------------------------------
00762 !
00763       IF(CALFLU) THEN
00764 !       NOW RETURNS TO REAL FLUXES, NOT FLUXES*DT
00765         FLUXF = FLUXF / DT
00766 !       PARALLEL MODE
00767         IF(NCSIZE.GT.1) FLUXF = P_DSUM(FLUXF)
00768       ENDIF
00769 !
00770 !-----------------------------------------------------------------------
00771 !
00772 !     RESTORES ORIGINAL LIFBOL FROM SECOND DIMENSION
00773 !
00774       DO IP=1,NPTFR3
00775         LIFBOL%I(IP)=LIFBOL%I(IP+NPTFR3)
00776       ENDDO
00777 !
00778 !     S0F%TYPR RESTORED (S0F MAY BE USED ELSEWHERE, E.G. IN A FURTHER
00779 !                        SUB ITERATION)
00780 !
00781       S0F%TYPR=S0FTYPR
00782 !
00783 !     RESTORING NUZ FOR FURTHER USES
00784 !
00785       VISCF%ADR(3)%P%TYPR=NUZTYPR
00786 !
00787 !-----------------------------------------------------------------------
00788 !
00789       RETURN
00790       END

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