propag.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\propag.f
00002 !
00198                      SUBROUTINE PROPAG
00199 !                    *****************
00200 !
00201      &(U,V,H,UCONV,VCONV,CONVV,H0,PATMOS,ATMOS,
00202      & HPROP,UN,VN,HN,UTILD,VTILD,HTILD,DH,DU,DV,DHN,VISC,VISC_S,FU,FV,
00203      & SMH,MESH,ZF,AM1,AM2,AM3,BM1,BM2,CM1,CM2,TM1,A23,A32,MBOR,
00204      & CV1,CV2,CV3,W1,UBOR,VBOR,AUBOR,HBOR,DIRBOR,
00205      & TE1,TE2,TE3,TE4,TE5,T1,T2,T3,T4,T5,T6,T7,T8,
00206      & LIMPRO,MASK,GRAV,ROEAU,CF,DIFVIT,IORDRH,IORDRU,LT,AT,DT,
00207      & TETAH,TETAHC,TETAU,TETAD,
00208      & AGGLOH,AGGLOU,KDIR,INFOGR,KFROT,ICONVF,
00209      & PRIVE,ISOUSI,BILMAS,MASSES,MASS_RAIN,YASMH,OPTBAN,CORCON,
00210      & OPTSUP,MSK,MASKEL,MASKPT,RO,ROVAR,
00211      & MAT,RHS,UNK,TB,S,BD,PRECCU,SOLSYS,CFLMAX,OPDVIT,OPTSOU,
00212      & NFRLIQ,SLVPRO,EQUA,VERTIC,ADJO,ZFLATS,TETAZCOMP,UDEL,VDEL,DM1,
00213      & ZCONV,COUPLING,FLBOR,BM1S,BM2S,CV1S,VOLU2D,V2DPAR,UNSV2D,
00214      & NDGA1,NDGB1,NWEIRS,NPSING,HFROT,FLULIM,YAFLULIM,RAIN,PLUIE,
00215      & MAXADV,OPTADV_VI)
00216 !
00217 !***********************************************************************
00218 ! TELEMAC2D   V7P0                                   21/08/2010
00219 !***********************************************************************
00220 !
00221 !
00222 !
00223 !
00224 !
00225 !
00226 !
00227 !
00228 !
00229 !
00230 !
00231 !
00232 !
00233 !
00234 !
00235 !
00236 !
00237 !
00238 !
00239 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00240 !| A23            |<->| MATRIX
00241 !| A32            |<->| MATRIX
00242 !| ADJO           |-->| IF YES : ADJOINT MODE
00243 !| AGGLOH         |-->| KEYWORD: 'MASS-LUMPING ON H'
00244 !| AGGLOU         |-->| KEYWORD: 'MASS-LUMPING ON VELOCITY'
00245 !| AM1            |<->| MATRIX APPLYING TO H
00246 !| AM2            |<->| MATRIX APPLYING TO U
00247 !| AM3            |<->| MATRIX APPLYING TO V
00248 !| AT             |-->| TIME IN SECONDS
00249 !| ATMOS          |-->| IF YES, ATMOSPHERIC PRESSURE IN PATMOS
00250 !| AUBOR          |<--| LAW OF FRICTION ON BOUNDARIES
00251 !|                |   | NUT*DU/DN=AUBOR*U+BUBOR
00252 !| BD             |---| ??????  NOT USED
00253 !| BILMAS         |-->| LOGICAL TRIGGERING A MASS BALANCE INFORMATION
00254 !| BM1S           |<->| MATRIX
00255 !| BM2            |<->| MATRIX
00256 !| BM2S           |<->| MATRIX
00257 !| CF             |<--| ADIMENSIONAL FRICTION COEFFICIENT
00258 !| CFLMAX         |<--| MAXIMUM CFL NUMBER (OBSERVED IN CURRENT TIME STEP)
00259 !| CM1            |<->| MATRIX
00260 !| CM2            |<->| MATRIX
00261 !| CONVV          |-->| ARRAY OF LOGICAL GIVING THE VARIABLES TO BE
00262 !|                |   | ADVECTED
00263 !|                |   | CONVV(1):U,V CONVV(2):H
00264 !| CORCON         |-->| CONTINUITY CORRECTION ON POINTS WITH
00265 !|                |   | IMPOSED DEPTH (COMPATIBLE FLUX IS COMPUTED)
00266 !| COUPLING       |-->| STRING WITH THE LIST OF COUPLED PROGRAMMES
00267 !| CV1            |<->| RIGHT-HAND SIDE OF LINEAR SYSTEM
00268 !| CV2            |<->| RIGHT-HAND SIDE OF LINEAR SYSTEM
00269 !| CV3            |<->| RIGHT-HAND SIDE OF LINEAR SYSTEM
00270 !| CV1S           |<->| RIGHT-HAND SIDE OF LINEAR SYSTEM
00271 !| DH             |<--| H(N+1)-H(N)
00272 !| DHN            |<--| H(N)-H(N-1)
00273 !| DIFVIT         |-->| IF YES, DIFFUSION OF VELOCITY
00274 !| DIRBOR         |<--| BLOCK WITH DIRICHLET BOUNDARY CONDITIONS
00275 !| DM1            |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
00276 !|                |   | IS DM1*GRAD(ZCONV), SEE SOLSYS.
00277 !| DT             |-->| TIME STEP
00278 !| DU             |<--| U(N+1)-U(N)
00279 !| DV             |<--| V(N+1)-V(N)
00280 !| EQUA           |-->| KEYWORD: 'EQUATIONS'
00281 !| FLBOR          |<--| FLUXES AT BOUNDARY POINTS
00282 !| FLULIM         |-->| FLUX LIMITATION
00283 !| FU             |<->| SOURCE TERMS ON VELOCITY U
00284 !| FV             |<->| SOURCE TERMS ON VELOCITY V
00285 !| GRAV           |-->| GRAVITY
00286 !| H0             |-->| REFERENCE DEPTH
00287 !| HBOR           |<--| PRESCRIBED BOUNDARY CONDITION ON DEPTH
00288 !| HN             |-->| DEPTH AT TIME T(N)
00289 !| HFROT          |-->| KEYWORD: 'DEPTH IN FRICTION TERMS'
00290 !| HPROP          |-->| PROPAGATION DEPTH
00291 !| HTILD          |-->| DEPTH AFTER ADVECTION
00292 !| ICONVF         |-->| TYPE OF ADVECTION: 4 INTEGERS
00293 !|                |   | ICONVF(1) : U AND V
00294 !|                |   | ICONVF(2) : H (MANDATORY VALUE = 5)
00295 !|                |   | ICONVF(3) : TRACERS
00296 !|                |   | ICONVF(4) : K AND EPSILON
00297 !| INFOGR         |-->| IF YES, INFORMATION ON GRADIENT
00298 !| IORDRH         |-->| ORDER OF INITIAL GUESS OF H
00299 !| IORDRU         |-->| ORDER OF INITIAL GUESS OF U
00300 !| ISOUSI         |-->| NUMBER OF SUB-ITERATION IN THE TIME-STEP
00301 !| KDIR           |-->| CONVENTION FOR DIRICHLET POINT
00302 !| KFROT          |-->| KEYWORD: 'LAW OF BOTTOM FRICTION'
00303 !| LIMPRO         |<->| BOUNDARY CONDITIONS FOR H, U V PER POINTS
00304 !|                |   | AND SEGMENTS
00305 !| LT             |-->| ITERATION NUMBER
00306 !| MASK           |-->| BLOCK OF MASKS FOR SEGMENTS :
00307 !|                |   | MASK(MSK1): 1. IF KDIR ON U 0. ELSE
00308 !|                |   | MASK(MSK2): 1. IF KDIR ON V 0. ELSE
00309 !|                |   | MASK(MSK3): 1. IF KDDL ON U 0. ELSE
00310 !|                |   | MASK(MSK4): 1. IF KDDL ON V 0. ELSE
00311 !|                |   | MASK(MSK6): 1. IF KNEU ON V 0. ELSE
00312 !|                |   | MASK(MSK7): 1. IF KOND 0. ELSE
00313 !|                |   | MASK(MSK9): 1. IF KDIR ON H (POINT)
00314 !| MASKEL         |-->| MASKING OF ELEMENTS
00315 !|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
00316 !| MASKPT         |-->| MASKING PER POINT.
00317 !|                |   | =1. : NORMAL   =0. : MASKED
00318 !| MASSES         |-->| MASS OF WATER ADDED BY SOURCE TERM
00319 !| MASS_RAIN      |-->| MASS ADDED BY RAIN OR EVAPORATION
00320 !| MAXADV         |-->| MAXIMUM NUMBER OF ITERATIONS FOR ADVECTION SCHEMES
00321 !| MAT            |<--| BLOCK OF MATRICES
00322 !| MBOR           |<--| BOUNDARY MATRIX
00323 !| MESH           |-->| MESH STRUCTURE
00324 !| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
00325 !| NFRLIQ         |-->| NUMBER OF LIQUID BOUNDARIES
00326 !| NPSING         |-->| NUMBER OF POINTS FOR EVERY SINGULARITY.
00327 !| NDGA1          |-->| NDGA1%ADR(I)%I(NP) : BOUNDARY NUMBER OF POINT NP
00328 !|                |   | OF WEIR I (side1)
00329 !| NDGB1          |-->| NDGB1%ADR(I)%I(NP) : BOUNDARY NUMBER OF POINT NP
00330 !|                |   | OF WEIR I (side2)
00331 !| NWEIRS         |-->| NUMBER OF SINGULARITIES
00332 !| OPDVIT         |-->| OPTION FOR DIFFUSION OF VELOCITIES
00333 !| OPTADV_VI      |-->| OPTION FOR THE ADVECTION SCHEME OF VELOCITIES
00334 !| OPTBAN         |-->| KEYWORD: 'OPTION FOR THE TREATMENT OF TIDAL FLATS'
00335 !| OPTSOU         |-->| KEYWORD: 'TYPE OF SOURCES'
00336 !| OPTSUP         |-->| KEYWORD: 'SUPG OPTION'
00337 !| PATMOS         |-->| ATMOSPHERIC PRESSURE
00338 !| PLUIE          |-->| RAIN OR EVAPORATION IN M/S IN A BIEF_OBJ
00339 !| PRECCU         |-->| KEYWORD: 'C-U PRECONDITIONING'
00340 !| PRIVE          |-->| BLOCK OF WORK BIEF_OBJ STRUCTURES
00341 !| RAIN           |-->| IF YES, RAIN OR EVAPORATION
00342 !| RHS            |<->| BLOCK OF PRIVATE BIEF_OBJ STRUCTURES
00343 !| RO             |-->| WATER DENSITY IF VARIABLE
00344 !| ROEAU          |-->| WATER DENSITY
00345 !| ROVAR          |-->| IF YES, VARIABLE WATER DENSITY.
00346 !| S              |-->| VOID STRUCTURE
00347 !| SLVPRO         |-->| SOLVER STRUCTURE FOR PROPAGATION
00348 !| SMH            |-->| SOURCE TERM IN CONTINUITY EQUATION
00349 !| SOLSYS         |-->| KEYWORD: 'TREATMENT OF THE LINEAR SYSTEM'
00350 !| T1             |<->| WORK BIEF_OBJ STRUCTURE
00351 !| T2             |<->| WORK BIEF_OBJ STRUCTURE
00352 !| T3             |<->| WORK BIEF_OBJ STRUCTURE
00353 !| T4             |<->| WORK BIEF_OBJ STRUCTURE
00354 !| T5             |<->| WORK BIEF_OBJ STRUCTURE
00355 !| T6             |<->| WORK BIEF_OBJ STRUCTURE
00356 !| T7             |<->| WORK BIEF_OBJ STRUCTURE
00357 !| T8             |<->| WORK BIEF_OBJ STRUCTURE
00358 !| TB             |<->| BLOCK WITH T1,T2,...
00359 !| TE1            |<->| WORK BIEF_OBJ STRUCTURE FOR ELEMENTS
00360 !| TE2            |<->| WORK BIEF_OBJ STRUCTURE FOR ELEMENTS
00361 !| TE3            |<->| WORK BIEF_OBJ STRUCTURE FOR ELEMENTS
00362 !| TE4            |<->| WORK BIEF_OBJ STRUCTURE FOR ELEMENTS
00363 !| TE5            |<->| WORK BIEF_OBJ STRUCTURE FOR ELEMENTS
00364 !| TETAD          |-->| IMPLICITATION ON DIFFUSION
00365 !| TETAH          |-->| IMPLICITATION OF H IN U EQUATION
00366 !| TETAHC         |-->| IMPLICITATION OF H IN CONTINUITY
00367 !| TETAU          |-->| IMPLICITATION OF U AND
00368 !| TM1            |<->| MATRIX
00369 !| UBOR           |<--| PRESCRIBED BOUNDARY CONDITION ON VELOCITY U
00370 !| VBOR           |<--| PRESCRIBED BOUNDARY CONDITION ON VELOCITY V
00371 !| UCONV          |-->| X-COMPONENT OF ADVECTION VELOCITY FIELD
00372 !| VCONV          |-->| Y-COMPONENT OF ADVECTION VELOCITY FIELD
00373 !| UDEL           |<--| COMPATIBLE X-COMPONENT OF ADVECTION VELOCITY FIELD
00374 !| UN             |<->| X-COMPONENT OF VELOCITY AT TIME T(N)
00375 !| VN             |<->| Y-COMPONENT OF VELOCITY AT TIME T(N)
00376 !| UNK            |<->| BLOCK OF UNKNOWNS
00377 !| UNSV2D         |-->| INVERSE OF INTEGRALS OF TEST FUNCTIONS
00378 !| UTILD          |-->| VELOCITY U IF ADVECTED BY CHARACTERISTICS
00379 !| V2DPAR         |-->| INTEGRAL OF TEST FUNCTIONS, ASSEMBLED IN PARALLEL
00380 !| VBOR           |-->| CONDITIONS AUX LIMITES SUR V.
00381 !| VDEL           |<--| COMPATIBLE Y-COMPONENT OF ADVECTION VELOCITY FIELD
00382 !| VERTIC         |-->| IF YES, THERE ARE VERTICAL STRUCTURES
00383 !| VISC           |-->| VISCOSITY COEFFICIENTS ALONG X,Y AND Z .
00384 !|                |   | IF P0 : PER ELEMENT
00385 !|                |   | IF P1 : PERR POINT
00386 !| VISC_S         |<->| WORK ARRAY FOR SAVING VISC
00387 !| VOLU2D         |-->| INTEGRAL OF TEST FUNCTIONS, NOT ASSEMBLED IN PARALLEL
00388 !| VTILD          |-->| VELOCITY V IF ADVECTED BY CHARACTERISTICS
00389 !| W1             |<->| WORK ARRAY
00390 !| YAFLULIM       |-->| IF, YES, FLULIM TAKEN INTO ACCOUNT
00391 !| YASMH          |-->| IF YES, SMH TAKEN INTO ACCOUNT
00392 !| ZCONV          |<--| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
00393 !|                |   | IS DM1*GRAD(ZCONV), SEE SOLSYS.
00394 !| ZF             |-->| ELEVATION OF BOTTOM
00395 !| ZFLATS         |<--| ELEVATION OF BOTTOM, MODIFIED FOR TIDAL FLATS
00396 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00397 !
00398       USE BIEF
00399       USE DECLARATIONS_TELEMAC, ONLY : ADV_CAR,ADV_SUP,ADV_NSC,ADV_PSI,
00400      &                                 ADV_PSI_NC,ADV_NSC_NC,ADV_LPO,
00401      &                                 ADV_NSC_TF,ADV_PSI_TF,ADV_LPO_TF,
00402      &                                 KDDL
00403       USE DECLARATIONS_TELEMAC2D, ONLY : TYPSEUIL,IT1,IT2
00404 !
00405       USE INTERFACE_TELEMAC2D, EX_PROPAG => PROPAG
00406 !
00407       IMPLICIT NONE
00408       INTEGER LNG,LU
00409       COMMON/INFO/LNG,LU
00410 !
00411 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00412 !
00413       INTEGER, INTENT(IN) :: LT,OPTSUP(4),KDIR,KFROT,ICONVF(4),NWEIRS
00414       INTEGER, INTENT(IN) :: IORDRH,IORDRU,ISOUSI,OPTBAN,OPTSOU,SOLSYS
00415       INTEGER, INTENT(IN) :: OPDVIT,NFRLIQ,HFROT,MAXADV,OPTADV_VI
00416       DOUBLE PRECISION, INTENT(IN)    :: TETAU,TETAD,TETAH,AGGLOH,AGGLOU
00417       DOUBLE PRECISION, INTENT(IN)    :: TETAHC,AT,DT,GRAV,ROEAU
00418       DOUBLE PRECISION, INTENT(IN)    :: TETAZCOMP
00419       DOUBLE PRECISION, INTENT(INOUT) :: CFLMAX,MASSES,MASS_RAIN
00420       LOGICAL, INTENT(IN) :: BILMAS,ATMOS,DIFVIT,INFOGR,CONVV(4),MSK
00421       LOGICAL, INTENT(IN) :: YASMH,ROVAR,PRECCU,VERTIC,ADJO,CORCON
00422       LOGICAL, INTENT(IN) :: YAFLULIM,RAIN
00423       TYPE(SLVCFG), INTENT(INOUT)     :: SLVPRO
00424       CHARACTER(LEN=20),  INTENT(IN)  :: EQUA
00425       CHARACTER(LEN=*) ,  INTENT(IN)  :: COUPLING
00426 !
00427 !  STRUCTURES OF VECTORS
00428 !
00429       TYPE(BIEF_OBJ), INTENT(IN)    :: NPSING,NDGA1,NDGB1
00430       TYPE(BIEF_OBJ), INTENT(IN)    :: UCONV,VCONV,SMH,UN,VN,HN
00431       TYPE(BIEF_OBJ), INTENT(IN)    :: VOLU2D,V2DPAR,UNSV2D,FLULIM
00432       TYPE(BIEF_OBJ), INTENT(INOUT) :: RO,UDEL,VDEL,DM1,ZCONV,FLBOR
00433       TYPE(BIEF_OBJ), INTENT(IN)    :: UTILD,VTILD,PATMOS,CF
00434       TYPE(BIEF_OBJ), INTENT(INOUT) :: U,V,H,CV1,CV2,CV3,PRIVE,DH,DHN
00435       TYPE(BIEF_OBJ), INTENT(INOUT) :: CV1S
00436       TYPE(BIEF_OBJ), INTENT(INOUT) :: DU,DV,FU,FV,VISC,VISC_S,HTILD
00437       TYPE(BIEF_OBJ), INTENT(INOUT) :: UBOR,VBOR,HBOR,AUBOR,LIMPRO
00438       TYPE(BIEF_OBJ), INTENT(IN)    :: MASKEL,MASKPT,ZF,PLUIE
00439       TYPE(BIEF_OBJ), INTENT(IN)    :: HPROP,H0
00440 !
00441 !     TE : BY ELEMENT               TE4,TE5 ONLY IF OPTBAN=3
00442       TYPE(BIEF_OBJ), INTENT(INOUT) :: TE1,TE2,TE3,TE4,TE5,ZFLATS
00443 !     T  : BY POINT
00444       TYPE(BIEF_OBJ), INTENT(INOUT) :: T1,T2,T3,T4,T5,T6,T7,T8
00445       TYPE(BIEF_OBJ), INTENT(INOUT) :: W1
00446 !     DUMMY STRUCTURE
00447       TYPE(BIEF_OBJ), INTENT(IN)    :: S
00448 !
00449 !-----------------------------------------------------------------------
00450 !
00451 !  STRUCTURES OF MATRICES
00452 !
00453       TYPE(BIEF_OBJ), INTENT(INOUT) :: AM1,AM2,AM3,BM1,BM2,CM1,CM2,TM1
00454       TYPE(BIEF_OBJ), INTENT(INOUT) :: A23,A32,MBOR,BM1S,BM2S
00455 !
00456 !-----------------------------------------------------------------------
00457 !
00458 !  STRUCTURES OF BLOCKS
00459 !
00460       TYPE(BIEF_OBJ), INTENT(INOUT) :: MASK,MAT,RHS,UNK,TB,BD,DIRBOR
00461 !
00462 !-----------------------------------------------------------------------
00463 !
00464 !  STRUCTURE OF MESH
00465 !
00466       TYPE(BIEF_MESH), INTENT(INOUT) :: MESH
00467 !
00468 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00469 !
00470       INTEGER I,IELMU,IELMH,UDIR,UDDL,UNEU,UNONNEU,IELEM,NELEM
00471       INTEGER I1,I2,I3,DIMLIM,DIMGLO,N,IOPT,DISCLIN
00472 !
00473       DOUBLE PRECISION Z(1),SL1,SL1U,C,FL1,FL2
00474 !
00475       DOUBLE PRECISION P_DSUM
00476       EXTERNAL         P_DSUM
00477 !
00478       LOGICAL MSKGRA
00479 !
00480       CHARACTER*16 FORMUL
00481 !
00482 !-----------------------------------------------------------------------
00483 !
00484       DIMLIM=LIMPRO%DIM1
00485       DIMGLO=MESH%GLOSEG%DIM1
00486 !
00487       IELMH=H%ELM
00488       IELMU=U%ELM
00489 !
00490 !  ADDRESSES OF THE ARRAYS IN THE MASKING BLOCK: MASK
00491 !
00492       UDIR = 1
00493 !     VDIR = 2
00494       UDDL = 3
00495 !     VDDL = 4
00496       UNEU = 5
00497 !     VNEU = 6
00498       UNONNEU = 8
00499 !
00500 !  CONVENTION FOR LINEAR DISCRETISATION
00501 !
00502       DISCLIN=11
00503 !
00504 !-----------------------------------------------------------------------
00505 !
00506 !  DIRICHLET BOUNDARY CONDITIONS FOR INCREASE IN H
00507 !
00508 !  HBOR = HBOR - HN (HBOR ON THE BOUNDARY, HN IN THE DOMAIN)
00509 !
00510       CALL OSBD( 'X=X-Y   ' , HBOR , HN , HN , C , MESH )
00511 !
00512 !=======================================================================
00513 !
00514 !    GRADIENT MATRICES FOR THE CONTINUITY EQUATION:
00515 !
00516 !    BM1 = - TETAU ( D(HN.F1)/DX * F2)
00517 !    BM2 = - TETAU ( D(HN.F1)/DY * F2)
00518 !
00519       IF(OPTBAN.EQ.3) THEN
00520 !       TAKES POROSITY INTO ACCOUNT
00521         CALL MATRIX(BM1,'M=N     ','MATFGR         X',IELMH,IELMU,
00522      &              TETAU,HPROP,S,S,S,S,S,MESH,.TRUE.,TE5)
00523         CALL MATRIX(BM2,'M=N     ','MATFGR         Y',IELMH,IELMU,
00524      &              TETAU,HPROP,S,S,S,S,S,MESH,.TRUE.,TE5)
00525 !
00526 !       MATRICES RESULTING FROM SUPG APPLIED TO THE ADVECTION TERM
00527 !
00528         IF(OPTSUP(2).EQ.1) THEN
00529           CALL KSUPG(TE1,TE2,1.D0,UCONV,VCONV,MESH)
00530           CALL MATRIX(BM1,'M=M+N   ','MATUGH         X',IELMH,IELMU,
00531      &                TETAU,HPROP,S,S,TE1,TE2,S,MESH,.TRUE.,TE5)
00532           CALL MATRIX(BM2,'M=M+N   ','MATUGH         Y',IELMH,IELMU,
00533      &                TETAU,HPROP,S,S,TE1,TE2,S,MESH,.TRUE.,TE5)
00534         ELSEIF(OPTSUP(2).EQ.2) THEN
00535           CALL MATRIX(BM1,'M=M+N   ','MATUGH         X',IELMH,IELMU,
00536      &                0.5*DT*TETAU,HPROP,S,S,UCONV,VCONV,S,
00537      &                MESH,.TRUE.,TE5)
00538           CALL MATRIX(BM2,'M=M+N   ','MATUGH         Y',IELMH,IELMU,
00539      &                0.5*DT*TETAU,HPROP,S,S,UCONV,VCONV,S,
00540      &                MESH,.TRUE.,TE5)
00541         ENDIF
00542 !
00543       ELSE
00544 !
00545         CALL MATRIX(BM1,'M=N     ','MATFGR         X',IELMH,IELMU,
00546      &              TETAU,HPROP,S,S,S,S,S,MESH,MSK,MASKEL)
00547         CALL MATRIX(BM2,'M=N     ','MATFGR         Y',IELMH,IELMU,
00548      &              TETAU,HPROP,S,S,S,S,S,MESH,MSK,MASKEL)
00549 !
00550 !       MATRICES RESULTING FROM SUPG APPLIED TO THE ADVECTION TERM
00551 !
00552         IF(OPTSUP(2).EQ.1) THEN
00553           CALL KSUPG(TE1,TE2,1.D0,UCONV,VCONV,MESH)
00554           CALL MATRIX(BM1,'M=M+N   ','MATUGH         X',IELMH,IELMU,
00555      &                TETAU,HPROP,S,S,TE1,TE2,S,MESH,MSK,MASKEL)
00556           CALL MATRIX(BM2,'M=M+N   ','MATUGH         Y',IELMH,IELMU,
00557      &                TETAU,HPROP,S,S,TE1,TE2,S,MESH,MSK,MASKEL)
00558         ELSEIF(OPTSUP(2).EQ.2) THEN
00559           CALL MATRIX(BM1,'M=M+N   ','MATUGH         X',IELMH,IELMU,
00560      &                0.5*DT*TETAU,HPROP,S,S,UCONV,VCONV,S,
00561      &                MESH,MSK,MASKEL)
00562           CALL MATRIX(BM2,'M=M+N   ','MATUGH         Y',IELMH,IELMU,
00563      &                0.5*DT*TETAU,HPROP,S,S,UCONV,VCONV,S,
00564      &                MESH,MSK,MASKEL)
00565         ENDIF
00566 !
00567       ENDIF
00568 !
00569 !=======================================================================
00570 !
00571 !   BUILDS THE DIFFUSION MATRIX :
00572 !
00573 !    TM1 =  VISC . (P03 * (DF1/DX * DF2/DX + DF1/DY * DF2/DY) )
00574 !
00575       IF(DIFVIT) THEN
00576 !
00577         IF(OPDVIT.EQ.2) THEN
00578 !         SAVES DIFFUSION
00579           CALL OS('X=Y     ',VISC_S,VISC,VISC,C)
00580 !         MULTIPLIES DIFFUSION BY HPROP
00581           CALL OV_2('X=XY    ',VISC%R,1,HPROP%R,1,Z,1,C,
00582      &                         VISC%MAXDIM1,VISC%DIM1)
00583           IF(VISC%DIM2.EQ.3) THEN
00584             CALL OV_2('X=XY    ',VISC%R,2,HPROP%R,1,Z,1,C,
00585      &                           VISC%MAXDIM1,VISC%DIM1)
00586             CALL OV_2('X=XY    ',VISC%R,3,HPROP%R,1,Z,1,C,
00587      &                           VISC%MAXDIM1,VISC%DIM1)
00588           ENDIF
00589         ENDIF
00590 !
00591         CALL MATRIX(TM1,'M=N     ','MATDIF          ',IELMU,IELMU,
00592      &              1.D0,S,S,S,VISC,S,S,MESH,MSK,MASKEL)
00593 !
00594         IF(OPDVIT.EQ.2) THEN
00595 !         MULTIPLIES THE MATRIX BY 1/HPROP
00596           CALL CPSTVC(HPROP,T4)
00597           DO I=1,HPROP%DIM1
00598 !           BEWARE: HIDDEN PARAMETER 1.D-2, NO DIFFUSION BELOW 1 CM DEPTH
00599 !                                           WITH THIS OPTION
00600             IF(HPROP%R(I).GT.1.D-2) THEN
00601               T4%R(I)=1.D0/HPROP%R(I)
00602             ELSE
00603               T4%R(I)=0.D0
00604             ENDIF
00605           ENDDO
00606           IF(T4%ELM.NE.IELMU) CALL CHGDIS(T4,T4%ELM,IELMU,MESH)
00607           CALL OM( 'M=X(M)  ' , TM1 , TM1 , S  , C , MESH )
00608           CALL OM( 'M=DM    ' , TM1 , TM1 , T4 , C , MESH )
00609 !         RESTORES DIFFUSION
00610           CALL OS('X=Y     ',VISC,VISC_S,VISC_S,C)
00611         ENDIF
00612 !
00613 !       'IF' ADDED ON 23/07/2002 BY JMH (MAY HAPPEN IN PARALLEL MODE)
00614 !
00615         IF(MESH%NELEB.GT.0) THEN
00616           CALL MATRIX(MBOR,'M=N     ','FMATMA          ',
00617      &                IELBOR(IELMU,1),IELBOR(IELMU,1),
00618      &               -1.D0,AUBOR,S,S,S,S,S,MESH,.TRUE.,MASK%ADR(UNEU)%P)
00619           CALL OM( 'M=M+N   ' , TM1 , MBOR , S , C , MESH )
00620         ENDIF
00621 !
00622 !       EXPLICIT PART DEALT WITH IN THE NEXT IF(DIFVIT...
00623 !
00624       ENDIF
00625 !
00626 !=======================================================================
00627 !
00628 !  COMPUTES THE FREE SURFACE ELEVATION (IN T8)
00629 !
00630       CALL OS( 'X=Y+Z   ' , T8 , HN , ZF , C )
00631 !
00632 !  OPTION 1 FOR THE TREATMENT OF TIDAL FLATS
00633 !  A MASK MAY BE USED FOR TIDAL FLATS ALTHOUGH THIS
00634 !  IS NOT THE MASKING OPTION.
00635 !  THIS MASK IS IN TE3
00636 !
00637       IF(OPTBAN.EQ.1.OR.OPTBAN.EQ.3) THEN
00638         CALL DECVRT(TE3,T8,ZF,MESH)
00639       ENDIF
00640 !
00641 !     FREE SURFACE GRADIENT
00642 !
00643       IF (OPTBAN.EQ.1.OR.OPTBAN.EQ.3) THEN
00644 !                   SL CORRECTED BY ELEMENT
00645         CALL CORRSL(ZFLATS,T8,ZF,MESH)
00646         CALL VECTOR(CV2,'=','GRADF          X',IELMU,
00647      &              -GRAV,ZFLATS,S,S,S,S,S,MESH,MSK,MASKEL)
00648         CALL VECTOR(CV3,'=','GRADF          Y',IELMU,
00649      &              -GRAV,ZFLATS,S,S,S,S,S,MESH,MSK,MASKEL)
00650 !       CORRSL DECLARES A DISCONTINUOUS ELEMENT, RESTORES BACK
00651       ELSE
00652         CALL VECTOR(CV2,'=','GRADF          X',IELMU,
00653      &              -GRAV,T8,S,S,S,S,S,MESH,MSK,MASKEL)
00654         CALL VECTOR(CV3,'=','GRADF          Y',IELMU,
00655      &              -GRAV,T8,S,S,S,S,S,MESH,MSK,MASKEL)
00656       ENDIF
00657 !
00658 !  ADDITIONAL GRADIENT TERMS : ATMOSPHERIC PRESSURE
00659 !                              VARIABLE DENSITY
00660 !
00661 !  THESE DRIVING TERMS SHOULD NOT BE ADDED IN TIDAL FLATS
00662 !
00663 !
00664       IF(ROVAR.OR.ATMOS) THEN
00665 !
00666 !  CASE WHERE THESE GRADIENTS SHOULD BE MASKED: MASKING REQUIRED OR
00667 !                                               OPTBAN=1
00668 !
00669       IF(MSK.OR.OPTBAN.EQ.1.OR.OPTBAN.EQ.3) THEN
00670         MSKGRA = .TRUE.
00671 !       IF OPTBAN=1 THE MASK USED HERE SHOULD HAVE VALUES
00672         IF(OPTBAN.NE.1.AND.OPTBAN.NE.3) THEN
00673           CALL OV('X=Y     ',TE3%R,MASKEL%R,MASKEL%R,C,TE3%DIM1)
00674         ENDIF
00675       ELSE
00676         MSKGRA = .FALSE.
00677       ENDIF
00678 !
00679 !     ATMOSPHERIC PRESSURE GRADIENT
00680 !
00681       IF(ATMOS) THEN
00682         CALL VECTOR(CV2,'+','GRADF          X',IELMU,
00683      &              -1.D0/ROEAU,PATMOS,S,S,S,S,S,MESH,MSKGRA,TE3)
00684         CALL VECTOR(CV3,'+','GRADF          Y',IELMU,
00685      &              -1.D0/ROEAU,PATMOS,S,S,S,S,S,MESH,MSKGRA,TE3)
00686       ENDIF
00687 !
00688 !     ADDITIONAL TERMS IF THE DENSITY IS VARIABLE
00689 !
00690       IF(ROVAR) THEN
00691 !
00692         CALL OS( 'X=X+C   ' , RO , S , S , -ROEAU )
00693 !
00694 !       PRESSURE BAROCLINE
00695 !
00696         CALL VECTOR(CV2,'+','GGRADF         X',IELMU,
00697      &              -0.5D0*GRAV/ROEAU,RO,HN,S,S,S,S,MESH,MSKGRA,TE3)
00698         CALL VECTOR(CV3,'+','GGRADF         Y',IELMU,
00699      &              -0.5D0*GRAV/ROEAU,RO,HN,S,S,S,S,MESH,MSKGRA,TE3)
00700         CALL OS( 'X=X+C   ' , RO , S , S , +ROEAU )
00701 !
00702       ENDIF
00703 !
00704       ENDIF
00705 !
00706 !=======================================================================
00707 !
00708 !    MASS MATRIX / DT
00709 !    AM1 WILL BE MODIFIED AT A LATER DATE IF SUPG IS USED ON H
00710 !
00711 !    AM1 = ( 1 / DT )  * (F1 * F2)
00712 !
00713       SL1 = 1.D0 / DT
00714       IF(OPTBAN.EQ.1.OR.OPTBAN.EQ.3) THEN
00715 !       LOCALLY LUMPED MASS MATRIX (TIDAL FLATS)
00716         FORMUL='MSLUMP          '
00717       ELSE
00718 !       NORMAL MASS MATRIX
00719         FORMUL='MATMAS          '
00720       ENDIF
00721       IF(OPTBAN.NE.3) THEN
00722         CALL MATRIX(AM1,'M=N     ',FORMUL,IELMH,IELMH,
00723      &              SL1,TE3,S,S,S,S,S,MESH,MSK,MASKEL)
00724       ELSE
00725         CALL MATRIX(AM1,'M=N     ',FORMUL,IELMH,IELMH,
00726      &              SL1,TE3,S,S,S,S,S,MESH,.TRUE.,TE5)
00727       ENDIF
00728 !
00729 !   MASS MATRIX FOR THE MOMENTUM EQUATION
00730 !
00731       IF(SOLSYS.EQ.1) THEN
00732 !
00733 !                           OPTBAN.NE.3 TO AVOID POROSITY IN AM2
00734       IF(IELMU.EQ.IELMH.AND.OPTBAN.NE.3) THEN
00735         CALL OM( 'M=N     ' , AM2 , AM1 , S , C , MESH )
00736       ELSE
00737         CALL MATRIX(AM2,'M=N     ',FORMUL,IELMU,IELMU,
00738      &              SL1,TE3,S,S,S,S,S,MESH,MSK,MASKEL)
00739       ENDIF
00740 !     MASS-LUMPING OF AM2 :
00741       IF(AGGLOU.GT.0.001D0) THEN
00742         CALL LUMP(T2,AM2,MESH,AGGLOU)
00743         CALL OM( 'M=CN    ' , AM2 , AM2 , S  , 1.D0-AGGLOU , MESH )
00744         CALL OM( 'M=M+D   ' , AM2 , AM2 , T2 , C           , MESH )
00745       ENDIF
00746 !
00747       ELSEIF(SOLSYS.EQ.2) THEN
00748 !
00749         IF(IELMU.NE.IELMH) THEN
00750           CALL VECTOR(AM2%D,'=','MASBAS          ',IELMU,
00751      &                SL1,S,S,S,S,S,S,MESH,MSK,MASKEL)
00752         ELSE
00753           CALL OS('X=CY    ',X=AM2%D,Y=VOLU2D,C=SL1)
00754         ENDIF
00755         AM2%TYPDIA='Q'
00756         AM2%TYPEXT='0'
00757 !
00758       ENDIF
00759 !
00760 ! MASS-LUMPING OF AM1 :
00761 !
00762       IF(AGGLOH.GT.0.001D0) THEN
00763         CALL LUMP(T1,AM1,MESH,AGGLOH)
00764         CALL OM( 'M=CN    ' , AM1 , AM1 , S , 1.D0-AGGLOH , MESH )
00765         CALL OM( 'M=M+D   ' , AM1 , AM1 , T1 , C          , MESH )
00766       ENDIF
00767 !
00768 ! END OF MASS-LUMPING
00769 !
00770 ! TEMPORARILY STORES THE LUMPED MASS MATRIX: AM2 IN AM3
00771 ! FOR THE COMPUTATION OF THE TERMS FROM TIME DERIVATIVES
00772 !
00773       IF(SOLSYS.EQ.1) THEN
00774         CALL OM( 'M=N     ' , AM3 , AM2 , S , C , MESH )
00775       ELSEIF(SOLSYS.EQ.2) THEN
00776         CALL OS('X=Y     ',X=AM3%D,Y=AM2%D)
00777         AM3%TYPDIA=AM2%TYPDIA
00778         AM3%TYPEXT=AM2%TYPEXT
00779       ENDIF
00780 !
00781 !=======================================================================
00782 !
00783 !   COMPUTES THE VECTORS IN THE SECOND MEMBERS
00784 !
00785 !   BEWARE : PAY A LOT OF ATTENTION TO PARAMETER LEGO (FALSE OR TRUE)
00786 !            IN CALLS TO MATVEC, VGRADF. IF LEGO = .FALSE. THE
00787 !            COMPUTATION RESULT GOES IN W1. IN THIS CASE, SHOULD ALWAYS
00788 !            USE OPERATIONS 'X=X+...' NOT TO ERASE WHAT WAS IN W1,
00789 !            AND END BY LEGO=.TRUE. FOR THE FINAL ASSEMBLY.
00790 !
00791 !
00792 !-----------------------------------------------------------------------
00793 !
00794 !     CV1 = CV1 + SL1U * ( BM1 * UN + BM2 * VN )
00795 !
00796       SL1U   = (TETAU-1.D0)/TETAU
00797 !                           TETAU : BM1 WAS BUILT WITH THIS COEFFICIENT
00798       IF(ABS(SL1U).GT.0.0001D0) THEN
00799         CALL MATVEC('X=CAY   ',CV1,BM1,UN,SL1U,MESH)
00800         CALL MATVEC('X=X+CAY ',CV1,BM2,VN,SL1U,MESH)
00801       ELSE
00802         CALL CPSTVC(H,CV1)
00803         CALL OS('X=0     ',X=CV1)
00804       ENDIF
00805 !
00806 !     SOURCE TERMS IN THE CONTINUITY EQUATION :
00807 !
00808       MASSES    = 0.D0
00809       MASS_RAIN = 0.D0
00810       IF(YASMH) THEN
00811         IF(OPTSOU.EQ.1) THEN
00812 !         STANDARD VERSION
00813 !         TAKING MASS-LUMPING INTO ACCOUNT
00814 !         COEFFICIENT DT TO COUNTERACT THE FACT THAT
00815 !         AM1 IS DONE WITH A COEFFICIENT 1/DT
00816           CALL MATVEC( 'X=CAY   ',T1,AM1,SMH,DT,MESH)
00817 !         WITHOUT TAKING MASS-LUMPING INTO ACCOUNT
00818 !         CALL VECTOR(T1,'=','MASVEC          ',IELMH,
00819 !    *                1.D0,SMH,S,S,S,S,S,MESH,MSK,MASKEL)
00820           CALL OS( 'X=X+Y   ' , X=CV1 , Y=T1 )
00821           IF(BILMAS) MASSES = DT * BIEF_SUM(T1)
00822         ELSE
00823 !         DIRAC VERSION
00824           CALL OS( 'X=X+Y   ' ,X=CV1,Y=SMH)
00825           IF(BILMAS) MASSES = DT * BIEF_SUM(SMH)
00826         ENDIF
00827 !       THE FOLLOWING LINE GOES IN BILAN
00828 !       IF (NCSIZE.GT.1) MASSES=P_DSUM (MASSES)
00829       ENDIF
00830 !
00831 !     SAME THING WITH RAIN-EVAPORATION (LIKE OPTSOU=1)
00832 !
00833       IF(RAIN) THEN
00834         CALL MATVEC( 'X=CAY   ',T1,AM1,PLUIE,DT,MESH)
00835         CALL OS( 'X=X+Y   ' , X=CV1 , Y=T1 )
00836         IF(BILMAS) THEN
00837           MASS_RAIN = DT * BIEF_SUM(T1)
00838           MASSES    = MASSES + MASS_RAIN
00839         ENDIF
00840       ENDIF
00841 !
00842 !  DEBUT DES CONVECTIONS DE U
00843 !
00844 !-----------------------------------------------------------------------
00845 !
00846 !     ADVECTION OF U AND V
00847 !
00848       IF(ICONVF(1).NE.ADV_LPO.AND.ICONVF(1).NE.ADV_LPO_TF.AND.
00849      &   ICONVF(1).NE.ADV_NSC.AND.ICONVF(1).NE.ADV_NSC_TF.AND.
00850      &   ICONVF(1).NE.ADV_PSI.AND.ICONVF(1).NE.ADV_PSI_TF     ) THEN
00851         CALL OS( 'X=CY    ' , X=T1 , Y=FU , C=DT )
00852         CALL OS( 'X=CY    ' , X=T2 , Y=FV , C=DT )
00853       ENDIF
00854 !
00855       IF(ICONVF(1).EQ.ADV_CAR.OR.(.NOT.CONVV(1))) THEN
00856 !
00857 !       ADVECTION WITH THE METHOD OF CHARACTERISTICS
00858 !
00859         CALL OS( 'X=X+Y   ' , X=T1 , Y=UTILD )
00860         CALL OS( 'X=X+Y   ' , X=T2 , Y=VTILD )
00861 !
00862 !------ SCHEMA SEMI-IMPLICITE CENTRE + S.U.P.G -----------------------
00863 !
00864       ELSEIF(ICONVF(1).EQ.ADV_SUP) THEN
00865 !
00866 !       CENTRED SEMI-IMPLICIT ADVECTION TERM: MATRIX
00867 !
00868         CALL MATRIX(CM2,'M=N     ','MATVGR          ',IELMU,IELMU,
00869      &              1.D0,S,S,S,UCONV,VCONV,VCONV,MESH,MSK,MASKEL)
00870 !
00871 !       SUPG CONTRIBUTION
00872 !
00873         IF(OPTSUP(1).EQ.1) THEN
00874 !         CLASSICAL SUPG
00875           CALL KSUPG(TE1,TE2,1.D0,UCONV,VCONV,MESH)
00876           CALL MATRIX(CM2,'M=M+N   ','MASUPG          ',IELMU,IELMU,
00877      &                1.D0,TE1,TE2,S,UCONV,VCONV,VCONV,
00878      &                MESH,MSK,MASKEL)
00879         ELSEIF(OPTSUP(1).EQ.2) THEN
00880 !         MODIFIED SUPG
00881           CALL MATRIX(CM2,'M=M+N   ','MAUGUG          ',IELMU,IELMU,
00882      &                0.5D0*DT,S,S,S,UCONV,VCONV,VCONV,
00883      &                MESH,MSK,MASKEL)
00884         ENDIF
00885 !
00886 !       END OF SUPG CONTRIBUTION
00887 !
00888 !       EXPLICIT SECOND MEMBER
00889 !
00890         IF(ABS(SL1U).GT.0.0001D0) THEN
00891           CALL MATVEC( 'X=X+CAY ',CV2,CM2,UN,TETAU-1.D0,MESH)
00892           CALL MATVEC( 'X=X+CAY ',CV3,CM2,VN,TETAU-1.D0,MESH)
00893         ENDIF
00894 !
00895 !       MATRIX : AM2 HAS A NON-SYMMETRICAL STRUCTURE
00896 !
00897         IF(AM2%TYPEXT.NE.'Q') THEN
00898           CALL OM( 'M=X(M)  ' , AM2 , AM2 , S , C , MESH )
00899         ENDIF
00900         CALL OM( 'M=M+CN  ' , AM2 , CM2 , S , TETAU , MESH )
00901 !
00902         CALL OS( 'X=X+Y   ' , X=T1 , Y=UN )
00903         CALL OS( 'X=X+Y   ' , X=T2 , Y=VN )
00904 !
00905 !-----------------------------------------------------------------
00906 !
00907 !  PSI SCHEME
00908 !
00909       ELSEIF(ICONVF(1).EQ.ADV_PSI_NC) THEN
00910 !
00911 !       STARTS COMPUTATION OF SECOND MEMBERS CV2 AND CV3
00912 !
00913         CALL VGFPSI(T6,IELMU,UCONV,VCONV,UN,DT,-1.D0,CFLMAX,
00914      &              T4,T5,MESH,MSK,MASKEL)
00915         CALL OS( 'X=X+Y   ' , X=CV2 , Y=T6 )
00916         CALL OS( 'X=X+Y   ' , X=T1  , Y=UN )
00917 !
00918         CALL VGFPSI(T6,IELMU,UCONV,VCONV,VN,DT,-1.D0,CFLMAX,
00919      &              T4,T5,MESH,MSK,MASKEL)
00920         CALL OS( 'X=X+Y   ' , X=CV3 , Y=T6 )
00921         CALL OS( 'X=X+Y   ' , X=T2  , Y=VN )
00922 !
00923 !------ SCHEMA SEMI-IMPLICITE N --------------------------------------
00924 !
00925       ELSEIF(ICONVF(1).EQ.ADV_NSC_NC) THEN
00926 !
00927 !       CENTRED SEMI-IMPLICIT ADVECTION TERM : MATRIX
00928 !
00929         CALL MATRIX(CM2,'M=N     ','MATVGR         N',IELMU,IELMU,
00930      &              1.D0,S,S,S,UCONV,VCONV,VCONV,MESH,MSK,MASKEL)
00931 !
00932 !       EXPLICIT SECOND MEMBER
00933 !
00934         IF(ABS(SL1U).GT.0.0001D0) THEN
00935           CALL MATVEC( 'X=X+CAY ',CV2,CM2,UN,TETAU-1.D0,MESH)
00936           CALL MATVEC( 'X=X+CAY ',CV3,CM2,VN,TETAU-1.D0,MESH)
00937         ENDIF
00938 !
00939 !       MATRIX : AM2 HAS A NON-SYMMETRICAL STRUCTURE
00940 !
00941         IF(AM2%TYPEXT.NE.'Q') THEN
00942           CALL OM( 'M=X(M)  ' , AM2 , AM2 , S , C , MESH )
00943         ENDIF
00944         CALL OM( 'M=M+CN  ' , AM2 , CM2 , S , TETAU , MESH )
00945         CALL OS( 'X=X+Y   ' , X=T1 , Y=UN )
00946         CALL OS( 'X=X+Y   ' , X=T2 , Y=VN )
00947 !
00948 !------ FINITE VOLUMES SCHEME --------------------------------------
00949 !
00950 !     NOTE: HERE THE CONTINUITY EQUATION IS NOT SOLVED
00951 !           BY H, HN AND UCONV,VCONV (UCONV, VCONV HAVE BEEN
00952 !           UPDATED AND H, HN ARE STILL UNCHANGED, SO THE FINAL H
00953 !           COMPUTED IN CVTRVF WILL NOT BE THE FINAL DEPTH OF THE
00954 !           TIMESTEP).
00955 !
00956       ELSEIF(ICONVF(1).EQ.ADV_LPO.OR.
00957      &       ICONVF(1).EQ.ADV_NSC.OR.
00958      &       ICONVF(1).EQ.ADV_PSI     ) THEN
00959 !
00960         IF(ICONVF(1).EQ.ADV_LPO.OR.ICONVF(1).EQ.ADV_NSC) IOPT=2
00961         IF(ICONVF(1).EQ.ADV_PSI) IOPT=3
00962 !       HERE YASMH=.FALSE. (SOURCES ACCOUNTED FOR IN FU)
00963         IF(TB%N.LT.22) THEN
00964           WRITE(LU,*) 'SIZE OF TB TOO SMALL IN PROPAG'
00965           CALL PLANTE(1)
00966           STOP
00967         ENDIF
00968 !       USING A COPY OF LIMPRO (IT MAY BE CHANGED BY CVTRVF)
00969         DO I=1,MESH%NPTFR
00970           IT1%I(I)=LIMPRO%I(DIMLIM+I)
00971           IT2%I(I)=LIMPRO%I(2*DIMLIM+I)
00972         ENDDO
00973         CALL CVTRVF(T1,UN,S,.FALSE.,.TRUE.,H,HN,
00974      &              HPROP,UCONV,VCONV,S,S,
00975      &              1,S,S,FU,S,.FALSE.,S,.FALSE.,UBOR,MASK,MESH,
00976      &              TB%ADR(13)%P,TB%ADR(14)%P,TB%ADR(15)%P,
00977      &              TB%ADR(16)%P,TB%ADR(17)%P,TB%ADR(18)%P,
00978      &              TB%ADR(19)%P,TB%ADR(20)%P,TB%ADR(21)%P,
00979      &              TB%ADR(22)%P,AGGLOH,TE1,DT,INFOGR,BILMAS,
00980      &              1,MSK,MASKEL,S,C,1,IT1%I,
00981      &              KDIR,KDDL,MESH%NPTFR,FLBOR,.FALSE.,
00982      &              VOLU2D,V2DPAR,UNSV2D,IOPT,TB%ADR(12)%P,MASKPT,
00983      &              RAIN,PLUIE,0.D0,OPTADV_VI)
00984         CALL CVTRVF(T2,VN,S,.FALSE.,.TRUE.,H,HN,
00985      &              HPROP,UCONV,VCONV,S,S,
00986      &              1,S,S,FV,S,.FALSE.,S,.FALSE.,VBOR,MASK,MESH,
00987      &              TB%ADR(13)%P,TB%ADR(14)%P,TB%ADR(15)%P,
00988      &              TB%ADR(16)%P,TB%ADR(17)%P,TB%ADR(18)%P,
00989      &              TB%ADR(19)%P,TB%ADR(20)%P,TB%ADR(21)%P,
00990      &              TB%ADR(22)%P,AGGLOH,TE1,DT,INFOGR,BILMAS,
00991      &              1,MSK,MASKEL,S,C,1,IT2%I,
00992      &              KDIR,KDDL,MESH%NPTFR,FLBOR,.FALSE.,
00993      &              VOLU2D,V2DPAR,UNSV2D,IOPT,TB%ADR(12)%P,MASKPT,
00994      &              RAIN,PLUIE,0.D0,OPTADV_VI)
00995         IF(IELMU.NE.11) THEN
00996           CALL CHGDIS(T1,DISCLIN,IELMU,MESH)
00997           CALL CHGDIS(T2,DISCLIN,IELMU,MESH)
00998         ENDIF
00999 !
01000 !------ SCHEMA VOLUMES FINIS AVEC BANCS DECOUVRANTS -------------------
01001 !
01002 !     NOTE: HERE THE CONTINUITY EQUATION IS NOT SOLVED
01003 !           BY H, HN AND UCONV,VCONV (UCONV, VCONV HAVE BEEN
01004 !           UPDATED AND H, HN ARE STILL UNCHANGED, SO THE FINAL H
01005 !           COMPUTED IN CVTRVF WILL NOT BE THE FINAL DEPTH OF THE
01006 !           TIMESTEP).
01007 !
01008       ELSEIF(ICONVF(1).EQ.ADV_LPO_TF.OR.
01009      &       ICONVF(1).EQ.ADV_NSC_TF.OR.
01010      &       ICONVF(1).EQ.ADV_PSI_TF     ) THEN
01011 !
01012         IF(ICONVF(1).EQ.ADV_LPO_TF) IOPT=2
01013         IF(ICONVF(1).EQ.ADV_NSC_TF) IOPT=2
01014         IF(ICONVF(1).EQ.ADV_PSI_TF) IOPT=3
01015 !       HERE YASMH=.FALSE. (SOURCES ACCOUNTED FOR IN FU)
01016         IF(TB%N.LT.22) THEN
01017           WRITE(LU,*) 'SIZE OF TB TOO SMALL IN PROPAG'
01018           CALL PLANTE(1)
01019           STOP
01020         ENDIF
01021 !       THIS IS EQUIVALENT TO TWO SUCCESSIVE CALLS TO CVTRVF_POS
01022 !       FOR U AND V
01023         CALL CVTRVF_POS_2(T1,UN,S,T2,VN,S,.FALSE.,.TRUE.,H,HN,
01024      &      HPROP,UCONV,VCONV,S,S,
01025      &      1,S,S,FU,FV,S,.FALSE.,S,S,.FALSE.,UBOR,VBOR,MASK,MESH,
01026      &      TB%ADR(13)%P,TB%ADR(14)%P,TB%ADR(15)%P,
01027      &      TB%ADR(16)%P,TB%ADR(17)%P,TB%ADR(18)%P,
01028      &      TB%ADR(19)%P,TB%ADR(20)%P,TB%ADR(21)%P,
01029      &      TB%ADR(22)%P,
01030      &      AGGLOH,TE1,DT,INFOGR,BILMAS,1,MSK,MASKEL,S,C,1,
01031      &      LIMPRO%I(1+DIMLIM:2*DIMLIM),
01032      &      LIMPRO%I(1+2*DIMLIM:3*DIMLIM),
01033      &      KDIR,3,MESH%NPTFR,FLBOR,.FALSE.,
01034      &      V2DPAR,UNSV2D,IOPT,TB%ADR(11)%P,TB%ADR(12)%P,MASKPT,
01035      &      MESH%GLOSEG%I(       1:  DIMGLO),
01036      &      MESH%GLOSEG%I(DIMGLO+1:2*DIMGLO),
01037      &      MESH%NBOR%I,2,FLULIM%R,YAFLULIM,RAIN,PLUIE,0.D0,0.D0,
01038      &      MAXADV)
01039 !                       2: HARDCODED OPTION
01040 !
01041         IF(IELMU.NE.11) THEN
01042           CALL CHGDIS(T1,DISCLIN,IELMU,MESH)
01043           CALL CHGDIS(T2,DISCLIN,IELMU,MESH)
01044         ENDIF
01045 !
01046 !-----------------------------------------------------------------
01047 !
01048       ELSE
01049 !
01050         IF(LNG.EQ.1) WRITE(LU,2002) ICONVF(1)
01051         IF(LNG.EQ.2) WRITE(LU,2003) ICONVF(1)
01052 2002    FORMAT(1X,'PROPAG : FORME DE LA CONVECTION DE U INCONNUE :',1I6)
01053 2003    FORMAT(1X,'PROPAG : UNKNOWN TYPE OF ADVECTION FOR U: ',1I6)
01054         CALL PLANTE(1)
01055         STOP
01056 !
01057       ENDIF
01058 !
01059 ! ADDS THE TERM AM3 * T1
01060 ! HERE AM3 MUST BE MASS/DT POSSIBLY WITH MASS-LUMPING ON U
01061 !
01062       IF(SOLSYS.EQ.1) THEN
01063         CALL MATVEC( 'X=X+AY  ',CV2,AM3,T1,C,MESH)
01064         CALL MATVEC( 'X=X+AY  ',CV3,AM3,T2,C,MESH)
01065       ELSEIF(SOLSYS.EQ.2) THEN
01066         CALL OS('X=X+YZ  ',X=CV2,Y=AM3%D,Z=T1)
01067         CALL OS('X=X+YZ  ',X=CV3,Y=AM3%D,Z=T2)
01068       ENDIF
01069 !
01070 !  END OF ADVECTION OF U AND V
01071 !
01072 !=======================================================================
01073 !
01074 !
01075 !     COMPUTES THE DIAGONAL MATRICES: - FU* (MASS OF THE BASES)
01076 !                                AND: - FV* (MASS OF THE BASES)
01077 !
01078       IF(KFROT.NE.0.OR.VERTIC) THEN
01079 !
01080 !       T3,T4 : BOTTOM FRICTION      T5,T6 : VERTICAL STRUCTURES
01081 !
01082         CALL FRICTI(T3,T4,T5,T6,UN,VN,HN,CF,MESH,T1,T2,VERTIC,
01083      &              UNSV2D,MSK,MASKEL,HFROT)
01084 !
01085 !       COMPUTES THE DIAGONAL MATRICES: - FU* (MASS OF THE BASES)
01086 !                                  AND: - FV* (MASS OF THE BASES)
01087 !
01088         CALL SLOPES(TE3,ZF,MESH)
01089         CALL VECTOR(T2,'=','MASBAS          ',IELMU,
01090      &              -1.D0,S,S,S,S,S,S,MESH,.TRUE.,TE3)
01091 !
01092         CALL OS( 'X=XY    ' , X=T3 , Y=T2 )
01093         CALL OS( 'X=XY    ' , X=T4 , Y=T2 )
01094 !
01095         IF(VERTIC) THEN
01096           CALL VECTOR(T2,'=','MASBAS          ',IELMU,
01097      &              -1.D0,S,S,S,S,S,S,MESH,.FALSE.,TE3)
01098           CALL OS( 'X=XY    ' , T5 , T2 , S , C )
01099           CALL OS( 'X=XY    ' , T6 , T2 , S , C )
01100           CALL OS( 'X=X+Y   ' , T3 , T5 , S , C )
01101           CALL OS( 'X=X+Y   ' , T4 , T6 , S , C )
01102         ENDIF
01103 !
01104       ELSE
01105 !
01106         CALL CPSTVC(U,T3)
01107         CALL CPSTVC(V,T4)
01108         CALL OS( 'X=0     ' , X=T3 )
01109         CALL OS( 'X=0     ' , X=T4 )
01110 !
01111         IF(OPTBAN.EQ.1) THEN
01112 !         SLOWING DOWN VELOCITIES ON TIDAL FLATS
01113 !         HAPPENS WHEN CALLED BY TELEMAC-3D WITH OPTT2D=1
01114           IF(IELMU.NE.IELMH) THEN
01115             CALL VECTOR(T2,'=','MASBAS          ',IELMU,
01116      &                  1.D0,S,S,S,S,S,S,MESH,.FALSE.,TE3)
01117             IF(NCSIZE.GT.1) CALL PARCOM(T2,2,MESH)
01118             CALL OS('X=Y     ',X=T5,Y=HPROP)
01119             IF(T5%ELM.NE.IELMU) CALL CHGDIS(T5,T5%ELM,IELMU,MESH)
01120             DO I=1,U%DIM1
01121 !             HIDDEN PARAMETER
01122               IF(T5%R(I).LT.1.D-3) THEN
01123                 T3%R(I)=10.D0*T2%R(I)/DT
01124                 T4%R(I)=T3%R(I)
01125               ENDIF
01126             ENDDO
01127           ELSE
01128             DO I=1,U%DIM1
01129 !             HIDDEN PARAMETER
01130               IF(HPROP%R(I).LT.1.D-3) THEN
01131                 T3%R(I)=10.D0*V2DPAR%R(I)/DT
01132                 T4%R(I)=T3%R(I)
01133               ENDIF
01134             ENDDO
01135           ENDIF
01136         ENDIF
01137 !
01138       ENDIF
01139 !
01140 !=======================================================================
01141 !
01142 !   COMPUTES THE MATRICES
01143 !
01144 !    AM1: ALREADY COMPUTED
01145 !
01146 !    AM2 = AM2 + TM1
01147 !
01148       IF(DIFVIT) THEN
01149 !
01150 !    TEST: IMPLICITATION OF TM1 DIAGONAL OF TM1 IN WAVE EQUATION
01151 !
01152 !
01153         IF(SOLSYS.EQ.2) THEN
01154 !
01155           CALL OS('X=X+CY  ',X=AM2%D,Y=TM1%D,C=TETAD)
01156           CALL OS('X=CX    ',X=TM1%D,C=1.D0-TETAD)
01157 !
01158           CALL MATVEC( 'X=X+CAY ',CV2,TM1,UN,-1.D0,MESH)
01159           CALL MATVEC( 'X=X+CAY ',CV3,TM1,VN,-1.D0,MESH)
01160 !
01161         ELSEIF(SOLSYS.EQ.1) THEN
01162 !
01163           IF(TETAD.LT.0.9999D0) THEN
01164             IF(TETAD.GT.0.0001D0) THEN
01165               IF(AM2%TYPEXT.EQ.'S'.AND.TM1%TYPEXT.EQ.'Q') THEN
01166                 CALL OM( 'M=X(M)  ' , AM2 , AM2 , S , C , MESH )
01167               ENDIF
01168               CALL OM( 'M=M+CN  ' , AM2 , TM1 , S , TETAD , MESH )
01169             ENDIF
01170 !           EXPLICIT PART :
01171             CALL MATVEC( 'X=X+CAY ',CV2,TM1,UN,TETAD-1.D0,MESH)
01172             CALL MATVEC( 'X=X+CAY ',CV3,TM1,VN,TETAD-1.D0,MESH)
01173           ELSE
01174 !           ENTIRELY IMPLICIT
01175             IF(AM2%TYPEXT.EQ.'S'.AND.TM1%TYPEXT.EQ.'Q') THEN
01176               CALL OM( 'M=X(M)  ' , AM2 , AM2 , S , C , MESH )
01177             ENDIF
01178             CALL OM( 'M=M+N   ' , AM2 , TM1 , S , C , MESH )
01179           ENDIF
01180 !
01181         ELSE
01182 !
01183           IF(LNG.EQ.1) WRITE(LU,*) 'PROPAG : MAUVAIS CHOIX POUR SOLSYS'
01184           IF(LNG.EQ.2) WRITE(LU,*) 'PROPAG : WRONG CHOICE FOR SOLSYS'
01185           CALL PLANTE(1)
01186           STOP
01187 !
01188         ENDIF
01189 !
01190       ENDIF
01191 !
01192 !    AM3 = AM2 (DIFFUSION HAS BEEN ADDED TO AM2, NOT TO AM3)
01193 !
01194       IF(SOLSYS.EQ.1) THEN
01195         CALL OM( 'M=N     ' , AM3 , AM2 , S , C , MESH )
01196       ELSEIF(SOLSYS.EQ.2) THEN
01197         CALL OS('X=Y     ',X=AM3%D,Y=AM2%D)
01198       ENDIF
01199 !
01200 !=======================================================================
01201 !
01202 !   DEFINA METHOD CORRECTED : RIGHT HAND SIDE MODIFIED
01203 !
01204 !     TM1 IS DONE AS AM1, BUT WITH TE4 INSTEAD OF TE5
01205 !
01206       IF(OPTBAN.EQ.3) THEN
01207         SL1 = 1.D0 / DT
01208 !       LOCALLY LUMPED MASS MATRIX (TIDAL FLATS)
01209         FORMUL='MSLUMP          '
01210         CALL MATRIX(TM1,'M=N     ',FORMUL,IELMH,IELMH,
01211      &              SL1,TE3,S,S,S,S,S,MESH,.TRUE.,TE4)
01212 !       MASS-LUMPING :
01213         IF(AGGLOH.GT.0.001D0) THEN
01214           CALL LUMP(T1,TM1,MESH,AGGLOH)
01215           CALL OM( 'M=CN    ' , TM1 , TM1 , S , 1.D0-AGGLOH , MESH )
01216           CALL OM( 'M=M+D   ' , TM1 , TM1 , T1 , C          , MESH )
01217         ENDIF
01218         CALL MATVEC('X=X+AY  ',CV1,TM1,HN,C,MESH)
01219       ENDIF
01220 !
01221 !=======================================================================
01222 !
01223 !  TAKES THE IMPLICIT SOURCE TERMS INTO ACCOUNT:
01224 !
01225       CALL OM( 'M=M+D   ' , AM2,AM2 , T3 , C , MESH )
01226       CALL OM( 'M=M+D   ' , AM3,AM3 , T4 , C , MESH )
01227 !
01228 !=======================================================================
01229 !
01230 !
01231 !  BOUNDARY TERMS
01232 !
01233 !     PARALLEL MODE : SUBDOMAINS MAY HAVE NO BOUNDARY POINT AT ALL
01234 !
01235       IF(MESH%NPTFR.GT.0) THEN
01236 !
01237 !  TAKES INTO ACCOUNT THE TERMS SUM(PSI H(N) U(N). N) AT THE BOUNDARY
01238 !  THESE TERMS SHOULD NOT BE TAKEN INTO ACCOUNT ON SOLID BOUNDARIES,
01239 !  HENCE THE USE OF MASK(*, 8) WHICH IS 0 FOR SEGMENTS OF TYPE KLOG.
01240 !
01241 !
01242       CALL VECTOR(FLBOR,'=','FLUBDF          ',IELBOR(IELMH,1),
01243      &            1.D0-TETAU,HPROP,S,S,UN,VN,S,
01244      &            MESH,.TRUE.,MASK%ADR(UNONNEU)%P)
01245 !
01246 !-----------------------------------------------------------------------
01247 !
01248 !  TAKES INTO ACCOUNT THE TERMS  SUM(PSI HN U(N+1) . N ) AT THE BOUNDARY
01249 !
01250 !  DIRICHLET CONDITIONS : U(N+1) = UBOR ; V(N+1) = VBOR
01251 !
01252 !  UDIR : 1 FOR A DIRICHLET SEGMENT FOR U, 0 OTHERWISE
01253 !
01254       CALL VECTOR(FLBOR,'+','FLUBDF          ',IELBOR(IELMH,1),
01255      &            TETAU,HPROP,S,S,UBOR,VBOR,S,
01256      &            MESH,.TRUE.,MASK%ADR(UDIR)%P)
01257 !
01258 !  TAKES INTO ACCOUNT THE TERMS  SUM(PSI HN U(N+1) . N ) AT THE BOUNDARY
01259 !
01260 !  FREE EXIT CONDITIONS
01261 !
01262 !  UDDL : 1 FOR A FREE EXIT SEGMENT FOR U, 0 OTHERWISE
01263 !
01264       CALL VECTOR(FLBOR,'+','FLUBDF          ',IELBOR(IELMH,1),
01265      &            TETAU,HPROP,S,S,UN,VN,S,
01266      &            MESH,.TRUE.,MASK%ADR(UDDL)%P)
01267 !
01268 !     WITH POROSITY
01269 !
01270       IF(OPTBAN.EQ.3) THEN
01271         DO I=1,MESH%NPTFR
01272           N=MESH%NELBOR%I(I)
01273 !         N MAY BE 0 IN PARALLELISM
01274           IF(N.GT.0) THEN
01275             FLBOR%R(I)=FLBOR%R(I)*TE5%R(N)
01276           ENDIF
01277         ENDDO
01278       ENDIF
01279 !
01280 !     END OF: IF(MESH%NPTFR.GT.0) THEN
01281       ENDIF
01282 !
01283 !     IMPOSING EQUALITY OF FLUXES ON EITHER SIDE OF A WEIR
01284 !     HERE ALL PROCESSORS DO ALL WEIRS...
01285 !
01286       IF(NWEIRS.GT.0.AND.TYPSEUIL.EQ.1) THEN
01287         DO N=1,NWEIRS
01288           DO I=1,NPSING%I(N)
01289             I1=NDGA1%ADR(N)%P%I(I)
01290             I2=NDGB1%ADR(N)%P%I(I)
01291             IF(I1.GT.0) THEN
01292               FL1=FLBOR%R(I1)
01293             ELSE
01294               FL1=0.D0
01295             ENDIF
01296             IF(I2.GT.0) THEN
01297               FL2=FLBOR%R(I2)
01298             ELSE
01299               FL2=0.D0
01300             ENDIF
01301 !           IN PARALLEL ASSEMBLED FLUXES
01302             IF(NCSIZE.GT.1) THEN
01303               FL1=P_DSUM(FL1)
01304               FL2=P_DSUM(FL2)
01305             ENDIF
01306 !           MULTIPLICATION FACTOR SO THAT BOTH FLUXES ARE EQUAL
01307             IF(ABS(FL1).GT.1.D-4.AND.ABS(FL2).GT.1.D-4) THEN
01308               IF(I1.GT.0) FLBOR%R(I1)= FLBOR%R(I1)*(FL1-FL2)*0.5D0/FL1
01309               IF(I2.GT.0) FLBOR%R(I2)=-FLBOR%R(I2)*(FL1-FL2)*0.5D0/FL2
01310             ELSE
01311               IF(I1.GT.0) FLBOR%R(I1)=0.D0
01312               IF(I2.GT.0) FLBOR%R(I2)=0.D0
01313             ENDIF
01314           ENDDO
01315         ENDDO
01316       ENDIF
01317 !
01318 !     BOUNDARY TERMS IN THE RIGHT HAND SIDE
01319 !
01320       IF(MESH%NPTFR.GT.0) THEN
01321         CALL OSDB( 'X=X-Y   ' , CV1 , FLBOR , FLBOR , C , MESH )
01322       ENDIF
01323 !
01324 ! END OF BOUNDARY TERMS
01325 !
01326 !-----------------------------------------------------------------------
01327 !
01328       IF(EQUA(1:10).EQ.'BOUSSINESQ') THEN
01329 !
01330 !       TAKES BOUSSINESQ TERMS INTO ACCOUNT
01331 !
01332         CALL ROTNE0(MESH,CM1,
01333      &              AM2,A23,A32,AM3,CV2,CV3,UN,VN,H0,MSK,MASKEL,S,DT)
01334 !
01335       ENDIF
01336 !
01337 !-----------------------------------------------------------------------
01338 !
01339 !  GRADIENT MATRICES
01340 !
01341 !     NOT USED WITH WAVE EQUATION
01342       IF(SOLSYS.EQ.1) THEN
01343 !
01344       CALL MATRIX(CM1,'M=N     ','MATGRA         X',IELMU,IELMH,
01345      &            TETAH*GRAV,S,S,S,S,S,S,MESH,MSK,MASKEL)
01346       CALL MATRIX(CM2,'M=N     ','MATGRA         Y',IELMU,IELMH,
01347      &            TETAH*GRAV,S,S,S,S,S,S,MESH,MSK,MASKEL)
01348 !
01349       ENDIF
01350 !
01351 !=======================================================================
01352 !
01353 ! INITIAL GUESS
01354 !
01355 ! FOR NOW U AND V ARE NOT MODIFIED, WHICH MEANS THAT U AND V STILL
01356 ! HOLD THE RESULT OF THE LAST SOLVED SYSTEM
01357 !
01358       IF(IORDRH.EQ.0) THEN
01359 !
01360         CALL OS( 'X=0     ' , X=DH )
01361 !
01362       ELSEIF(IORDRH.EQ.1) THEN
01363 !
01364         IF(LT.EQ.1.AND.ISOUSI.EQ.1) THEN
01365           CALL OS( 'X=0     ' , X=DH )
01366         ENDIF
01367 !
01368       ELSEIF(IORDRH.EQ.2) THEN
01369 !
01370         IF(LT.EQ.1.AND.ISOUSI.EQ.1) THEN
01371           CALL OS( 'X=0     ' , X=DH )
01372           CALL OS( 'X=0     ' , X=DHN)
01373         ENDIF
01374         IF (LT.GT.2) CALL OS( 'X=CX    ' , DH , S , S , 2.D0 )
01375         CALL OS( 'X=X-Y   ' , DH , DHN , S , C )
01376 !       STORES DH(N) IN DH(N-1)
01377         CALL OS( 'X=X+Y   ' , DHN , DH   , S , C     )
01378         CALL OS( 'X=CX    ' , DHN , DHN  , S , 0.5D0 )
01379 !
01380       ELSE
01381 !
01382         IF(LNG.EQ.1) WRITE(LU,30) IORDRH
01383         IF(LNG.EQ.2) WRITE(LU,31) IORDRH
01384 30      FORMAT(1X,'PROPAG : IORDRH=',1I6,' VALEUR NON PREVUE')
01385 31      FORMAT(1X,'PROPAG : IORDRH=',1I6,' VALUE OUT OF RANGE')
01386         CALL PLANTE(1)
01387         STOP
01388 !
01389       ENDIF
01390 !
01391       IF(IORDRU.EQ.0) THEN
01392 !
01393         CALL OS( 'X=0     ' , X=U )
01394         CALL OS( 'X=0     ' , X=V )
01395 !
01396       ELSEIF(IORDRU.EQ.1) THEN
01397 !
01398 !       U = UN AND V = VN ALREADY DONE
01399 !
01400       ELSEIF(IORDRU.EQ.2) THEN
01401 !
01402         IF(LT.EQ.1.AND.ISOUSI.EQ.1) THEN
01403           CALL OS( 'X=0     ' , X=DU )
01404           CALL OS( 'X=0     ' , X=DV )
01405         ENDIF
01406         IF(ISOUSI.EQ.1) THEN
01407           CALL OS( 'X=Y+Z   ' , X=U , Y=UN , Z=DU )
01408           CALL OS( 'X=Y+Z   ' , X=V , Y=VN , Z=DV )
01409         ENDIF
01410 !
01411       ELSE
01412 !
01413         IF(LNG.EQ.1) WRITE(LU,32) IORDRU
01414         IF(LNG.EQ.2) WRITE(LU,33) IORDRU
01415 32      FORMAT(1X,'PROPAG : IORDRU=',1I6,' VALEUR NON PREVUE')
01416 33      FORMAT(1X,'PROPAG : IORDRU=',1I6,' VALUE OUT OF RANGE')
01417         CALL PLANTE(1)
01418         STOP
01419 !
01420       ENDIF
01421 !
01422 !=======================================================================
01423 !
01424       IF(SOLSYS.EQ.2) THEN
01425 !
01426         IF(NCSIZE.GT.1) THEN
01427           CALL PARCOM(AM2%D,2,MESH)
01428           CALL PARCOM(AM3%D,2,MESH)
01429         ENDIF
01430 !       INVERSION OF AM2%D AND AM3%D (WILL BE USED AGAIN AT A LATER STAGE)
01431         CALL OS( 'X=1/Y   ' , AM2%D , AM2%D , AM2%D , C ,2,0.D0,1.D-6)
01432         CALL OS( 'X=1/Y   ' , AM3%D , AM3%D , AM3%D , C ,2,0.D0,1.D-6)
01433 !
01434 !       ADDS THE "DIFFUSION" MATRIX TO AM1
01435 !
01436 !       HERE AM2%D HAS ALREADY BEEN INVERSED
01437         IF(IELMH.EQ.IELMU) THEN
01438 !         WANT TO DIVIDE BY (1/DT + FROT) WHICH IS IN AM2%D EXCEPT
01439 !         THAT IT IS PROJECTED ON THE BASES (IN AM2%D); THEREFORE HAS
01440 !         TO MULTIPLY BY THE MASS OF THE BASES
01441           CALL OS('X=CYZ   ',X=DM1,Y=V2DPAR,Z=AM2%D,C=-TETAU*TETAH*GRAV)
01442         ELSE
01443 !         TAKE HERE THE MASS OF THE BASES FOR U
01444           CALL VECTOR(T4,'=','MASBAS          ',IELMU,
01445      &                1.D0,S,S,S,S,S,S,MESH,.FALSE.,TE3)
01446           IF(NCSIZE.GT.1) CALL PARCOM(T4,2,MESH)
01447           CALL OS('X=CYZ   ',X=DM1,Y=T4,Z=AM2%D,C=-TETAU*TETAH*GRAV)
01448           CALL CHGDIS(DM1,IELMU,IELMH,MESH)
01449         ENDIF
01450 !
01451         CALL MATRIX(AM1,'M=M+N   ','MATDIFUV        ',IELMH,IELMH,
01452      &              -1.D0,S,S,S,DM1,HPROP,S,MESH,MSK,MASKEL)
01453 !
01454 !       NEW SECOND MEMBER CV1
01455 !
01456         IF(ABS(TETAZCOMP-1.D0).LT.1.D-6) THEN
01457           CALL OS( 'X=YZ    ' , X=T2 , Y=CV2 , Z=AM2%D )
01458           CALL OS( 'X=YZ    ' , X=T3 , Y=CV3 , Z=AM3%D )
01459         ELSE
01460           CALL OS( 'X=Y     ' , X=T4 , Y=CV2 )
01461           CALL OS( 'X=Y     ' , X=T5 , Y=CV3 )
01462 !         FREE SURFACE GRADIENT
01463           IF(OPTBAN.EQ.1.OR.OPTBAN.EQ.3) THEN
01464             CALL VECTOR(T4,'+','GRADF          X',IELMU,
01465      &      (1.D0-TETAZCOMP)*GRAV,ZFLATS,S,S,S,S,S,MESH,MSK,MASKEL)
01466             CALL VECTOR(T5,'+','GRADF          Y',IELMU,
01467      &      (1.D0-TETAZCOMP)*GRAV,ZFLATS,S,S,S,S,S,MESH,MSK,MASKEL)
01468           ELSE
01469             CALL VECTOR(T4,'+','GRADF          X',IELMU,
01470      &      (1.D0-TETAZCOMP)*GRAV,T8,S,S,S,S,S,MESH,MSK,MASKEL)
01471             CALL VECTOR(T5,'+','GRADF          Y',IELMU,
01472      &      (1.D0-TETAZCOMP)*GRAV,T8,S,S,S,S,S,MESH,MSK,MASKEL)
01473           ENDIF
01474           CALL OS( 'X=YZ    ' , X=T2 , Y=T4 , Z=AM2%D )
01475           CALL OS( 'X=YZ    ' , X=T3 , Y=T5 , Z=AM3%D )
01476         ENDIF
01477 !
01478         IF(NCSIZE.GT.1) THEN
01479           CALL PARCOM(T2,2,MESH)
01480           CALL PARCOM(T3,2,MESH)
01481         ENDIF
01482 !
01483 !       TAKES THE BOUNDARY CONDITIONS INTO ACCOUNT
01484 !       ERROR IN GRAD(DH)
01485         DO I=1,MESH%NPTFR
01486           IF(LIMPRO%I(I+DIMLIM).EQ.KDIR) THEN
01487             T2%R(MESH%NBOR%I(I)) = UBOR%R(I)
01488           ENDIF
01489           IF(LIMPRO%I(I+2*DIMLIM).EQ.KDIR) THEN
01490             T3%R(MESH%NBOR%I(I)) = VBOR%R(I)
01491           ENDIF
01492         ENDDO
01493 !
01494 !       REMEMBER THAT COEFFICIENT TETAU IS IN BM1 AND BM2
01495 !       AND THAT SUPG UPWINDING IS ALSO IN BM1 AND BM2
01496 !       OTHERWISE COULD BE INCLUDED IN HUGRADP BELOW
01497         CALL MATVEC('X=X-AY  ',CV1,BM1,T2,C,MESH)
01498         CALL MATVEC('X=X-AY  ',CV1,BM2,T3,C,MESH)
01499 !
01500         IF(ABS(TETAZCOMP-1.D0).GT.1.D-6) THEN
01501           FORMUL='HUGRADP3        '
01502 !                        3: U AND V, HERE T2 AND T3 NOT TAKEN
01503           CALL OS('X=CY    ',X=T1,Y=DM1,
01504      &                       C=(1.D0-TETAZCOMP)/TETAH/TETAU)
01505           IF(OPTBAN.EQ.1.OR.OPTBAN.EQ.3) THEN
01506             CALL VECTOR(CV1,'+',FORMUL,IELMH,TETAU,
01507      &                  HPROP,T1,ZFLATS,T2,T3,T3,MESH,MSK,MASKEL)
01508           ELSE
01509             CALL VECTOR(CV1,'+',FORMUL,IELMH,TETAU,
01510      &                  HPROP,T1,T8    ,T2,T3,T3,MESH,MSK,MASKEL)
01511           ENDIF
01512         ENDIF
01513 !
01514         CALL OS('X=CY    ',X=UDEL,Y=T2,C=TETAU)
01515         CALL OS('X=CY    ',X=VDEL,Y=T3,C=TETAU)
01516         CALL OS('X=X+CY  ',X=UDEL,Y=UN,C=1.D0-TETAU)
01517         CALL OS('X=X+CY  ',X=VDEL,Y=VN,C=1.D0-TETAU)
01518 !
01519       ENDIF
01520 !
01521 !=======================================================================
01522 !
01523 !     AT THIS STAGE, A23 AND A32 EQUAL 0
01524 !     THE MATRICES HAVE A VALUE ONLY AFTER A DIAGONAL-BLOCK
01525 !     PRECONDITIONING OU WITH BOUSSINESQ
01526 !
01527       IF(EQUA(1:10).NE.'BOUSSINESQ') THEN
01528         A23%TYPDIA='0'
01529         A32%TYPDIA='0'
01530         A23%TYPEXT='0'
01531         A32%TYPEXT='0'
01532       ENDIF
01533 !
01534 !=======================================================================
01535 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
01536 !
01537 !     IN ADJOINT MODE : THE SYSTEM IS NOT SOLVED, RETURN HERE
01538 !
01539 !
01540       IF(ADJO) RETURN
01541 !
01542 !
01543 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
01544 !=======================================================================
01545 !
01546 !   DIRICHLET BOUNDARY CONDITIONS:
01547 !
01548       IF(SOLSYS.EQ.1) THEN
01549 !
01550         IF(CORCON.AND.NFRLIQ.GT.0) THEN
01551 !
01552 !         SAVES THE CONTINUITY EQUATION
01553           CALL OM( 'M=N     ' , TM1  , AM1 , S , C , MESH )
01554           CALL OM( 'M=N     ' , BM1S , BM1 , S , C , MESH )
01555           CALL OM( 'M=N     ' , BM2S , BM2 , S , C , MESH )
01556           CALL OS( 'X=Y     ' ,X=CV1S,Y=CV1)
01557 !
01558         ENDIF
01559 !
01560         CALL DIRICH(UNK,MAT,RHS,DIRBOR,LIMPRO%I,TB,MESH,KDIR,MSK,MASKPT)
01561 !
01562       ENDIF
01563 !
01564 !  SOLVES THE OBTAINED SYSTEM:
01565 !
01566 !+++++++++++++++++++++++++++++++++++
01567 !  SPECIAL PRECONDITIONING H-U     +
01568 !+++++++++++++++++++++++++++++++++++
01569 !
01570       IF(PRECCU) THEN
01571 !
01572 !     PREPARES THE DIAGONALS FOR PRECONDITIONING:
01573 !
01574 !     HTILD: D1 (MUST BE SIMPLIFIED BY GRAV THERE? )
01575       CALL OS('X=+(Y,C)' , X=HTILD , Y=HN , C=0.D0 )
01576       CALL OS('X=CX    ' , X=HTILD , C=4.D0/GRAV )
01577       CALL OS('X=SQR(Y)' , X=HTILD , Y=HTILD )
01578       CALL OS('X=+(Y,C)' , X=HTILD , Y=HTILD , C=2.D0/GRAV )
01579 !     T1: D2 (NOT KEPT)
01580       CALL OS('X=1/Y   ' , X=T1 , Y=HTILD )
01581       CALL OS('X=CX    ' , X=T1 , C=4.D0*TETAH/TETAU )
01582 !
01583 !     MODIFIES THE SECOND MEMBER
01584 !
01585       CALL OS('X=XY    ' , X=CV1 , Y=T1 )
01586 !
01587 !     MODIFIES THE VARIABLE DH
01588 !
01589       CALL OS('X=Y/Z   ' , X=DH , Y=DH , Z=HTILD )
01590 !
01591 !     PRECONDITIONING FOR AM1
01592 !
01593       IF(AM1%TYPEXT.EQ.'S') CALL OM( 'M=X(M)  ' , AM1,AM1 ,S,C,MESH )
01594       CALL OM( 'M=DM    ' , AM1 , AM1 , T1    , C , MESH )
01595       CALL OM( 'M=MD    ' , AM1 , AM1 , HTILD , C , MESH )
01596 !
01597 !     PRECONDITIONING FOR BM1 AND BM2 :
01598 !
01599       CALL OM( 'M=DM    ' , BM1 , BM1 , T1 , C , MESH )
01600       CALL OM( 'M=DM    ' , BM2 , BM2 , T1 , C , MESH )
01601 !
01602 !     PRECONDITIONING FOR CM1 AND CM2 :
01603 !
01604       CALL OM( 'M=MD    ' , CM1 , CM1 , HTILD , C , MESH )
01605       CALL OM( 'M=MD    ' , CM2 , CM2 , HTILD , C , MESH )
01606 !
01607       ENDIF
01608 !
01609 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01610 !  END OF SPECIAL PRECONDITIONING H-U                               +
01611 !  EXCEPT FOR RECOVERY OF THE DH VARIABLE (SEE AFTER CALL TO SOLV09)+
01612 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01613 !
01614 !
01615       IF(SOLSYS.EQ.1) THEN
01616 !
01617 !       CLASSICAL METHOD
01618 !
01619 !       CASE OF THE BLOCK-DIAGONAL PRECONDITIONING
01620 !       A23 AND A32 WILL BE USED, THEY ARE INITIALISED.
01621 !       ALREADY DONE WITH COUPLED BOUSSINESQ
01622         IF(EQUA(1:10).NE.'BOUSSINESQ') THEN
01623           IF(3*(SLVPRO%PRECON/3).EQ.SLVPRO%PRECON) THEN
01624             A23%TYPDIA='Q'
01625             A32%TYPDIA='Q'
01626             A23%TYPEXT='Q'
01627             A32%TYPEXT='Q'
01628             CALL OM( 'M=0     ' , A23,A23 , S,C , MESH )
01629             CALL OM( 'M=0     ' , A32,A32 , S,C , MESH )
01630             IF (AM2%TYPEXT.EQ.'S') THEN
01631              CALL OM( 'M=X(M)  ' , AM2,AM2 , S,C , MESH )
01632             ENDIF
01633             IF (AM3%TYPEXT.EQ.'S') THEN
01634              CALL OM( 'M=X(M)  ' , AM3,AM3 , S,C , MESH )
01635             ENDIF
01636           ENDIF
01637         ENDIF
01638 !
01639 !       FOR CHECKING PROBLEMS IN PARALLEL
01640 !       NON ASSEMBLED VECTORS MUST BE ASSEMBLED FIRST
01641 !       AND THEIR DOT PRODUCT (THAT SHOULD BE EQUAL IN SCALAR AND PARALLEL)
01642 !       COMPUTED WITH P_DOTS. ALL THIS IS DONE BY APPDOTS
01643 !
01644 !       CALL APPDOTS(AM1%D,MESH)
01645 !       CALL APPDOTS(AM2%D,MESH)
01646 !       CALL APPDOTS(AM3%D,MESH)
01647 !       CALL APPDOTS(BM1%D,MESH)
01648 !       CALL APPDOTS(BM2%D,MESH)
01649 !       CALL APPDOTS(CM1%D,MESH)
01650 !       CALL APPDOTS(CM2%D,MESH)
01651 !       CALL APPDOTS(CV1,MESH)
01652 !       CALL APPDOTS(CV2,MESH)
01653 !       CALL APPDOTS(CV3,MESH)
01654 !
01655         CALL SOLVE(UNK,MAT,RHS,TB,SLVPRO,INFOGR,MESH,TM1)
01656 !
01657       ELSEIF(SOLSYS.EQ.2) THEN
01658 !
01659 !       GENERALISED WAVE EQUATION
01660 !
01661 !       SYSTEM IN H
01662 !
01663 !       STORES AM1 IN BM1 AND CV1 IN BM2%D
01664 !
01665         IF(CORCON.AND.NFRLIQ.GT.0) THEN
01666           CALL OM('M=N     ',BM1,AM1,S,C,MESH)
01667           CALL OS('X=Y     ',X=BM2%D,Y=CV1)
01668         ENDIF
01669 !
01670         CALL DIRICH(DH,AM1,CV1,HBOR,LIMPRO%I,TB,MESH,KDIR,MSK,MASKPT)
01671         CALL SOLVE(DH,AM1,CV1,TB,SLVPRO,INFOGR,MESH,TM1)
01672 !
01673         NELEM=MESH%NELEM
01674         DO IELEM=1,NELEM
01675           ZCONV%R(IELEM        )=DH%R(MESH%IKLE%I(IELEM        ))
01676           ZCONV%R(IELEM+  NELEM)=DH%R(MESH%IKLE%I(IELEM+  NELEM))
01677           ZCONV%R(IELEM+2*NELEM)=DH%R(MESH%IKLE%I(IELEM+2*NELEM))
01678         ENDDO
01679         IF(ABS(1.D0-TETAZCOMP).GT.1.D-6) THEN
01680           C=(1.D0-TETAZCOMP)/TETAH
01681           IF(OPTBAN.EQ.1.OR.OPTBAN.EQ.3) THEN
01682 !           FREE SURFACE PIECE-WISE LINEAR IN ZFLATS
01683             CALL OS('X=X+CY  ',X=ZCONV,Y=ZFLATS,C=C)
01684           ELSE
01685 !           FREE SURFACE LINEAR
01686             DO IELEM=1,NELEM
01687               I1=MESH%IKLE%I(IELEM        )
01688               I2=MESH%IKLE%I(IELEM+  NELEM)
01689               I3=MESH%IKLE%I(IELEM+2*NELEM)
01690               ZCONV%R(IELEM        )=ZCONV%R(IELEM        )+
01691      &        C*(T8%R(I1))
01692               ZCONV%R(IELEM+  NELEM)=ZCONV%R(IELEM+  NELEM)+
01693      &        C*(T8%R(I2))
01694               ZCONV%R(IELEM+2*NELEM)=ZCONV%R(IELEM+2*NELEM)+
01695      &        C*(T8%R(I3))
01696             ENDDO
01697           ENDIF
01698         ENDIF
01699 !
01700 !       SYSTEMS IN U AND V
01701 !
01702         CALL VECTOR(CV2,'+','GRADF          X',IELMU,
01703      &              -GRAV*TETAH,DH,S,S,S,S,S,MESH,MSK,MASKEL)
01704         CALL VECTOR(CV3,'+','GRADF          Y',IELMU,
01705      &              -GRAV*TETAH,DH,S,S,S,S,S,MESH,MSK,MASKEL)
01706         IF(NCSIZE.GT.1) THEN
01707           CALL PARCOM(CV2,2,MESH)
01708           CALL PARCOM(CV3,2,MESH)
01709         ENDIF
01710 !                                      AM2%D AND AM3%D ALREADY INVERSED
01711         CALL OS('X=YZ    ',X=U,Y=CV2,Z=AM2%D)
01712         CALL OS('X=YZ    ',X=V,Y=CV3,Z=AM3%D)
01713 !
01714         DO I=1,MESH%NPTFR
01715           IF(LIMPRO%I(I+DIMLIM).EQ.KDIR) THEN
01716             U%R(MESH%NBOR%I(I)) = UBOR%R(I)
01717           ENDIF
01718           IF(LIMPRO%I(I+2*DIMLIM).EQ.KDIR) THEN
01719             V%R(MESH%NBOR%I(I)) = VBOR%R(I)
01720           ENDIF
01721         ENDDO
01722 !
01723 !       FINAL CORRECTION OF BOUNDARY FLUXES FOR ELEVATION IMPOSED
01724 !       BOUNDARIES (TO SOLVE THE CONTINUITY EQUATION)
01725 !
01726         IF(CORCON.AND.NFRLIQ.GT.0) THEN
01727           CALL MATVEC('X=X+CAY ',BM2%D,BM1,DH,-1.D0,MESH)
01728           DO I=1,MESH%NPTFR
01729             IF(LIMPRO%I(I).EQ.KDIR) THEN
01730               FLBOR%R(I)=FLBOR%R(I)+BM2%D%R(MESH%NBOR%I(I))
01731             ENDIF
01732           ENDDO
01733         ENDIF
01734 !
01735       ENDIF
01736 !
01737 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01738 ! SPECIAL PRECONDITIONING H-U : RECOVERY OF THE DH VARIABLE
01739 !
01740       IF(PRECCU) CALL OS('X=XY    ' , X=DH , Y=HTILD )
01741 !
01742       IF(CORCON.AND.SOLSYS.EQ.1.AND.NFRLIQ.GT.0) THEN
01743 !
01744 !       FINAL CORRECTION OF BOUNDARY FLUXES FOR ELEVATION IMPOSED
01745 !       BOUNDARIES (TO SOLVE THE CONTINUITY EQUATION)
01746 !
01747         CALL MATVEC('X=X+CAY ',CV1S,TM1,DH,-1.D0,MESH)
01748         CALL MATVEC('X=X+CAY ',CV1S,BM1S,U,-1.D0,MESH)
01749         CALL MATVEC('X=X+CAY ',CV1S,BM2S,V,-1.D0,MESH)
01750         DO I=1,MESH%NPTFR
01751           IF(LIMPRO%I(I).EQ.KDIR) THEN
01752             FLBOR%R(I)=FLBOR%R(I)+CV1S%R(MESH%NBOR%I(I))
01753           ENDIF
01754         ENDDO
01755       ENDIF
01756 !
01757 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01758 !
01759 !  FROM DH TO H
01760 !
01761       CALL OS( 'X=Y+Z   ' , X=H , Y=DH , Z=HN )
01762 !
01763 !  HBOR = HBOR + HN
01764 !  HBOR IS USED AGAIN IN SUB-ITERATIONS
01765 !
01766       CALL OSBD( 'X=X+Y   ' , HBOR , HN , S , C , MESH )
01767 !
01768 !  STORES THE RELATIVE CHANGE IN SPEEDS
01769 !
01770       IF(IORDRU.EQ.2) THEN
01771         CALL OS( 'X=Y-Z   ' , X=DU , Y=U , Z=UN )
01772         CALL OS( 'X=Y-Z   ' , X=DV , Y=V , Z=VN )
01773       ENDIF
01774 !
01775 !  COMPATIBLE VELOCITY FIELD IN CONTINUITY EQUATION
01776 !
01777       IF(SOLSYS.EQ.1) THEN
01778         CALL OS ('X=CY    ',X=UDEL,Y=U ,C=     TETAU)
01779         CALL OS ('X=CY    ',X=VDEL,Y=V ,C=     TETAU)
01780         CALL OS ('X=X+CY  ',X=UDEL,Y=UN,C=1.D0-TETAU)
01781         CALL OS ('X=X+CY  ',X=VDEL,Y=VN,C=1.D0-TETAU)
01782       ENDIF
01783 !
01784 !-----------------------------------------------------------------------
01785 !
01786       RETURN
01787       END

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