flux3d.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac3d\flux3d.f
00002 !
00089                      SUBROUTINE FLUX3D
00090 !                    *****************
00091 !
00092      &(FLUINT,FLUEXT,FLUEXTPAR,UCONV,VCONV,TRA01,TRA02,TRA03,
00093      & W1,NETAGE,NPLAN,NELEM3,IELM3,IELM2H,IELM2V,SVIDE,MESH3,
00094      & MSK,MASKEL,MASK_3D,LIHBOR,KENT,NPTFR,DT,VOLU,VOLUN,
00095      & MESH2,GRAPRD,SIGMAG,TRAV2,NPOIN2,NPOIN3,DM1,GRAZCO,
00096      & FLBOR,PLUIE,RAIN,FLODEL,OPT_HNEG,FLULIM,YACVVF,LT,BYPASS,
00097      & N_ADV,WEL)
00098 !
00099 !***********************************************************************
00100 ! TELEMAC3D   V6P2                                   21/08/2010
00101 !***********************************************************************
00102 !
00103 !
00104 !
00105 !
00106 !
00107 !
00108 !
00109 !
00110 !
00111 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00112 !| BYPASS         |---| IF YES, BYPASS VOID VOLUMES
00113 !| DM1            |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
00114 !|                |   | IS DM1*GRAD(ZCONV), SEE SOLSYS. GRAD(ZCONV) IS
00115 !|                |   | HERE GRAZCO.
00116 !| DT             |-->| TIME STEP
00117 !| FLBOR          |<->| FLUXES AT BOUNDARIES
00118 !| FLODEL         |<->| FLUXES BY SEGMENT
00119 !| FLUEXT         |<--| OUTPUT FLUX BY NODE
00120 !| FLUEXTPAR      |<--| OUTPUT FLUX BY NODE, IN PARALLEL
00121 !| FLUINT         |<--| INPUT FLUX BY NODE
00122 !| FLULIM         |<->| LIMITATION OF FLUXES
00123 !| GRAPRD         |-->| GRAPHIC PRINTOUT PERIOD
00124 !| GRAZCO         |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
00125 !|                |   | IS DM1*GRAD(ZCONV), GRAZCO=GRAD(ZCONV)
00126 !| IELM2H         |-->| TYPE OF ELEMENT
00127 !| IELM2V         |-->| TYPE DE DISCRETISATION 2DV
00128 !| IELM3          |-->| TYPE DE DISCRETISATION 3D
00129 !| KENT           |-->| CONVENTION FOR PRESCRIBED DEPTH
00130 !| LIHBOR         |-->| BOUNDARY CONDITIONS ON DEPTH
00131 !| LT             |-->| CURRENT TIME STEP N
00132 !| MASK_3D        |-->| 3D MASKS ON LATERAL BOUNDARIES
00133 !| MASKEL         |<->| MASKING OF 3D ELEMENTS
00134 !|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
00135 !| MESH2          |<->| 2D MESH
00136 !| MESH3          |<->| 3D MESH
00137 !| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
00138 !| NELEM3         |-->| NUMBER OF ELEMENTS IN 3D
00139 !| NETAGE         |-->| NUMBER OF LAYERS I.E. NUMBER OF PLANES - 1
00140 !| NPLAN          |-->| NUMBER OF PLANES IN THE 3D MESH OF PRISMS
00141 !| NPOIN2         |-->| NUMBER OF 2D POINTS
00142 !| NPOIN3         |-->| NUMBER OF 3D POINTS
00143 !| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
00144 !| OPT_HNEG       |---| OPTION FOR NEGATIVE DEPTHS
00145 !| PLUIE          |-->| RAIN IN M/S MULTIPLIED BY VOLU2D
00146 !| RAIN           |-->| IF YES, THERE IS RAIN OR EVAPORATION
00147 !| SIGMAG         |-->| LOGICAL FOR GENERALISED SIGMA TRANSFORMATION
00148 !| SVIDE          |<->| VOID STRUCTURE
00149 !| TRA01          |<->| WORK ARRAY
00150 !| TRA02          |<->| WORK ARRAY
00151 !| TRA03          |<->| WORK ARRAY
00152 !| TRAV2          |<->| WORK ARRAY
00153 !| VCONV          |-->| ADVECTION VELOCITY FIELD
00154 !| VOLU           |-->| VOLUME AROUND POINTS AT TIME N+1
00155 !| VOLUN          |-->| VOLUME AROUND POINTS AT TIME N
00156 !| W1             |<->| WORK ARRAY IN 3D
00157 !| WEL            |<->| BIEF WORK ARRAY
00158 !| YACVVF         |-->| THERE IS AN ADVECTION WITH FINITE VOLUMES
00159 !|                |   | (HENCE COMPUTATION OF FLUXES REQUIRED)
00160 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00161 !
00162       USE BIEF
00163       USE DECLARATIONS_TELEMAC,   ONLY : ADV_NSC,ADV_PSI,ADV_NSC_TF
00164       USE DECLARATIONS_TELEMAC3D, ONLY : LV
00165 !
00166       IMPLICIT NONE
00167       INTEGER LNG,LU
00168       COMMON/INFO/LNG,LU
00169 !
00170 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00171 !
00172       INTEGER, INTENT(IN) :: NETAGE,NPLAN,NELEM3,NPOIN2,NPOIN3,LT
00173       INTEGER, INTENT(IN) :: IELM3,IELM2H,IELM2V,OPT_HNEG
00174       INTEGER, INTENT(IN) :: KENT,NPTFR,GRAPRD
00175       INTEGER, INTENT(IN) :: LIHBOR(NPTFR),N_ADV(0:15)
00176 !
00177       TYPE(BIEF_OBJ), INTENT(INOUT) :: FLUINT,FLUEXT,FLUEXTPAR
00178       TYPE(BIEF_OBJ), INTENT(IN)    :: VOLU,VOLUN,DM1,GRAZCO
00179       TYPE(BIEF_OBJ), INTENT(IN)    :: UCONV,VCONV,PLUIE,MASK_3D
00180       TYPE(BIEF_OBJ), INTENT(INOUT) :: MASKEL,WEL
00181       TYPE(BIEF_OBJ), INTENT(INOUT) :: FLULIM
00182       TYPE(BIEF_OBJ), INTENT(INOUT), TARGET :: FLODEL
00183       TYPE(BIEF_OBJ), INTENT(INOUT) :: FLBOR
00184       TYPE(BIEF_OBJ), INTENT(INOUT) :: SVIDE,TRA01,TRA02,TRA03,W1,TRAV2
00185       TYPE(BIEF_MESH), INTENT(INOUT):: MESH3, MESH2
00186 !
00187       DOUBLE PRECISION, INTENT(IN)  :: DT
00188       LOGICAL, INTENT(IN)           :: MSK,SIGMAG,RAIN,YACVVF,BYPASS
00189 !
00190 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00191 !
00192       INTEGER IPLAN,IPTFR,I,IOPT
00193       CHARACTER(LEN=16) FORMUL
00194 !     DOUBLE PRECISION SUM_FLUEXT
00195 !
00196 !***********************************************************************
00197 !
00198       CALL OS('X=0     ',X=FLUEXT)
00199 !
00200 !=======================================================================
00201 !
00202 !   INTERNAL ADVECTION FLUXES
00203 !
00204 !=======================================================================
00205 !
00206 !        /            D(PSII*)           D(PSII*)
00207 !       /     H * U * -------- + H * V * --------   D(OMEGA*)
00208 !      /OMEGA*           DX                 DY
00209 !
00210 !
00211       FORMUL = 'VGRADP 2     HOR'
00212       CALL VECTOR(FLUINT,'=',FORMUL,IELM3,1.D0,DM1,GRAZCO,GRAZCO,
00213      &            UCONV,VCONV,SVIDE,MESH3,MSK,MASKEL)
00214 !
00215 !     STORING NON-ASSEMBLED FLUINT IN WEL
00216 !     WHICH IS NOT TO BE USED BEFORE BUILDING MATRICES MMURD OR
00217 !     MURD_TF (SEE PRECON AND MT14PP)
00218 !
00219       IF(N_ADV(ADV_NSC).GT.0.OR.N_ADV(ADV_PSI   ).GT.0
00220      &                      .OR.N_ADV(ADV_NSC_TF).GT.0) THEN
00221         IF(IELM3.EQ.41) THEN
00222           DO I=1,6*MESH3%NELEM
00223             WEL%R(I)=MESH3%W%R(I)
00224           ENDDO
00225           IF(OPT_HNEG.EQ.2) THEN
00226 !           LIMITATION OF NON-ASSEMBLED FLUXES STORED IN WEL
00227             CALL NA_FLUX3D_LIM(WEL%R,FLULIM%R,
00228      &                         MESH2%NSEG,
00229      &                         MESH3%NELEM,MESH3%NELMAX,
00230      &                         MESH2%NELEM,MESH2%NELMAX,
00231      &                         MESH2%ELTSEG%I,MESH2%ORISEG%I)
00232 !           REASSEMBLING FLUINT
00233             CALL ASSVEC(FLUINT%R,MESH3%IKLE%I,NPOIN3,MESH3%NELEM,
00234      &                MESH3%NELMAX,IELM3,WEL%R,.TRUE.,
00235      &                LV,MSK,MASKEL%R,6)
00236           ENDIF
00237         ELSEIF(IELM3.EQ.51) THEN
00238           DO I=1,4*MESH3%NELEM
00239             WEL%R(I)=MESH3%W%R(I)
00240           ENDDO
00241         ELSE
00242           WRITE(LU,*) 'FLUX3D: ELEMENT ',IELM3,' NOT IMPLEMENTED'
00243           CALL PLANTE(1)
00244           STOP
00245         ENDIF
00246       ENDIF
00247 !
00248 !     COMPUTING POINT TO POINT FLUXES:
00249 !     FOR ADVECTION SCHEMES ADV_LEO OR ADV_LEO_TF, TO HAVE FLODEL
00250 !     OR IN CASE THERE IS A LIMITATION WITH FLULIM (OPT_HNEG=2) BECAUSE
00251 !     IN THAT CASE A NEW FLUINT CONSISTENT WITH THE NEW CONTINUITY
00252 !     EQUATION MUST BE BUILT
00253 !
00254 !     HERE THE CONVENTION FOR SEGMENTS, DUE TO THE CHOICE OF FLUINT, IS
00255 !     THAT A SEGMENT WITH POSITIVE FLUX IS MEANT WITH A FLOW FROM POINT 2
00256 !     TO POINT 1.
00257 !
00258 !     NOT YET DONE FOR TETRAHEDRA
00259       IF(IELM3.EQ.41) THEN
00260 !
00261       IF(YACVVF) THEN
00262 !
00263         IOPT=2
00264 !       BEWARE, WITH IELM3=51 RETURNS FLODEL IN 2D
00265 !       BUT ADV_LEO AND ADV_LEO_PT NEVER DONE WITH IELM3=51
00266         CALL FLUX_EF_VF_3D(FLODEL%R,MESH2%W%R,MESH3%W%R,
00267      &                     MESH2%NSEG,MESH3%NSEG,MESH2%NELEM,
00268      &                     MESH3%NELEM,MESH2,.TRUE.,IOPT,1,
00269      &                     IELM3,NPLAN,MESH3%IKLE%I,MESH3%NELMAX,
00270      &                     MESH3%KNOLG%I)
00271 !
00272         IF(OPT_HNEG.EQ.2) THEN
00273 !         LIMITATION OF 3D FLUXES WITH 2D LIMITATIONS
00274           CALL FLUX3DLIM(FLODEL%R,FLULIM%R,NPLAN,MESH2%NSEG,NPOIN2,1)
00275 !         NEW ASSEMBLY OF FLUINT (IN THIS CASE ASSEMBLING FLUINT
00276 !                                 IN VECTOR ABOVE IS USELESS)
00277           CALL ASSEG_3D(FLODEL%R,FLUINT%R,NPOIN3,NPLAN,MESH2%NSEG,
00278      &                  MESH3%GLOSEG%I,MESH3%GLOSEG%DIM1,.TRUE.)
00279         ENDIF
00280 !
00281       ENDIF
00282 !
00283       ENDIF
00284 !
00285 !=======================================================================
00286 !
00287 !   COMPUTES THE ADVECTIVE FLUXES ON THE LATERAL LIQUID BOUNDARIES
00288 !
00289 !=======================================================================
00290 !
00291 !     /        ->  ->
00292 !    /     H * U . N  PSII*  D(OMEGA*)
00293 !   /
00294 !  /LIQUID BOUNDARIES*
00295 !
00296       FORMUL = 'FLUBOR          '
00297 !
00298 !     MASK_3D%ADR(8)%P: MASK ON LIQUID LATERAL BOUNDARIES
00299 !
00300       CALL VECTOR
00301      & (TRA02, '=', FORMUL, IELM2V, 1.D0, SVIDE, SVIDE, SVIDE,
00302      &  UCONV, VCONV, SVIDE, MESH3,.TRUE.,MASK_3D%ADR(8)%P)
00303 !
00304       CALL OSDB( 'X=Y     ' , FLUEXT , TRA02 , TRA02 , 0.D0 , MESH3 )
00305 !
00306 !
00307 !     IF(RAIN) THEN
00308 !       DO I=1,NPOIN2
00309 !         FLUEXT%R((NPLAN-1)*NPOIN2+I)=
00310 !    &    FLUEXT%R((NPLAN-1)*NPOIN2+I)+PLUIE%R(I)
00311 !       ENDDO
00312 !     ENDIF
00313 !
00314 !-----------------------------------------------------------------------
00315 !
00316 !     COMPUTATION OF FLUEXT ON POINTS WITH PRESCRIBED DEPTH
00317 !     SO THAT CONTINUITY IS ENSURED.
00318 !
00319 !     EXCEPT AT THE FIRST CALL (BY THE FIRST CALL TO PRECON) FLBOR
00320 !     HAS ALREADY BEEN COMPUTED (WAVE_EQUATION AND POSSIBLY
00321 !     POSITIVE_DEPTHS). IT SHOULD GIVE HERE THE SAME VALUE
00322 !
00323       IF(NPTFR.GT.0) THEN
00324         DO IPTFR = 1,NPTFR
00325           IF(LIHBOR(IPTFR).EQ.KENT) THEN
00326             FLBOR%R(IPTFR)=0.D0
00327             DO IPLAN = 1,NPLAN-1
00328               I=(IPLAN-1)*NPOIN2+MESH2%NBOR%I(IPTFR)
00329 !             FLUEXT COMPUTED TO SOLVE CONTINUITY IN 3D
00330 !             WITH ASSUMPTION THAT W* IS ZERO.
00331               FLUEXT%R(I)=FLUINT%R(I)+(VOLUN%R(I)-VOLU%R(I))/DT
00332               FLBOR%R(IPTFR)=FLBOR%R(IPTFR)+FLUEXT%R(I)
00333             ENDDO
00334 !           LAST PLANE, RAIN MUST BE TAKEN INTO ACCOUNT
00335             I=(NPLAN-1)*NPOIN2+MESH2%NBOR%I(IPTFR)
00336             IF(RAIN) THEN
00337               FLUEXT%R(I)=FLUINT%R(I)+(VOLUN%R(I)-VOLU%R(I))/DT
00338      &                   +PLUIE%R(MESH2%NBOR%I(IPTFR))
00339             ELSE
00340               FLUEXT%R(I)=FLUINT%R(I)+(VOLUN%R(I)-VOLU%R(I))/DT
00341             ENDIF
00342             FLBOR%R(IPTFR)=FLBOR%R(IPTFR)+FLUEXT%R(I)
00343 !
00344 !           CHECKING THAT SUM OF FLUEXT IS STILL EQUAL TO FLBOR
00345 !           IN THIS CASE DO NOT COMPUTE FLBOR ABOVE
00346 !
00347 !           SUM_FLUEXT=0.D0
00348 !           DO IPLAN = 1,NPLAN
00349 !             I=(IPLAN-1)*NPOIN2+MESH2%NBOR%I(IPTFR)
00350 !             SUM_FLUEXT=SUM_FLUEXT+FLUEXT%R(I)
00351 !           ENDDO
00352 !           IF(ABS(SUM_FLUEXT-FLBOR%R(IPTFR)).GT.1.D-10) THEN
00353 !             PRINT*,'PROBLEM AT POINT ',IPTFR
00354 !             PRINT*,'FLBOR= ',FLBOR%R(IPTFR),' SUM_FLUEXT=',SUM_FLUEXT
00355 !             DO IPLAN = 1,NPLAN
00356 !               I=(IPLAN-1)*NPOIN2+MESH2%NBOR%I(IPTFR)
00357 !               FLUEXT%R(I)=FLUEXT%R(I)*(FLBOR%R(IPTFR)/SUM_FLUEXT)
00358 !             ENDDO
00359 !             CALL PLANTE(1)
00360 !             STOP
00361 !           ENDIF
00362           ENDIF
00363         ENDDO
00364 !
00365 !       SPECIFIC TREATMENT OF POINTS THAT REMAIN WITHOUT VOLUME
00366 !       THIS DOES NOT CHANGE FLBOR
00367 !
00368         IF((OPT_HNEG.EQ.2.OR.SIGMAG).AND.BYPASS) THEN
00369           DO IPTFR = 1,NPTFR
00370             I=MESH2%NBOR%I(IPTFR)
00371             DO IPLAN = 1,NPLAN-1
00372               IF(VOLUN%R(I).LT.1.D-14.AND.VOLU%R(I).LT.1.D-14) THEN
00373 !               FLUEXT GIVEN TO UPPER LAYER
00374                 FLUEXT%R(I+NPOIN2)=FLUEXT%R(I+NPOIN2)+FLUEXT%R(I)
00375                 FLUEXT%R(I)=0.D0
00376               ENDIF
00377               I=I+NPOIN2
00378             ENDDO
00379           ENDDO
00380         ENDIF
00381 !
00382       ENDIF
00383 !
00384 !     ASSEMBLED VERSION OF FLUEXT
00385 !
00386       IF(NCSIZE.GT.1) THEN
00387         CALL OS('X=Y     ',X=FLUEXTPAR,Y=FLUEXT)
00388         CALL PARCOM(FLUEXTPAR,2,MESH3)
00389 !     ELSE
00390 !       FLUEXTPAR%R=>FLUEXT%R   ! DONE ONCE FOR ALL IN POINT_TELEMAC3D
00391       ENDIF
00392 !
00393 !=======================================================================
00394 !
00395 !   COMPUTES THE ADVECTIVE FLUXES THROUGH THE BOTTOM AND FREE SURFACE
00396 !
00397 !=======================================================================
00398 !
00399 !  DIRICHLET TERMS AT THE BOTTOM AND FREE SURFACE
00400 !
00401 !     /
00402 !    /     H * W*  PSII*  D(OMEGA*)
00403 !   /FREE SURFACE AND BOTTOM (IN THE TRANSFPORMED MESH)
00404 !
00405 !   BOTTOM :
00406 !
00407 !     CALL VECTOR
00408 !    &(TRAV2, '=', 'FLUBOR          ', IELM2H, -1.D0, SVIDE, SVIDE,
00409 !    & SVIDE, SVIDE, SVIDE, WSBORF, MESH2,MSK,MASKEL)
00410 !
00411 !     CALL OV ( 'X=X+Y   ' ,FLUEXT%R(1:NPOIN2),
00412 !    &                      TRAV2%R, TRAV2%R, 0.D0, NPOIN2)
00413 !
00414 !   SURFACE :
00415 !
00416 !     CALL VECTOR
00417 !    &(TRAV2, '=', 'FLUBOR          ', IELM2H, 1.D0, SVIDE, SVIDE,
00418 !    & SVIDE, SVIDE, SVIDE, WSBORS, MESH2,MSK,MASKEL)
00419 !
00420 !     CALL OV ( 'X=X+Y   ' ,FLUEXT%R((NPOIN3-NPOIN2+1):NPOIN3),
00421 !    &                      TRAV2%R, TRAV2%R, 0.D0, NPOIN2)
00422 !
00423 !-----------------------------------------------------------------------
00424 !
00425       RETURN
00426       END

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