wave_equation.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac3d\wave_equation.f
00002 !
00248                      SUBROUTINE WAVE_EQUATION
00249 !                    ************************
00250 !
00251      &(LT,ISOUSI)
00252 !
00253 !***********************************************************************
00254 ! TELEMAC3D   V7P0                                   07/07/2014
00255 !***********************************************************************
00256 !
00257 !
00258 !
00259 !
00260 !
00261 !
00262 !
00263 !
00264 !
00265 !
00266 !
00267 !
00268 !
00269 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00270 !| ISOUSI         |-->| RANK OF CURRENT SUB-ITERATION
00271 !| LT             |-->| CURRENT TIME STEP NUMBER
00272 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00273 !
00274       USE BIEF
00275       USE INTERFACE_TELEMAC3D, EX_WAVE_EQUATION => WAVE_EQUATION
00276       USE DECLARATIONS_TELEMAC
00277       USE DECLARATIONS_TELEMAC3D
00278 !
00279       IMPLICIT NONE
00280       INTEGER LNG,LU
00281       COMMON/INFO/LNG,LU
00282 !
00283 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00284 !
00285       INTEGER, INTENT(IN) :: LT,ISOUSI
00286 !
00287 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00288 !
00289       CHARACTER(LEN=16) :: FORMUL
00290       INTEGER           :: I,IPTFR3,IPOIN2,IPLAN,I1,I2,I3,NP,IELEM
00291       INTEGER           :: I3D,IP,IAD1,IAD2,IAD3
00292       LOGICAL           :: MSKGRA
00293       DOUBLE PRECISION  :: C
00294 !
00295 !-----------------------------------------------------------------------
00296 !
00297 !     DEFINES POINTERS TO RENAME PARTS OF MEMORY
00298 !     MTRA2%X%R HAS AT LEAST THE SIZE 30*NELMAX (STORAGE 1)
00299 !                                  OR  2*NSEG3D (STORAGE 3)
00300 !
00301 !     WE NEED: 8*NPOIN
00302 !           OR 8*NPOIN
00303 !
00304 !
00305       DOUBLE PRECISION, POINTER :: TRIC(:),TRID(:),TRIE(:)
00306       DOUBLE PRECISION, POINTER :: UAUX(:),VAUX(:)
00307       TRIC=>MTRA2%X%R(          1:   NPOIN3)
00308       TRID=>MTRA2%X%R(   NPOIN3+1: 2*NPOIN3)
00309       TRIE=>MTRA2%X%R( 2*NPOIN3+1: 3*NPOIN3)
00310 !         =>MTRA2%X%R( 3*NPOIN3+1: 4*NPOIN3) USED IN TRID3D
00311 !         =>MTRA2%X%R( 4*NPOIN3+1: 5*NPOIN3) USED IN TRID3D
00312       UAUX=>MTRA2%X%R( 5*NPOIN3+1: 6*NPOIN3)
00313       VAUX=>MTRA2%X%R( 6*NPOIN3+1: 7*NPOIN3)
00314 !
00315       IF(7*NPOIN3.GT.30*MESH3D%NELMAX.OR.
00316      &   7*NPOIN3.GT.2*MESH3D%NSEG) THEN
00317         WRITE(LU,*) 'PROBLEME DE PLACE MEMOIRE DANS WAVE_EQUATION'
00318         WRITE(LU,*) 'NPOIN3=',NPOIN3
00319         WRITE(LU,*) 'NELMAX=',MESH3D%NELMAX
00320         WRITE(LU,*) 'NSEG=',MESH3D%NSEG
00321         CALL PLANTE(1)
00322         STOP
00323       ENDIF
00324 !
00325 !=======================================================================
00326 !     COMPUTES THE DIFFUSION TERMS DIFF
00327 !
00328 !       AND THEN UC + DT(F - DIFF -G GRAD(Z))
00329 !
00330 !       STORED IN T3_01 AND T3_02
00331 !
00332 !=======================================================================
00333 !
00334       FORMUL='MATDIF          '
00335       IF(INCHYD) FORMUL(7:7)='2'
00336 !     NOTE: COULD BE OPTIMISED IF SCHDVI=0
00337       CALL MATRIX(MDIFF,'M=N     ',FORMUL,IELM3,IELM3,1.D0,
00338      &            VISCVI%ADR(1)%P,VISCVI%ADR(2)%P,VISCVI%ADR(3)%P,
00339      &            SVIDE,SVIDE,SVIDE,MESH3D,MSK,MASKEL)
00340 !
00341 !     IMPLICITATION OF DIAGONAL TERMS
00342 !     BUILDS A TRIDIAGONAL MATRIX IN OFF-DIAGONAL TERMS OF MTRA2
00343 !
00344       CALL GETTRI(MTRA2%X%R,MDIFF,TETADI,MESH3D,NPLAN,MESH2D%NPOIN,
00345      &            MESH2D%NSEG,IELM3,NELEM2)
00346 !
00347       DO I=1,U%DIM1
00348         TRIC(I)=TRIC(I)*UNSV3D%R(I)*DT
00349         TRID(I)=TRID(I)*UNSV3D%R(I)*DT
00350         TRIE(I)=TRIE(I)*UNSV3D%R(I)*DT
00351       ENDDO
00352 !
00353 !     EXPLICIT DIFFUSION TERMS
00354 !
00355       CALL MATVEC ('X=AY     ',T3_01,MDIFF,UN,0.D0,MESH3D)
00356       CALL MATVEC ('X=AY     ',T3_02,MDIFF,VN,0.D0,MESH3D)
00357 !
00358 !     EXPLICIT STRESS TERMS
00359 !
00360       CALL T3D_STRESS(T3_01,'X=X-Y   ',T2_02,T3_04,
00361      &                BUBORL,BUBORF,BUBORS,NPOIN2,NPOIN3,MESH2D,
00362      &                MESH3D,IELM3,IELM2H,IELM2V,SVIDE,
00363      &                MSK,MASKBR,MASKEL,IPBOT%I,SIGMAG,OPTBAN,NPLAN)
00364       CALL T3D_STRESS(T3_02,'X=X-Y   ',T2_02,T3_04,
00365      &                BVBORL,BVBORF,BVBORS,NPOIN2,NPOIN3,MESH2D,
00366      &                MESH3D,IELM3,IELM2H,IELM2V,SVIDE,
00367      &                MSK,MASKBR,MASKEL,IPBOT%I,SIGMAG,OPTBAN,NPLAN)
00368 !
00369 !     REQUIRES REAL VALUES IN PARALLEL MODE
00370 !
00371       IF(NCSIZE.GT.1) THEN
00372         CALL PARCOM(T3_01,2,MESH3D)
00373         CALL PARCOM(T3_02,2,MESH3D)
00374       ENDIF
00375 !
00376 !     DIVIDES BY VOLUME OF BASES
00377 !
00378       CALL OS('X=XY    ',X=T3_01,Y=UNSV3D)
00379       CALL OS('X=XY    ',X=T3_02,Y=UNSV3D)
00380 !
00381 !     COMPUTES UC + DT(F - DIFF -G GRAD(Z))
00382 !     STARTS THE COMPUTATION OF UAUX AND VAUX
00383 !
00384 !     NEW IN VERSION 5.8 : SUPG SCHEME IS POSSIBLE
00385 !                          AND RESULT OF ADVECTION IS THEN IN UD AND VD
00386 !                          AFTER CALL TO CVDF3D (BUT DIFFUSION AND
00387 !                          SOURCES TERMS NOT DONE IN CVDF3D, SEE
00388 !                          SCHDVI_HOR AND YAS0U, YAS1U IN TELEMAC3D.F)
00389 !
00390       IF(SCHCVI.EQ.ADV_SUP) THEN
00391         CALL OS('X=Y     ',X=UC,Y=UD)
00392         CALL OS('X=Y     ',X=VC,Y=VD)
00393       ENDIF
00394 !
00395       IF(  (SCHCVI.NE.ADV_NSC.AND.
00396      &      SCHCVI.NE.ADV_PSI.AND.
00397      &      SCHCVI.NE.ADV_LPO.AND.
00398      &      SCHCVI.NE.ADV_NSC_TF.AND.
00399      &      SCHCVI.NE.ADV_LPO_TF).AND.S0U%TYPR.NE.'0') THEN
00400 !       CASES WHERE S0U MUST BE TAKEN INTO ACCOUNT:
00401 !       IT IS NOT 0 AND IT HAS NOT BEEN TREATED BY THE ADVECTION SCHEME
00402 !       (SO FAR ACTUALLY CASES ADV_SUPG AND ADV_CAR)
00403         DO I=1,U%DIM1
00404           I2=MOD(I-1,NPOIN2)+1
00405           T3_01%R(I)=UC%R(I)
00406      &    +DT*(S0U%R(I)-T3_01%R(I)-TETAZCOMP*GRAV*GRADZN%ADR(1)%P%R(I2))
00407           T3_02%R(I)=VC%R(I)
00408      &    +DT*(S0V%R(I)-T3_02%R(I)-TETAZCOMP*GRAV*GRADZN%ADR(2)%P%R(I2))
00409         ENDDO
00410       ELSE
00411 !       FORMULA WITHOUT S0U
00412         DO I=1,U%DIM1
00413           I2=MOD(I-1,NPOIN2)+1
00414           T3_01%R(I)=UC%R(I)
00415      &    +DT*(-T3_01%R(I)-TETAZCOMP*GRAV*GRADZN%ADR(1)%P%R(I2))
00416           T3_02%R(I)=VC%R(I)
00417      &    +DT*(-T3_02%R(I)-TETAZCOMP*GRAV*GRADZN%ADR(2)%P%R(I2))
00418         ENDDO
00419       ENDIF
00420 !
00421 !=======================================================================
00422 !     COMPUTES 1/(1+DT*(FROT3D+S1U))      IN T3_04
00423 !
00424 !       IT IS ASSUMED THAT S1V=S1U         (IT SHOULD FOR TENSORIALITY)
00425 !                     THAT AVBORL=AUBORL
00426 !                     THAT AVBORF=AUBORF
00427 !                     THAT AVBORS=AUBORS
00428 !
00429 !=======================================================================
00430 !
00431 !     COMPUTES THE FRICTION TERM FROT3D + IMPLICIT SOURCES
00432 !
00433       CALL CPSTVC(S1U,T3_04)
00434 !     ERASES ALL VALUES (ONLY BOUNDARY VALUES WILL BE CHANGED BELOW)
00435       CALL OS('X=0     ',X=T3_04)
00436 !
00437 !     BOTTOM (MASS-LUMPED FORM AS IN 2D):
00438 !
00439       IF(AUBORF%TYPR.NE.'0') THEN
00440 !       VERSION WITH 1/COS(SLOPE) (OTHERWISE USE VOLU2D INSTEAD OF T2_01)
00441         CALL SLOPES(TE1,ZF,MESH2D)
00442         CALL VECTOR(T2_01,'=','MASBAS          ',IELM2H,1.D0,SVIDE,
00443      &              SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH2D,.TRUE.,TE1)
00444 !       IF(NCSIZE.GT.1) DONE ON T3_04 AT THE END
00445 !
00446         IF(SIGMAG.OR.OPTBAN.EQ.1) THEN
00447 !         TREATS CRUSHED LAYERS AND TIDAL FLATS IN THE SAME WAY
00448           DO IPOIN2=1,NPOIN2
00449             DO NP=0,IPBOT%I(IPOIN2)
00450               I=NP*NPOIN2+IPOIN2
00451               T3_04%R(I)=-AUBORF%R(IPOIN2)*T2_01%R(IPOIN2)
00452             ENDDO
00453           ENDDO
00454         ELSE
00455           DO I=1,NPOIN2
00456             T3_04%R(I)=-AUBORF%R(I)*T2_01%R(I)
00457           ENDDO
00458         ENDIF
00459       ENDIF
00460 !
00461 !     LATERAL FACES (MASS-LUMPED FORM)
00462 !
00463       IF(AUBORL%TYPR.NE.'0') THEN
00464         CALL VECTOR(T3_06, '=','MASBAS          ',IELM2V,+1.D0,SVIDE,
00465      &              SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH3D,MSK, MASKEL)
00466 !       IF(NCSIZE.GT.1) : DONE ON THE FINAL T3_04
00467         CALL OSDB( 'X=X-YZ  ' ,T3_04,AUBORL,T3_06,C,MESH3D)
00468       ENDIF
00469 !
00470 !     SURFACE (MASS-LUMPED FORM):
00471 !
00472       IF(AUBORS%TYPR.NE.'0') THEN
00473         CALL OV('X=X-YZ  ',T3_04%R(NPOIN3-NPOIN2+1:NPOIN3),
00474      &                     AUBORS%R,VOLU2D%R,C,NPOIN2)
00475       ENDIF
00476 !
00477       IF(NCSIZE.GT.1) CALL PARCOM(T3_04,2,MESH3D)
00478 !
00479 !     COMPUTES THE INVERSE OF THE DENOMINATOR FOR U(N+1) AND V(N+1)
00480 !
00481       IF(S1U%TYPR.NE.'0') THEN
00482         DO I=1,U%DIM1
00483           TRID(I)=TRID(I)+1.D0+DT*(S1U%R(I)+T3_04%R(I)*UNSV3D%R(I))
00484         ENDDO
00485       ELSE
00486         DO I=1,U%DIM1
00487           TRID(I)=TRID(I)+1.D0+DT*(T3_04%R(I)*UNSV3D%R(I))
00488         ENDDO
00489       ENDIF
00490 !
00491 !     COMPUTES THE SOLUTION OF TRI * X = UNITY VECTOR EVERYWHERE
00492 !     PUT IN INV1
00493 !
00494       CALL OS('X=C     ',X=T3_04,C=1.D0)
00495       CALL TRID3D(MTRA2%X%R,DM1%R,T3_04%R,NPOIN3,NPOIN2)
00496 !
00497 !     LATERAL BOUNDARY CONDITION: CANCELS DM1 FOR THE VELOCITY DIRICHLET
00498 !
00499       DO IPTFR3 = 1,NPTFR3
00500         IF(LIUBOL%I(IPTFR3).EQ.KENT.OR.
00501      &     LIUBOL%I(IPTFR3).EQ.KENTU.OR.LIUBOL%I(IPTFR3).EQ.KADH) THEN
00502           DM1%R(MESH3D%NBOR%I(IPTFR3)) = 0.D0
00503         ENDIF
00504       ENDDO
00505 !
00506 !=======================================================================
00507 !     COMPUTES THE NEW DEPTH WITH WAVE EQUATION
00508 !=======================================================================
00509 !
00510 !     STARTS COMPUTATION OF THE SECOND MEMBER (IN SEM2D%ADR(1)%P)
00511 !
00512       CALL OS('X=Y     ',X=SEM2D%ADR(1)%P,Y=SMH)
00513 !
00514 !     PSEUDO-VISCOSITY IN THE WAVE EQUATION (IN NUWAVE, P0 FUNCTION)
00515 !
00516       CALL NUWAVE_P0(NUWAVE%R,DM1%R,Z,T3_03%R,MESH3D%IKLE%I,
00517      &               NPOIN2,NPLAN,MESH3D%NELEM,MESH3D%NELMAX,NELEM2,
00518      &               GRAV*TETAH*TETAU*DT,IELM3,MESH3D%X%R,MESH3D%Y%R,
00519      &               MESH2D%SURFAC%R)
00520 !
00521 !     CORRESPONDING DIFFUSION MATRIX
00522 !
00523       CALL MATRIX(MAT2D%ADR(1)%P,'M=N     ','MATDIF          ',
00524      &            IELM2H,IELM2H,1.D0,SVIDE,SVIDE,SVIDE,
00525      &            NUWAVE,SVIDE,SVIDE,MESH2D,MSK,MASKEL)
00526 !
00527 !     STORES THIS MATRIX FOR THE COMPUTATION OF FLINT2
00528 !
00529       CALL OM('M=N     ',MAT2D%ADR(2)%P,MAT2D%ADR(1)%P,
00530      &        SVIDE,0.D0,MESH2D)
00531 !
00532 !     SEM2D%ADR(1)%P = SEM2D%ADR(1)%P + INTEGRAL ON OMEGA3D
00533 !
00534 !     3D VECTOR TO INTEGRATE (IN UCONV, VCONV)
00535 !
00536 !     COMPUTES UAUX BY SOLVING TRIDIAGONAL SYSTEMS
00537 !
00538       CALL TRID3D(MTRA2%X%R,UAUX,T3_01%R,NPOIN3,NPOIN2)
00539       CALL TRID3D(MTRA2%X%R,VAUX,T3_02%R,NPOIN3,NPOIN2)
00540 !
00541 !     TAKES THE PRESSURE GRADIENT INTO ACCOUNT
00542 !
00543       IF(NONHYD.AND.DPWAVEQ) THEN
00544 !
00545 !       COMPUTES AN ESTIMATE OF THE DYNAMIC PRESSURE WITH UAUX TAKEN
00546 !       AS U(N+1). THIS ESTIMATE WILL BE THE RESULT GIVEN IN THE RESULT
00547 !       FILE, AS DP IN THE SECOND CALL TO PREDIV IS (ONLY) INCREMENTAL.
00548 !
00549         CALL CPSTVC(UN,T3_04)
00550         CALL CPSTVC(VN,T3_05)
00551         CALL CPSTVC(WD,T3_06)
00552         CALL OV('X=Y     ',T3_04%R,UAUX,UAUX,0.D0,NPOIN3)
00553         CALL OV('X=Y     ',T3_05%R,VAUX,VAUX,0.D0,NPOIN3)
00554         CALL OS('X=Y     ',X=T3_06,Y=WD)
00555 !       NON COMPATIBLE PART OF FREE SURFACE GRADIENT
00556 !      (WHICH IS NOT IN UAUX, SEE ALSO FINAL COMPUTATION OF U AND V)
00557         IF(ABS(1.D0-TETAZCOMP).GT.1.D-6) THEN
00558           C=-DT*GRAV*(1.D0-TETAZCOMP)
00559           DO I=1,U%DIM1
00560             T3_04%R(I)=T3_04%R(I)+
00561      &      C*GRADZN%ADR(1)%P%R(MOD(I-1,NPOIN2)+1)*DM1%R(I)
00562             T3_05%R(I)=T3_05%R(I)+
00563      &      C*GRADZN%ADR(2)%P%R(MOD(I-1,NPOIN2)+1)*DM1%R(I)
00564           ENDDO
00565         ENDIF
00566 !
00567 !       TAKES DH INTO ACCOUNT, IF KNOWN (I.E. FROM 2ND SUBITERATIONS ON)
00568 !       TAKING DH FROM PREVIOUS TIMESTEP IS NOT A GOOD IDEA
00569 !
00570         IF(ISOUSI.GT.1) THEN
00571           CALL VECTOR(T2_02,'=','GRADF          X',
00572      &                IELM2H,-GRAV*TETAH,DH,
00573      &                SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,
00574      &                MESH2D,MSK,MASKEL)
00575           CALL VECTOR(T2_03,'=','GRADF          Y',
00576      &                IELM2H,-GRAV*TETAH,DH,
00577      &                SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,
00578      &                MESH2D,MSK,MASKEL)
00579           IF(NCSIZE.GT.1) THEN
00580             CALL PARCOM(T2_02,2,MESH2D)
00581             CALL PARCOM(T2_03,2,MESH2D)
00582           ENDIF
00583           CALL OS('X=XY    ',X=T2_03,Y=UNSV2D)
00584           CALL OS('X=XY    ',X=T2_02,Y=UNSV2D)
00585           DO I=1,U%DIM1
00586             T3_04%R(I)=T3_04%R(I)+DT*T2_02%R(MOD(I-1,NPOIN2)+1)*DM1%R(I)
00587             T3_05%R(I)=T3_05%R(I)+DT*T2_03%R(MOD(I-1,NPOIN2)+1)*DM1%R(I)
00588           ENDDO
00589         ENDIF
00590 !
00591 !       APPLIES LATERAL BOUNDARY CONDITIONS
00592 !       NOT VERY SIGNIFICANT
00593 !
00594 !       DO IPTFR3 = 1,NPTFR3
00595 !         IF(LIUBOL%I(IPTFR3).EQ.KENT .OR.
00596 !    *       LIUBOL%I(IPTFR3).EQ.KENTU.OR.
00597 !    *       LIUBOL%I(IPTFR3).EQ.KADH) THEN
00598 !            T3_04%R(MESH3D%NBOR%I(IPTFR3)) = UBORL%R(IPTFR3)
00599 !         ENDIF
00600 !         IF(LIVBOL%I(IPTFR3).EQ.KENT .OR.
00601 !    *       LIVBOL%I(IPTFR3).EQ.KENTU.OR.
00602 !    *       LIVBOL%I(IPTFR3).EQ.KADH) THEN
00603 !            T3_05%R(MESH3D%NBOR%I(IPTFR3)) = VBORL%R(IPTFR3)
00604 !         ENDIF
00605 !       ENDDO
00606 !
00607 !       BEWARE: PREDIV WILL ERASE ALL T3_** WORK ARRAYS BECAUSE CALLS SOLVE
00608         CALL PREDIV(DP,T3_04,T3_05,T3_06,INFOGR,.TRUE.,1,
00609      &              .TRUE.,.TRUE.,.TRUE.)
00610 !       APPLIES CORRECTION TO UAUX
00611         CALL VELRES(UAUX,VAUX,WD%R,DP,
00612      &              T3_08,T3_09,T3_10,MSK,MASKEL,MESH3D,
00613      &              SVIDE,IELM3,NPLAN,OPTBAN,UNSV3D,.FALSE.,
00614      &              NPOIN3,NPOIN2,
00615      &              SIGMAG,IPBOT%I,AGGLOH)
00616 !
00617       ENDIF
00618 !
00619 !     END OF 'TAKES THE PRESSURE GRADIENT INTO ACCOUNT'
00620 !
00621 !
00622 !     LATERAL BOUNDARY CONDITION IN UAUX AND VAUX
00623 !
00624       DO IPTFR3 = 1,NPTFR3
00625         IF(LIUBOL%I(IPTFR3).EQ.KENT.OR.
00626      &     LIUBOL%I(IPTFR3).EQ.KENTU.OR.LIUBOL%I(IPTFR3).EQ.KADH) THEN
00627           UAUX(MESH3D%NBOR%I(IPTFR3)) = UBORL%R(IPTFR3)
00628         ENDIF
00629         IF(LIVBOL%I(IPTFR3).EQ.KENT.OR.
00630      &     LIVBOL%I(IPTFR3).EQ.KENTU.OR.LIVBOL%I(IPTFR3).EQ.KADH) THEN
00631           VAUX(MESH3D%NBOR%I(IPTFR3)) = VBORL%R(IPTFR3)
00632         ENDIF
00633       ENDDO
00634 !
00635 !     COMPUTES TETAU * UAUX + (1-TETAU) * UN
00636 !
00637       DO I=1,U%DIM1
00638         UCONV%R(I)=TETAU*UAUX(I)+(1.D0-TETAU)*UN%R(I)
00639         VCONV%R(I)=TETAU*VAUX(I)+(1.D0-TETAU)*VN%R(I)
00640       ENDDO
00641 !
00642 !     SEM2D%ADR(1)%P = SEM2D%ADR(1)%P - FLUX2D
00643 !
00644 !     UNONNEU=8 : 1 IF NOT A WALL
00645 !     CALL EXTMSK(MASKBR,MASK%ADR(8)%P%R,MESH2D%NPTFR,NETAGE)
00646 !     CALL VECTOR(T3_06,'=','FLUBOR          ',IELBOR(IELM3,2),
00647 !    &            1.D0,SVIDE,SVIDE,SVIDE,UCONV,VCONV,SVIDE,
00648 !    &            MESH3D,.TRUE.,MASKBR)
00649       CALL VECTOR(T3_06,'=','FLUBOR          ',IELBOR(IELM3,2),
00650      &            1.D0,SVIDE,SVIDE,SVIDE,UCONV,VCONV,SVIDE,
00651      &            MESH3D,.TRUE.,MASK_3D%ADR(8)%P)
00652 !
00653       CALL SUMVER(FLBOR%R,T3_06%R,NPOIN2,NPLAN,MESH2D%NPTFR)
00654 !
00655       CALL OSDB( 'X=X-Y   ',SEM2D%ADR(1)%P,FLBOR,FLBOR,C,MESH2D)
00656 !
00657 !     MULTIPLIES BY THE GRADIENT OF THE 3D BASES
00658 !
00659       FORMUL = 'VGRADP       HOR'
00660       CALL VECTOR(T3_01,'=',FORMUL,IELM3,1.D0,SVIDE,SVIDE,SVIDE,
00661      &            UCONV,VCONV,SVIDE,MESH3D,MSK,MASKEL)
00662 !
00663 !     SUM ON THE VERTICAL
00664 !     FLINT2 WILL BE ADDED TO SEM2D, BUT MAY BE USED TO CHECK IN TRIDW2
00665 !     THAT THE SUM ON THE VERTICAL OF FLUINT = FLINT2
00666 !
00667       CALL OS('X=0     ',X=FLINT2)
00668       DO IPLAN=1,NPLAN
00669         DO I=1,NPOIN2
00670           FLINT2%R(I)=FLINT2%R(I)+T3_01%R(I+(IPLAN-1)*NPOIN2)
00671         ENDDO
00672       ENDDO
00673 !
00674 !=======================================================================
00675 !     CONTRIBUTION OF NON COMPATIBLE LAPLACIAN
00676 !     SEE ALSO MODIFICATION OF ZCONV LATER ON
00677 !=======================================================================
00678 !
00679       IF(ABS(1.D0-TETAZCOMP).GT.1.D-6) THEN
00680 !
00681 !       ADDS NON COMPATIBLE LAPLACIAN
00682 !       FOR THE CONTRIBUTION OF EXPLICIT FREE-SURFACE
00683 !       TO THE VELOCITY
00684 !       TETAZCOMP=1 : COMPATIBLE
00685 !       TETAZCOMP=0 : NON COMPATIBLE
00686 !
00687 !       ADDS THE NON COMPATIBLE EXPLICIT FREE SURFACE GRADIENT
00688         CALL CPSTVC(H,T2_04)
00689         CALL CPSTVC(H,T2_05)
00690         IF(OPTBAN.EQ.1) THEN
00691 !         FREE SURFACE PIECE-WISE LINEAR IN ZFLATS
00692           CALL VECTOR(FLINT2,'+','VGRADP 2        ',IELM2H,
00693      &                -(1.D0-TETAZCOMP)/TETAH,
00694      &                SVIDE,SVIDE,SVIDE,NUWAVE,ZFLATS,SVIDE,
00695      &                MESH2D,MSK,MASKEL)
00696         ELSE
00697 !         FREE SURFACE LINEAR IN T2_04
00698           DO I=1,NPOIN2
00699             T2_04%R(I)=HN%R(I)+ZF%R(I)
00700           ENDDO
00701           IF(ATMOS) THEN
00702 !         ADDING ATMOSPHERIC PRESSURE
00703             C=1.D0/(RHO0*GRAV)
00704             DO I=1,NPOIN2
00705               T2_04%R(I)=T2_04%R(I)+C*PATMOS%R(I)
00706             ENDDO
00707           ENDIF
00708           CALL VECTOR(FLINT2,'+','VGRADP 2        ',IELM2H,
00709      &                -(1.D0-TETAZCOMP)/TETAH,
00710      &                SVIDE,SVIDE,SVIDE,NUWAVE,T2_04,SVIDE,
00711      &                MESH2D,MSK,MASKEL)
00712         ENDIF
00713 !
00714       ENDIF
00715 !
00716       CALL OS('X=X+Y   ',X=SEM2D%ADR(1)%P,Y=FLINT2)
00717 !
00718 !=======================================================================
00719 !     ADDS THE MASS MATRIX (LUMPED OR NOT) TO THE SYSTEM MATRIX
00720 !     SOLVES THE EQUATION
00721 !=======================================================================
00722 !
00723 !     COMPUTES THE PARTIALLY LUMPED 2D MASS MATRIX
00724 !
00725       FORMUL='MATMAS          '
00726 !     NOTE: THERE IS LOCAL LUMPING IN PROPAG
00727 !           ON THE TIDAL FLATS
00728       CALL MATRIX(MAT2D%ADR(3)%P,'M=N     ',FORMUL,IELM2H,IELM2H,
00729      &            1.D0/DT,
00730      &            SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,
00731      &            MESH2D,MSK,MASKEL)
00732       CALL LUMP(T2_01,MAT2D%ADR(3)%P,MESH2D,AGGLOH)
00733       CALL OM('M=M+CN  ',MAT2D%ADR(1)%P,MAT2D%ADR(3)%P,
00734      &        SVIDE,1.D0-AGGLOH,MESH2D)
00735       CALL OM('M=M+D   ',MAT2D%ADR(1)%P,MAT2D%ADR(1)%P,
00736      &        T2_01,C,MESH2D)
00737 !
00738 !     INITIAL GUESS FOR DH
00739 !
00740       IF(IORDRH.EQ.0) THEN
00741         CALL OS('X=0     ',X=DH)
00742       ELSEIF(IORDRH.EQ.1) THEN
00743 !       DONE IN TELEMAC3D.F
00744 !       IF(LT.EQ.1.AND.ISOUSI.EQ.1) CALL OS('X=0     ',X=DH)
00745       ELSE
00746         IF(LNG.EQ.1) WRITE(LU,30) IORDRH
00747         IF(LNG.EQ.2) WRITE(LU,31) IORDRH
00748 30      FORMAT(1X,'WAVE_EQUATION : IORDRH=',1I6,' VALEUR NON PREVUE')
00749 31      FORMAT(1X,'WAVE_EQUATION: IORDRH=',1I6,' VALUE OUT OF RANGE')
00750         CALL PLANTE(1)
00751         STOP
00752       ENDIF
00753 !
00754 !     SAVES THE ORIGINAL MATRIX, BEFORE DIRICHLETS AND PRECONDITIONING
00755 !
00756       IF(OPTBAN.EQ.1.AND.OPT_HNEG.EQ.2.AND.NPTFR2.GT.0) THEN
00757         CALL OM('M=N     ',MAT2D%ADR(3)%P,MAT2D%ADR(1)%P,
00758      &          SVIDE,0.D0,MESH2D)
00759       ENDIF
00760 !
00761 !     DIRICHLET + SOLVES THE SYSTEM
00762 !
00763       CALL OSBD( 'X=X-Y   ' , HBOR , HN , HN , C , MESH2D )
00764       CALL DIRICH(DH,MAT2D%ADR(1)%P,SEM2D%ADR(1)%P,
00765      &            HBOR,LIHBOR%I,TRAV2,MESH2D,KENT,MSK,MASKPT)
00766       CALL SOLVE(DH,MAT2D%ADR(1)%P,SEM2D%ADR(1)%P,
00767      &           TRAV2,SLVPRO,INFOGR,MESH2D,MAT2D%ADR(2)%P)
00768       CALL OSBD( 'X=X+Y   ' , HBOR , HN , HN , C , MESH2D )
00769 !
00770 !     RECOVERS THE NEW DEPTH
00771 !
00772       CALL OS('X=Y+Z   ',X=H,Y=HN,Z=DH)
00773 !
00774 !     COMPLETES THE 2D INTERNAL FLUXES
00775 !     COMPATIBLE WITH CONTINUITY EQUATION
00776 !
00777 !     BOUNDARY FLUXES THAT SOLVE THE CONTINUITY EQUATION
00778 !     WHEN DEPTHS ARE PRESCRIBED
00779 !
00780       IF(OPTBAN.EQ.1.AND.OPT_HNEG.EQ.2.AND.NPTFR2.GT.0) THEN
00781         CALL MATVEC('X=AY    ',T2_01,MAT2D%ADR(3)%P,DH,1.D0,MESH2D)
00782         DO I=1,NPTFR2
00783           IF(LIHBOR%I(I).EQ.KENT) THEN
00784             FLBOR%R(I)=FLINT2%R(NBOR2%I(I))-T2_01%R(NBOR2%I(I))
00785           ENDIF
00786         ENDDO
00787       ENDIF
00788 !
00789 !     UNCOMMENT THIS LINE
00790 !     TO CHECK THAT SUM OF FLUINT = FLINT2 (IN TRIDW2)
00791 !
00792 !     CALL MATVEC ('X=X+CAY  ',FLINT2,MAT2D%ADR(2)%P,DH,-1.D0,MESH2D)
00793 !
00794 !=======================================================================
00795 !     PREPARING ZCONV AND GRAZCO
00796 !=======================================================================
00797 !
00798       CALL MAKE_ZCONV(ZCONV,GRAZCO,ZFLATS,DH,HN,ZF,TETAZCOMP,TETAH,
00799      &                NELEM2,OPTBAN,MESH2D%IKLE%I,MESH2D)
00800 !
00801 !=======================================================================
00802 !     COMPUTES NEW VELOCITIES
00803 !=======================================================================
00804 !
00805 !     GRADIENT OF DH (IN T2_02 AND T2_03)
00806 !
00807 !     COMPONENT X (IN T2_02)
00808 !
00809       CALL VECTOR
00810      & (T2_02,'=','GRADF          X',IELM2H,-GRAV*TETAH,DH,SVIDE,SVIDE,
00811      &  SVIDE,SVIDE,SVIDE, MESH2D, MSK, MASKEL)
00812       IF (NCSIZE.GT.1) CALL PARCOM(T2_02,2,MESH2D)
00813       CALL OS ( 'X=XY    ' ,X=T2_02,Y=UNSV2D)
00814 !
00815 !     COMPONENT Y (IN T2_03)
00816 !
00817       CALL VECTOR
00818      & (T2_03,'=','GRADF          Y',IELM2H,-GRAV*TETAH,DH,SVIDE,SVIDE,
00819      &  SVIDE,SVIDE,SVIDE, MESH2D, MSK, MASKEL)
00820       IF(NCSIZE.GT.1) CALL PARCOM(T2_03,2,MESH2D)
00821       CALL OS ( 'X=XY    ' ,X=T2_03,Y=UNSV2D)
00822 !
00823 !     THE NON COMPATIBLE PART OF THE LAPLACIAN FOR THE FREE SURFACE
00824 !     GRADIENT HAS NOT BEEN PUT IN UAUX, IT IS ADDED HERE IN A
00825 !     COMPATIBLE WAY TO COMPUTE U AND V
00826 !
00827       IF(ABS(1.D0-TETAZCOMP).GT.1.D-6) THEN
00828         C=-GRAV*(1.D0-TETAZCOMP)
00829         CALL OS('X=X+CY  ',X=T2_02,Y=GRADZN%ADR(1)%P,C=C)
00830         CALL OS('X=X+CY  ',X=T2_03,Y=GRADZN%ADR(2)%P,C=C)
00831       ENDIF
00832 !
00833 !     FINAL COMPUTATION OF U AND V
00834 !
00835       DO I=1,U%DIM1
00836         U%R(I)=UAUX(I)+DT*T2_02%R(MOD(I-1,NPOIN2)+1)*DM1%R(I)
00837         V%R(I)=VAUX(I)+DT*T2_03%R(MOD(I-1,NPOIN2)+1)*DM1%R(I)
00838       ENDDO
00839 !
00840 !     MODIFIES DM1 FOR USE IN PRECON, FLUX3D, ETC
00841 !
00842       CALL OS('X=CX    ',X=DM1,C=-DT*GRAV*TETAH*TETAU)
00843 !
00844 !     DIRICHLET TYPE BOUNDARY CONDITIONS
00845 !
00846 !     LATERAL BOUNDARY CONDITION
00847 !
00848       DO IPTFR3 = 1,NPTFR3
00849         IF(LIUBOL%I(IPTFR3).EQ.KENT.OR.
00850      &     LIUBOL%I(IPTFR3).EQ.KENTU.OR.LIUBOL%I(IPTFR3).EQ.KADH) THEN
00851           U%R(MESH3D%NBOR%I(IPTFR3)) = UBORL%R(IPTFR3)
00852         ENDIF
00853         IF(LIVBOL%I(IPTFR3).EQ.KENT.OR.
00854      &     LIVBOL%I(IPTFR3).EQ.KENTU.OR.LIVBOL%I(IPTFR3).EQ.KADH) THEN
00855           V%R(MESH3D%NBOR%I(IPTFR3)) = VBORL%R(IPTFR3)
00856         ENDIF
00857       ENDDO
00858 !
00859 !     BOTTOM AND SURFACE BOUNDARY CONDITION
00860 !
00861       DO IPOIN2 = 1,NPOIN2
00862         IF(LIUBOF%I(IPOIN2).EQ.KENT.OR.LIUBOF%I(IPOIN2).EQ.KADH) THEN
00863           U%R(IPOIN2) = UBORF%R(IPOIN2)
00864         ENDIF
00865         IF(LIVBOF%I(IPOIN2).EQ.KENT.OR.LIVBOF%I(IPOIN2).EQ.KADH) THEN
00866           V%R(IPOIN2) = VBORF%R(IPOIN2)
00867         ENDIF
00868         IF(LIUBOS%I(IPOIN2).EQ.KENT.OR.LIUBOS%I(IPOIN2).EQ.KADH) THEN
00869           U%R(NPOIN3-NPOIN2+IPOIN2) = UBORS%R(IPOIN2)
00870         ENDIF
00871         IF(LIVBOS%I(IPOIN2).EQ.KENT.OR.LIVBOS%I(IPOIN2).EQ.KADH) THEN
00872           V%R(NPOIN3-NPOIN2+IPOIN2) = VBORS%R(IPOIN2)
00873         ENDIF
00874       ENDDO
00875 !
00876 !     PROJECTION ON SOLID BOUNDARIES (KLOG)
00877 !
00878       IF(VELPROLAT) THEN
00879         CALL AIRWIK3(LIHBOR%I,U%R,V%R,MESH2D%XNEBOR%R,MESH2D%YNEBOR%R,
00880      &               NBOR2%I,NPTFR2,NPLAN,NPOIN2,KLOG)
00881       ENDIF
00882 !
00883 !     THIS SEQUENCE WILL BE DONE AFTER IF DYNAMIC PRESSURE HAS NOT BEEN
00884 !     COMPUTED HERE
00885 !
00886       IF(NONHYD) CALL OS ('X=Y     ', X=W , Y=WD  )
00887       IF(NONHYD.AND.DPWAVEQ) THEN
00888 !       BOUNDARY CONDITION ON FREE SURFACE STRONGLY ENFORCED
00889         IF(CLDYN) THEN
00890           CALL OV('X=Y     ',W%R(NPOIN3-NPOIN2+1:NPOIN3),DSSUDT%R,
00891      &                       DSSUDT%R,0.D0,NPOIN2)
00892           CALL OV('X=X+YZ  ',W%R(NPOIN3-NPOIN2+1:NPOIN3),
00893      &                       GRADZS%ADR(1)%P%R,
00894      &                       U%R(NPOIN3-NPOIN2+1:NPOIN3),0.D0,NPOIN2)
00895           CALL OV('X=X+YZ  ',W%R(NPOIN3-NPOIN2+1:NPOIN3),
00896      &                       GRADZS%ADR(2)%P%R,
00897      &                       V%R(NPOIN3-NPOIN2+1:NPOIN3),0.D0,NPOIN2)
00898         ENDIF
00899 !       BOUNDARY CONDITION ON BOTTOM STRONGLY ENFORCED
00900         IF(VELPROBOT) THEN
00901           IF(SIGMAG.OR.OPTBAN.EQ.1) THEN
00902             DO I=1,NPOIN2
00903               DO IP=0,IPBOT%I(I)
00904                 I3D=IP*NPOIN2+I
00905                 W%R(I3D)=GRADZF%ADR(1)%P%R(I)*U%R(I3D)
00906      &                  +GRADZF%ADR(2)%P%R(I)*V%R(I3D)
00907               ENDDO
00908             ENDDO
00909           ELSE
00910             DO I=1,NPOIN2
00911               W%R(I)=GRADZF%ADR(1)%P%R(I)*U%R(I)
00912      &              +GRADZF%ADR(2)%P%R(I)*V%R(I)
00913             ENDDO
00914           ENDIF
00915         ENDIF
00916       ENDIF
00917 !
00918 !     VELOCITIES OF FIRST FREE POINT COPIED BELOW
00919 !
00920       IF(SIGMAG.OR.OPTBAN.EQ.1) THEN
00921         DO I2=1,NPOIN2
00922           IF(IPBOT%I(I2).GT.0) THEN
00923             I1=I2+IPBOT%I(I2)*NPOIN2
00924 !           VALUE OF THE FIRST FREE POINT IS COPIED BELOW
00925             DO NP=0,IPBOT%I(I2)-1
00926               I3=I2+NP*NPOIN2
00927               U%R(I3)=U%R(I1)
00928               V%R(I3)=V%R(I1)
00929             ENDDO
00930           ENDIF
00931         ENDDO
00932       ENDIF
00933 !
00934 !     COMPUTING WCONV, BEFORE MODIFICATION OF W BY PRESSURE,
00935 !     IN ACCORDANCE WITH UCONV AND VCONV DONE BEFORE
00936 !
00937       IF(NONHYD) THEN
00938         DO I=1,NPOIN3
00939           WCONV%R(I)=TETAU*W%R(I)+(1.D0-TETAU)*WN%R(I)
00940         ENDDO
00941       ENDIF
00942 !
00943 !-----------------------------------------------------------------------
00944 !
00945       CALL VERMOY(U2D%R,V2D%R,U%R,V%R,2,Z,
00946      &            T3_01%R,T3_02%R,T3_03%R,1,NPLAN,NPOIN2,NPLAN,OPTBAN)
00947 !
00948 !-----------------------------------------------------------------------
00949 !
00950 !     CLASSICAL ADVECTION FIELD IS USED FOR CHARACTERISTICS
00951 !     IT IS REBUILT HERE IN UCONVC AND VCONVC FOR USE IN PRECON
00952 !
00953       IF(N_ADV(ADV_CAR).GT.0) THEN
00954         CALL OS( 'X=CY    ' , X=UCONVC, Y=UN , C=1.D0-TETAU )
00955         CALL OS( 'X=X+CY  ' , X=UCONVC, Y=U  , C=     TETAU )
00956         CALL OS( 'X=CY    ' , X=VCONVC, Y=VN , C=1.D0-TETAU )
00957         CALL OS( 'X=X+CY  ' , X=VCONVC, Y=V  , C=     TETAU )
00958       ENDIF
00959 !
00960 !     MERCATOR PROJECTION, ADVECTING FIELD MODIFIED FOR CHARACTERISTICS
00961 !
00962       IF(SPHERI) THEN
00963         DO IPLAN=1,NPLAN
00964           DO I=1,NPOIN2
00965             I3D=(IPLAN-1)*NPOIN2+I
00966             UCONVC%R(I3D)=UCONVC%R(I3D)/MESH3D%COSLAT%R(I)
00967             VCONVC%R(I3D)=VCONVC%R(I3D)/MESH3D%COSLAT%R(I)
00968           ENDDO
00969         ENDDO
00970       ENDIF
00971 !
00972 !-----------------------------------------------------------------------
00973 !
00974       RETURN
00975       END

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