preadv.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac3d\preadv.f
00002 !
00203                      SUBROUTINE PREADV
00204 !                    *****************
00205 !
00206      &(W,WS,ZPROP,ISOUSI,LT,VOLU,VOLUN)
00207 !
00208 !***********************************************************************
00209 ! TELEMAC3D   V6P2                                   21/08/2010
00210 !***********************************************************************
00211 !
00212 !
00213 !
00214 !
00215 !
00216 !
00217 !
00218 !
00219 !
00220 !
00221 !
00222 !
00223 !
00224 !
00225 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00226 !| ISOUSI         |-->| RANK OF CURRENT SUB-ITERATION
00227 !| LT             |-->| CURRENT TIME STEP NUMBER
00228 !| VOLU           |-->| VOLUMES (INTEGRAL OF BASES) AT NEW TIME STEP
00229 !| VOLUN          |-->| VOLUMES (INTEGRAL OF BASES) AT OLD TIME STEP
00230 !| W              |<->| W VELOCITY
00231 !| WS             |<->| W VELOCITY IN TRANSFORMED MESH
00232 !| ZPROP          |<->| TRANSFORMED ZPROP, TEMPORARILY PUT IN MESH3D%Z
00233 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00234 !
00235       USE BIEF
00236       USE INTERFACE_TELEMAC3D, EX_PREADV => PREADV
00237       USE DECLARATIONS_TELEMAC
00238       USE DECLARATIONS_TELEMAC3D, ONLY : MESH3D,FLUINT,FLUEXT,
00239      &                                   FLUEXTPAR,UCONV,VCONV,T3_01,
00240      &                                   T3_02,T3_03,T3_04,T3_05,
00241      &                                   NETAGE,NPLAN,Z,OPTSUP,FN3D,
00242      &                                   NELEM3,IELM3,IELM2H,IELM2V,
00243      &                                   SVIDE,MSK,MASKEL,MASK_3D,
00244      &                                   LIHBOR,NPTFR2,DT,MESH2D,GRAPRD,
00245      &                                   SIGMAG,T2_01,NPOIN2,NPOIN3,
00246      &                                   DM1,GRAZCO,FLBOR,PLUIE,RAIN,
00247      &                                   FLODEL,FLOPAR,FLULIM,MTRA1,
00248      &                                   N_ADV,BYPASS,WEL,MMURD,MURD_TF,
00249      &                                   WSCONV,VOLU2D,NONHYD,GRADEB,
00250      &                                   ITURBV,MTRA2,FC3D,IT1,IT2,IT3,
00251      &                                   IT4,TRAV3,ZCONV,IKLE2,UCONVC,
00252      &                                   VCONVC,OPT_HNEG,WCONV,ZCHAR,
00253      &                                   NELEM2,MSUPG,UNSV2D,NSCE,
00254      &                                   SOURCES,SEM2D,UNSV3D,U,GRADZF,
00255      &                                   SEM3D,DSSUDT,OPTBAN,INFOGR,
00256      &                                   SLVPRO,AGGLOW,NGAUSS,OPTCHA
00257 !
00258       IMPLICIT NONE
00259       INTEGER LNG,LU
00260       COMMON/INFO/LNG,LU
00261 !
00262 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00263 !
00264       TYPE(BIEF_OBJ), INTENT(INOUT) :: W,WS,ZPROP
00265       TYPE(BIEF_OBJ), INTENT(IN)    :: VOLU,VOLUN
00266 !
00267       INTEGER, INTENT(IN) :: ISOUSI,LT
00268 !
00269 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00270 !
00271       INTEGER ::I,IS,IP,IWS,NSEG3D,OPT_TRID
00272       CHARACTER(LEN=16) FORMUL
00273       CHARACTER(LEN=8) OPER
00274       DOUBLE PRECISION, POINTER, DIMENSION(:) :: SAVEZ
00275 !
00276 !=======================================================================
00277 !
00278 !     MESH MODIFIED TO BE EQUIVALENT TO THE DEPTH USED IN THE 2D
00279 !     CONTINUITY EQUATION, TO CALL : FLUX3D, VECTOR, MATRIX
00280 !
00281 !     ZPROP IS TEMPORARILY PUT IN MESH3D%Z
00282       SAVEZ=>MESH3D%Z%R
00283       MESH3D%Z%R=>ZPROP%R
00284       NSEG3D=MESH3D%NSEG
00285 !
00286 !=======================================================================
00287 !
00288 ! COMPUTES INTERNAL AND EXTERNAL FLUXES AND ADVECTION FIELDS
00289 !
00290 !=======================================================================
00291 !
00292       CALL FLUX3D
00293      & (FLUINT,FLUEXT,FLUEXTPAR,UCONV,VCONV,T3_01,T3_02,T3_03,MESH3D%W,
00294      &  NETAGE,NPLAN,NELEM3,IELM3,IELM2H,IELM2V,SVIDE,MESH3D,
00295      &  MSK,MASKEL,MASK_3D,LIHBOR%I,KENT,NPTFR2,DT,VOLU,VOLUN,MESH2D,
00296      &  GRAPRD,SIGMAG,T2_01,NPOIN2,NPOIN3,DM1,GRAZCO,FLBOR,
00297      &  PLUIE,RAIN,FLODEL,OPT_HNEG,FLULIM,
00298      &  (N_ADV(ADV_LPO).GT.0.OR.N_ADV(ADV_LPO_TF).GT.0),
00299      &  LT,BYPASS,N_ADV,WEL)
00300 !
00301 !=======================================================================
00302 !   COMPUTES (DZW*)JH,IV+1/2 AND ACCUMULATES IN WSCONV
00303 !=======================================================================
00304 !
00305 !     HARDCODED OPTION !!!!!!!!!!!!
00306 !
00307 !     2: DIVERGENCE-FREE FLUXES OBTAINED BY MODIFYING VERTICAL FLUXES
00308 !
00309 !     3: DIVERGENCE-FREE FLUXES OBTAINED BY MODIFYING ALL FLUXES
00310 !        WITH THE HELP OF WCONV AND A PRESSURE EQUATION
00311 !
00312       OPT_TRID=2
00313 !
00314       IF(OPT_TRID.EQ.2) THEN
00315         CALL TRIDW2(WSCONV,VOLU,VOLUN,SEM2D,FLUINT,FLUEXT,SOURCES,
00316      &              NSCE,NETAGE,NPOIN2,DT,UNSV2D,MESH2D)
00317       ELSEIF(OPT_TRID.EQ.3) THEN
00318 !       OTHERWISE WCONV DONE IN WAVE_EQUATION
00319         IF(LT.EQ.0) CALL OS('X=Y     ',X=WCONV,Y=W)
00320         CALL TRIDW3(WSCONV,T3_01,T3_02,T3_03,T3_04,T3_05,MTRA1%D,LT,
00321      &              VOLU,VOLUN,U,UCONV,VCONV,WCONV,DT,NPOIN3,SIGMAG,
00322      &              OPTBAN,MESH3D,MTRA1,MASKEL,NPOIN2,T2_01,NPLAN,
00323      &              FLUEXT,NSCE,SOURCES,RAIN,PLUIE,FLUINT,NETAGE,UNSV2D,
00324      &              SVIDE,NELEM2,MSK,N_ADV,VOLU2D,INFOGR,DSSUDT,IELM3,
00325      &              DM1,GRAZCO,MESH2D,SEM3D,NELEM3,GRADZF)
00326       ENDIF
00327 !
00328 !=======================================================================
00329 !     FOR DEBUGGING: SUMMARY OF ADVECTED VARIABLES AND THEIR SCHEME
00330 !=======================================================================
00331 !
00332 !     DO I=1,15
00333 !       IF(N_ADV(I).GT.0) THEN
00334 !         DO IS=1,N_ADV(I)
00335 !           WRITE(LU,*) 'ADVECTION OF ',
00336 !    &                  BL_FN%ADR(LIST_ADV(IS,I))%P%NAME,
00337 !    &                  ' BY SCHEME ',I
00338 !         ENDDO
00339 !       ENDIF
00340 !     ENDDO
00341 !
00342 !=======================================================================
00343 !     PREPARES ADVECTION BY MURD METHOD
00344 !     STORAGE IS ALWAYS EBE
00345 !=======================================================================
00346 !
00347       IF(N_ADV(ADV_NSC).GT.0.OR.N_ADV(ADV_PSI).GT.0) THEN
00348 !
00349 !       NOTE: THE MATRIX IS THE SAME IN N-SCHEME OR PSI-SCHEME BUT
00350 !             WITH PSI SCHEME THE DIAGONAL IS NOT ASSEMBLED BECAUSE
00351 !             IT IS ASSEMBLED IN MURD3D.
00352 !
00353         IF(N_ADV(ADV_NSC).GT.0) THEN
00354           FORMUL = 'MAMURD 2     N  '
00355         ELSE
00356           FORMUL = 'MAMURD 2     PSI'
00357         ENDIF
00358         CALL MATRIX
00359 !                                                           !!!
00360      &  (MMURD,'M=N     ',FORMUL,IELM3,IELM3,1.D0,DM1,ZCONV,WEL,
00361      &   UCONV,VCONV,WSCONV,MESH3D,MSK,MASKEL)
00362 !       HERE THE BYPASS IS NOT OPTIONAL, OTHERWISE
00363 !       THE SCHEMES ARE NOT MASS-CONSERVATIVE
00364 !       IF(BYPASS) THEN
00365         IF((OPT_HNEG.EQ.2.OR.SIGMAG).AND.IELM3.EQ.41) THEN
00366 !         NOT PROGRAMMED YET FOR TETRAHEDRA !!!!!!!!!!
00367           CALL BYPASS_CRUSHED_POINTS_EBE(VOLU%R,VOLU,VOLUN%R,VOLUN,
00368      &                                   MMURD%X%R,T3_01,MESH2D,MESH3D,
00369      &                                   NPOIN3,NELEM2,NELEM3,NPLAN,
00370      &                                   MESH3D%IKLE%I)
00371 !         OFF-DIAGONAL TERMS MODIFIED, SO DIAGONAL CHANGED ACCORDINGLY
00372           CALL DIAG_MURD(MMURD%D%R,MMURD%X%R,NELEM3,MESH3D%NELMAX,
00373      &                   NPOIN3,MESH3D%IKLE%I,IELM3,
00374      &                   BIEF_DIM2_EXT(IELM3,IELM3,1,'Q',MESH3D))
00375         ENDIF
00376 !       ENDIF
00377 !
00378       ENDIF
00379 !
00380 !=======================================================================
00381 !     PREPARES ADVECTION BY MURD METHOD IN EDGE-BASED FORM
00382 !     STORAGE IS ALWAYS EDGE-BASED
00383 !=======================================================================
00384 !
00385       IF(N_ADV(ADV_NSC_TF).GT.0) THEN
00386 !
00387 !       NOTE: THE MATRIX IS THE SAME IN BOTH CASES BUT
00388 !             WITH PSI SCHEME THE DIAGONAL IS NOT ASSEMBLED
00389 !             IT IS WHAT WE WANT HERE
00390         FORMUL = 'MAMURD 2     PSI'
00391         CALL MATRIX
00392 !                                                             !!!
00393      &  (MURD_TF,'M=N     ',FORMUL,IELM3,IELM3,1.D0,DM1,ZCONV,WEL,
00394      &   UCONV,VCONV,WSCONV,MESH3D,MSK,MASKEL)
00395 !
00396 !       FROM 30 SEGMENTS WITH POSITIVE FLUXES, WE GO TO 15 WITH
00397 !       POSITIVE OR NEGATIVE FLUXES (CASE OF PRISMS, OR 12 TO 6
00398 !       FOR TETRAHEDRA)
00399         DO I=1,NSEG3D
00400           MURD_TF%X%R(I) = MURD_TF%X%R(I) - MURD_TF%X%R(I+NSEG3D)
00401         ENDDO
00402 !       CALL BYPASS: OPTIONAL BUT SAVES ITERATIONS
00403         IF((OPT_HNEG.EQ.2.OR.SIGMAG).AND.BYPASS) THEN
00404           CALL BYPASS_CRUSHED_POINTS_SEG(VOLU%R,VOLU,VOLUN%R,VOLUN,
00405      &                                   MURD_TF%X%R,
00406      &                                   T3_01,MESH2D,MESH3D,
00407      &                                   NPOIN3,ADV_NSC_TF,NPOIN2,
00408      &                                   MESH3D%GLOSEG%I,
00409      &                                   MESH3D%GLOSEG%DIM1,
00410      &                                   MESH2D%NSEG,NPLAN)
00411         ENDIF
00412         IF(NCSIZE.GT.1) THEN
00413 !         ASSEMBLED FORM OF FLUXES STORED IN SECOND PART
00414 !         OF MATRIX WHICH OTHERWISE IS NOT USED
00415           CALL OV('X=Y     ',MURD_TF%X%R(NSEG3D+1:2*NSEG3D),
00416      &                       MURD_TF%X%R(       1:  NSEG3D),
00417      &                       MURD_TF%X%R(       1:  NSEG3D),
00418      &                       0.D0,NSEG3D)
00419           CALL PARCOM2_SEG(MURD_TF%X%R(NSEG3D+1:2*NSEG3D),
00420      &                     MURD_TF%X%R(NSEG3D+1:2*NSEG3D),
00421      &                     MURD_TF%X%R(NSEG3D+1:2*NSEG3D),
00422      &                     MESH2D%NSEG,NPLAN,2,1,MESH2D,2,IELM3)
00423         ENDIF
00424 !
00425       ENDIF
00426 !
00427 !=======================================================================
00428 !     PREPARES LEO POSTMA ADVECTION SCHEMES (PRISMS ONLY)
00429 !=======================================================================
00430 !
00431 !     RETRIEVES VERTICAL FLUXES FROM WSCONV
00432 !     VERTICAL FLUXES ARE STORED IN FLODEL AFTER
00433 !     THE HORIZONTAL FLUXES (THERE ARE NSEG*NPLAN HORIZONTAL FLUXES)
00434 !     USEFUL SIZE OF WSCONV IS (NPOIN2,NPLAN-1)
00435 !
00436 !     THE HORIZONTAL FLUXES HAVE ALREADY BEEN DONE IN FLUX3D
00437 !
00438       IF(N_ADV(ADV_LPO).GT.0.OR.N_ADV(ADV_LPO_TF).GT.0) THEN
00439         IS=MESH2D%NSEG*NPLAN
00440         DO IP=1,NPLAN-1
00441           DO I=1,NPOIN2
00442             IWS=I+(IP-1)*NPOIN2
00443 !           NOTE 1: WSCONV IS ALREADY ASSEMBLED
00444 !                   USING VOLU2D FLODEL WILL BE THE NON ASSEMBLED FORM
00445 !           NOTE 2: WE COULD KEEP THE ORIGINAL RIGHT HAND SIDE IN
00446 !                   TRIDW2
00447 !           NOTE 3: AGAIN CONVENTION REVERSED, HERE FLOW FROM
00448 !                   POINT 2 (UP) TO POINT 1 (DOWN)
00449             FLODEL%R(IS+IWS)=-WSCONV%R(IWS)*VOLU2D%R(I)
00450           ENDDO
00451         ENDDO
00452 !       CALL BYPASS: OPTIONAL BUT SAVES ITERATIONS
00453         IF((OPT_HNEG.EQ.2.OR.SIGMAG).AND.BYPASS) THEN
00454           CALL BYPASS_CRUSHED_POINTS_SEG(VOLU%R,VOLU,VOLUN%R,VOLUN,
00455      &                                   FLODEL%R,
00456      &                                   T3_01,MESH2D,MESH3D,
00457      &                                   NPOIN3,ADV_LPO_TF,NPOIN2,
00458      &                                   MESH3D%GLOSEG%I,
00459      &                                   MESH3D%GLOSEG%DIM1,
00460      &                                   MESH2D%NSEG,NPLAN)
00461         ENDIF
00462 !       FLOPAR = FLODEL ASSEMBLED IN PARALLEL MODE
00463         IF(NCSIZE.GT.1) THEN
00464           CALL OS('X=Y     ',X=FLOPAR,Y=FLODEL)
00465           CALL PARCOM2_SEG(FLOPAR%R,FLOPAR%R,FLOPAR%R,
00466      &                     MESH2D%NSEG,NPLAN,2,1,MESH2D,1,IELM3)
00467         ELSE
00468           FLOPAR%R=>FLODEL%R
00469         ENDIF
00470       ENDIF
00471 !
00472 !=======================================================================
00473 !     PREPARES ADVECTION BY SUPG METHOD
00474 !=======================================================================
00475 !
00476       IF(N_ADV(ADV_SUP).GT.0) THEN
00477 !
00478         IF(OPTSUP(1).EQ.2) THEN
00479 !         HORIZONTAL UPWIND (HERE UPWIND COEFFICIENT=CFL)
00480           FORMUL = 'MAUGUG2         '
00481           CALL MATRIX
00482      &    (MSUPG,'M=N     ',FORMUL,IELM3,IELM3,0.5D0*DT,SVIDE,SVIDE,
00483      &     SVIDE,UCONV,VCONV,WSCONV,MESH3D,MSK,MASKEL)
00484 !         MSUPG IS SYMMETRICAL
00485         ELSEIF(OPTSUP(1).EQ.1) THEN
00486 !         HORIZONTAL UPWIND (HERE UPWIND COEFFICIENT=1)
00487           FORMUL = 'MAUGUG1         '
00488           CALL MATRIX
00489      &    (MSUPG,'M=N     ',FORMUL,IELM3,IELM3,1.D0,SVIDE,SVIDE,
00490      &     SVIDE,UCONV,VCONV,WSCONV,MESH3D,MSK,MASKEL)
00491 !         MSUPG IS NOT SYMMETRICAL
00492         ELSEIF(OPTSUP(1).NE.0) THEN
00493           IF (LNG.EQ.1) WRITE(LU,*) 
00494 'VALEUR NON PREVU DE OPTSUP     &                               AND PREADV'
00495           IF (LNG.EQ.2) WRITE(LU,*) 
00496 'UNEXPECTED VALUE OF OPTSUP     &                               ET PREADV'
00497           CALL PLANTE(1)
00498           STOP
00499         ENDIF
00500 !
00501 !       MSUPG TRANSFORMED INTO NON SYMMETRICAL MATRIX
00502         IF(OPTSUP(1).EQ.2) THEN
00503           CALL OM('M=X(M)  ',MSUPG,MSUPG,SVIDE,0.D0,MESH3D)
00504           OPER='M=M+N   '
00505         ELSEIF(OPTSUP(1).EQ.1) THEN
00506           OPER='M=M+N   '
00507         ELSE
00508           OPER='M=N     '
00509         ENDIF
00510 !
00511 !       ADDS CENTRED ADVECTION TERM
00512 !
00513         FORMUL = 'MATVGR          '
00514         FORMUL(8:8) = '2'
00515         CALL MATRIX
00516      &  (MSUPG,OPER,FORMUL,IELM3,IELM3,1.D0,DM1,ZCONV,SVIDE,
00517      &   UCONV,VCONV,WSCONV,MESH3D,MSK,MASKEL)
00518 !
00519 !       VERTICAL UPWIND (SUBROUTINE UPWIND EXPECTS SYMMETRICAL MATRICES)
00520 !       HERE UPWIND COEFFICIENT = 1, BUT WSCONV USED INSTEAD OF W.
00521 !       WITH TETRAHEDRONS UPWINDING IS DONE IN 3D, SO NO NEED
00522 !       TO CALL THIS UPWIND ON THE VERTICAL.
00523 !
00524         IF(IELM3.EQ.41) THEN
00525           CALL UPWIND(MSUPG,WSCONV,1.D0,MESH2D,MESH3D,NPLAN)
00526         ENDIF
00527 !
00528       ENDIF
00529 !
00530 !=======================================================================
00531 !
00532 !     RESTORES MESH3D%Z
00533 !
00534       MESH3D%Z%R=>SAVEZ
00535 !
00536 !=======================================================================
00537 !
00538 !     COMPUTES DELTAZ*WSTAR (IN WS) AT NODES
00539 !
00540       CALL WSTAR(WS,WSCONV,Z,NPOIN2,NPLAN)
00541 !
00542 !=======================================================================
00543 !
00544 !     COMPUTES W FROM  (DZW*)JH,IV+1/2
00545 !
00546 !        (WITH HYDROSTATIC ASSUMPTION, W IS NEVER USED,
00547 !                  IT IS DONE HERE FOR OUTPUTS)
00548 !        HOWEVER IT IS ALWAYS USED WITH THE K-EPSILON OR K-OMEGA MODELS
00549 !
00550       IF(.NOT.NONHYD) THEN
00551         IF(((LT/GRAPRD)*GRAPRD.EQ.LT.AND.LT.GE.GRADEB).OR.
00552      &      (ITURBV.EQ.3.OR.ITURBV.EQ.7)) THEN
00553           CALL WSTARW(W,WSCONV,T3_03%R,T3_04%R,T3_05%R)
00554         ENDIF
00555       ENDIF
00556 !
00557 !=======================================================================
00558 !
00559 ! ADVECTION BY METHOD OF CHARACTERISTICS
00560 !
00561 !=======================================================================
00562 !
00563       IF(N_ADV(ADV_CAR).GT.0) THEN
00564 !
00565 !       NOTES:
00566 !
00567 !       IN BLOCK FN3D THERE IS U,V,W INSTEAD OF UN,VN,WN
00568 !       BECAUSE ADVECTION IS DONE FOR THE NEXT TIME STEP
00569 !
00570 !       FN3D IS THE BLOCK OF VARIABLES ADVECTED WITH
00571 !       SCHEME ADV_CAR (SEE POINT_TELEMAC3D)
00572 !
00573 !       WITH MERCATOR PROJECTION, UCONVC AND VCONVC HAVE BEEN DIVIDED
00574 !       BY THE COSINUS OF LATITUDE, BECAUSE WE WORK HERE ON X AND Y
00575 !       WHICH ARE DISTORTED BY THE PROJECTION. UNLIKE IN 2D
00576 !       UCONVC AND VCONVC ARE NOT MULTIPLIED BY THE COSINUS AFTER
00577 !       USE BECAUSE THE SINGLE USE IS THERE.
00578 !
00579         CALL CHARAC(FN3D,FC3D,FC3D%N,UCONVC,VCONVC,WS,WS,ZCHAR,ZCHAR,
00580      &              DT,MESH3D%IFABOR,IELM3,NPOIN2,NPLAN,1,1,
00581      &              MSK,MASKEL,MTRA2%X,MTRA2%D,MTRA2%D,TRAV3,
00582      &              IT1%I,IT2%I,IT2%I,IT3%I,IT4%I,IT2%I,
00583      &              MESH3D,NELEM2,MESH2D%NELMAX,IKLE2,MESH2D%SURDET,
00584      &              MTRA1,SEM3D,SLVPRO,AGGLOW,INFOGR,NGAUSS,UNSV3D,
00585      &              OPTCHA,SIGMA=.TRUE.)
00586 !
00587       ENDIF
00588 !
00589 !-----------------------------------------------------------------------
00590 !
00591       RETURN
00592       END

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