mt02pp_star.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\mt02pp_star.f
00002 !
00117                      SUBROUTINE MT02PP_STAR
00118 !                    **********************
00119 !
00120      &(T,XM,XMUL,SF,SG,SH,F,G,H,X,Y,Z,SURFAC,IKLE,NELEM,NELMAX,INCHYD,
00121      & FORMUL,NPLAN)
00122 !
00123 !***********************************************************************
00124 ! BIEF   V6P3                                   21/08/2010
00125 !***********************************************************************
00126 !
00127 !
00128 !
00129 !
00130 !
00131 !
00132 !
00133 !
00134 !
00135 !
00136 !
00137 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00138 !| F              |-->| FUNCTION USED IN THE FORMULA
00139 !| FORMUL         |-->| FORMULA DESCRIBING THE RESULTING MATRIX
00140 !| G              |-->| FUNCTION USED IN THE FORMULA
00141 !| H              |-->| FUNCTION USED IN THE FORMULA
00142 !| IKLE           |-->| CONNECTIVITY TABLE.
00143 !| INCHYD         |-->| IF YES, TREATS HYDROSTATIC INCONSISTENCIES
00144 !| NELEM          |-->| NUMBER OF ELEMENTS
00145 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00146 !| NPLAN          |-->| NUMBER OF PLANES IN THE MESH OF PRISMS
00147 !| SF             |-->| STRUCTURE OF FUNCTIONS F
00148 !| SG             |-->| STRUCTURE OF FUNCTIONS G
00149 !| SH             |-->| STRUCTURE OF FUNCTIONS H
00150 !| SURFAC         |-->| AREA OF 2D ELEMENTS
00151 !| T              |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
00152 !| X              |-->| ABSCISSAE OF POINTS, PER ELEMENT
00153 !| Y              |-->| ORDINATES OF POINTS, PER ELEMENT
00154 !| Z              |-->| ELEVATIONS OF POINTS, PER POINT
00155 !| XM             |<->| OFF-DIAGONAL TERMS
00156 !| XMUL           |-->| COEFFICIENT FOR MULTIPLICATION
00157 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00158 !
00159       USE BIEF, EX_MT02PP_STAR => MT02PP_STAR
00160 !
00161       IMPLICIT NONE
00162       INTEGER LNG,LU
00163       COMMON/INFO/LNG,LU
00164 !
00165 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00166 !
00167       INTEGER, INTENT(IN) :: NELEM,NELMAX,NPLAN
00168       INTEGER, INTENT(IN) :: IKLE(NELMAX,6)
00169 !
00170 !                                                              15 OR 30
00171       DOUBLE PRECISION, INTENT(INOUT) :: T(NELMAX,6),XM(NELMAX,*)
00172       DOUBLE PRECISION, INTENT(IN)    :: SURFAC(NELMAX)
00173 !
00174       DOUBLE PRECISION, INTENT(IN)    :: XMUL
00175       DOUBLE PRECISION, INTENT(IN)    :: F(*),G(*),H(*)
00176 !
00177 !     STRUCTURES DE F,G,H
00178 !
00179       TYPE(BIEF_OBJ), INTENT(IN)      :: SF,SG,SH
00180 !
00181       DOUBLE PRECISION, INTENT(IN)    :: X(NELMAX,6),Y(NELMAX,6),Z(*)
00182 !
00183       LOGICAL, INTENT(IN)             :: INCHYD
00184       CHARACTER(LEN=16), INTENT(IN)   :: FORMUL
00185 !
00186 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00187 !
00188 !     DECLARATIONS SPECIFIQUES A CE SOUS-PROGRAMME
00189 !
00190       DOUBLE PRECISION H1,H2,H3
00191       DOUBLE PRECISION NF1,NF2,NF3,NG1,NG2,NG3,NH1,NH2,NH3
00192       DOUBLE PRECISION XS06,XS24,NPX,NPY,NPXL,NPXU,XS03
00193       DOUBLE PRECISION NPYL,NPYU,GX(3),GY(3),NUXMOY,NUYMOY
00194 !
00195       DOUBLE PRECISION X2,X3,Y2,Y3,SOMVX,X2X3,X2AUX,X3AUX,NPX2,NPY2
00196       DOUBLE PRECISION SOMVY,Y2Y3,Y2AUX,Y3AUX,AUX
00197 !
00198       INTEGER I1,I2,I3,I4,I5,I6,IELEM,II4,II5,II6,NPOU0
00199 !
00200       DOUBLE PRECISION EPSI
00201       DATA EPSI/1.D-4/
00202 !
00203 !***********************************************************************
00204 !
00205       XS03  =XMUL/3.D0
00206       XS06  =XMUL/6.D0
00207       XS24  =XMUL/24.D0
00208 !
00209       IF(SF%ELM.NE.41) THEN
00210         IF (LNG.EQ.1) WRITE(LU,1000) SF%ELM
00211         IF (LNG.EQ.2) WRITE(LU,1001) SF%ELM
00212 1000    FORMAT(1X,'MT02PP_STAR (BIEF) : TYPE DE F NON PREVU : ',I6)
00213 1001    FORMAT(1X,'MT02PP_STAR (BIEF) : TYPE OF F NOT IMPLEMENTED: ',I6)
00214         CALL PLANTE(1)
00215         STOP
00216       ENDIF
00217       IF(SG%ELM.NE.41) THEN
00218         IF (LNG.EQ.1) WRITE(LU,2000) SG%ELM
00219         IF (LNG.EQ.2) WRITE(LU,2001) SG%ELM
00220 2000    FORMAT(1X,'MT02PP_STAR (BIEF) : TYPE DE G NON PREVU : ',I6)
00221 2001    FORMAT(1X,'MT02PP_STAR (BIEF) : TYPE OF G NOT IMPLEMENTED: ',I6)
00222         CALL PLANTE(1)
00223         STOP
00224       ENDIF
00225       IF(SH%ELM.NE.41) THEN
00226         IF (LNG.EQ.1) WRITE(LU,3000) SH%ELM
00227         IF (LNG.EQ.2) WRITE(LU,3001) SH%ELM
00228 3000    FORMAT(1X,'MT02PP_STAR (BIEF) : TYPE DE H NON PREVU : ',I6)
00229 3001    FORMAT(1X,'MT02PP_STAR (BIEF) : TYPE OF H NOT IMPLEMENTED: ',I6)
00230         CALL PLANTE(1)
00231         STOP
00232       ENDIF
00233 !
00234 !     POUR LE TRAITEMENT DE LA VISCOSITE VERTICALE P0 SUR LA VERTICALE
00235 !     VOIR VISCLM DE TELEMAC-3D
00236 !
00237       IF(SH%DIMDISC.EQ.0) THEN
00238 !       VISCOSITE VERTICALE P1
00239         NPOU0=SH%DIM1/NPLAN
00240       ELSEIF(SH%DIMDISC.EQ.4111) THEN
00241 !       VISCOSITE VERTICALE P0 SUR LA VERTICALE (VOIR II4,5,6)
00242         NPOU0=0
00243       ELSE
00244         IF (LNG.EQ.1) WRITE(LU,4000) SH%DIMDISC
00245         IF (LNG.EQ.2) WRITE(LU,4001) SH%DIMDISC
00246 4000    FORMAT(1X,'MT02PP_STAR (BIEF) : DIMDISC NON PREVU : ',I6)
00247 4001    FORMAT(1X,'MT02PP_STAR (BIEF): DIMDISC NOT IMPLEMENTED: ',I6)
00248         CALL PLANTE(1)
00249         STOP
00250       ENDIF
00251 !
00252       IF(FORMUL(11:13).EQ.'MON') THEN
00253 !
00254         IF(XMUL.LT.0.D0) THEN
00255           IF (LNG.EQ.1) WRITE(LU,4002) XMUL
00256           IF (LNG.EQ.2) WRITE(LU,4003) XMUL
00257 4002      FORMAT(1X,'MT02PP (BIEF) : XMUL NEGATIF NON PREVU : ',G16.7)
00258 4003      FORMAT(1X,'MT02PP (BIEF): NEGATIVE XMUL EXCLUDED: ',G16.7)
00259           CALL PLANTE(1)
00260           STOP
00261         ENDIF
00262 !
00263 !-----------------------------------------------------------------------
00264 !
00265 !     VERSION WITH MONOTONICITY
00266 !
00267 !-----------------------------------------------------------------------
00268 !
00269       DO IELEM=1,NELEM
00270 !
00271         I1=IKLE(IELEM,1)
00272         I2=IKLE(IELEM,2)
00273         I3=IKLE(IELEM,3)
00274         I4=IKLE(IELEM,4)
00275         I5=IKLE(IELEM,5)
00276         I6=IKLE(IELEM,6)
00277 !       SUIVANT NPOU0, II4 SERA I4 OU I1, ETC.
00278         II4=I1+NPOU0
00279         II5=I2+NPOU0
00280         II6=I3+NPOU0
00281 !
00282         H1=Z(I4)-Z(I1)
00283         H2=Z(I5)-Z(I2)
00284         H3=Z(I6)-Z(I3)
00285 !
00286         IF((INCHYD.AND.MAX(Z(I1),Z(I2),Z(I3)).GT.
00287      &                 MIN(Z(I4),Z(I5),Z(I6)))    .OR.
00288      &      H1.LT.EPSI.OR.H2.LT.EPSI.OR.H3.LT.EPSI ) THEN
00289           NH1=0.D0
00290           NH2=0.D0
00291           NH3=0.D0
00292         ELSE
00293 !         SUIVANT LES CAS (II4=I1 OU I4, ETC.)
00294 !         VARIANTE AVEC VISCOSITE VERTICALE LINEAIRE (II4=I4)
00295 !         VARIANTE AVEC VISCOSITE VERTICALE P0 SUR LA VERTICALE (II4=I1)
00296           NH1=0.5D0*(H(I1)+H(II4))/H1
00297           NH2=0.5D0*(H(I2)+H(II5))/H2
00298           NH3=0.5D0*(H(I3)+H(II6))/H3
00299         ENDIF
00300 !
00301 !       DIFFUSION ALONG PLANES
00302 !
00303 !       X2 = X(I2)-X(I1)
00304 !       X3 = X(I3)-X(I1)
00305 !       Y2 = Y(I2)-Y(I1)
00306 !       Y3 = Y(I3)-Y(I1)
00307         X2 = X(IELEM,2)
00308         X3 = X(IELEM,3)
00309         Y2 = Y(IELEM,2)
00310         Y3 = Y(IELEM,3)
00311 !
00312         X2X3  = X2*X3
00313         Y2Y3  = Y2*Y3
00314         X2AUX = X2X3-X2**2
00315         X3AUX = X3**2-X2X3
00316         Y2AUX = Y2Y3-Y2**2
00317         Y3AUX = Y3**2-Y2Y3
00318 !
00319         AUX = XS24 / SURFAC(IELEM)
00320 !
00321 !       2D DIFFUSION, LOWER LEVEL (SEE MT02AA)
00322 !       LIKE IN 2D BUT MULTIPLIED BY DELTA(Z)/2
00323         SOMVX = ( F(I1)*H1+F(I2)*H2+F(I3)*H3 ) * AUX
00324         SOMVY = ( G(I1)*H1+G(I2)*H2+G(I3)*H3 ) * AUX
00325 !
00326 !       OFF-DIAGONAL TERMS FOR POINTS OF LOWER LEVEL
00327 !       WITH MONOTONICITY ENSURED
00328 !
00329         XM(IELEM, 1) = MIN(0.D0,- SOMVY * X3AUX - SOMVX * Y3AUX)
00330         XM(IELEM, 2) = MIN(0.D0,  SOMVY * X2AUX + SOMVX * Y2AUX)
00331         XM(IELEM, 6) = MIN(0.D0,- SOMVY * X2X3  - SOMVX * Y2Y3 )
00332 !
00333 !       2D DIFFUSION, UPPER LEVEL
00334         SOMVX = ( F(I4)*H1+F(I5)*H2+F(I6)*H3 ) * AUX
00335         SOMVY = ( G(I4)*H1+G(I5)*H2+G(I6)*H3 ) * AUX
00336 !
00337 !       OFF-DIAGONAL TERMS FOR POINTS OF UPPER LEVEL
00338 !       WITH MONOTONICITY ENSURED
00339 !
00340         XM(IELEM,13) = MIN(0.D0,- SOMVY * X3AUX - SOMVX * Y3AUX)
00341         XM(IELEM,14) = MIN(0.D0,  SOMVY * X2AUX + SOMVX * Y2AUX)
00342         XM(IELEM,15) = MIN(0.D0,- SOMVY * X2X3  - SOMVX * Y2Y3 )
00343 !
00344 !       HERE TERMS 2 AND 3 HAVE BEEN REMOVED FOR MONOTONICITY
00345 !       CROSSED TERMS CANCELLED
00346 !
00347         XM(IELEM,04)=0.D0
00348         XM(IELEM,05)=0.D0
00349         XM(IELEM,07)=0.D0
00350         XM(IELEM,09)=0.D0
00351         XM(IELEM,10)=0.D0
00352         XM(IELEM,11)=0.D0
00353 !
00354 !       TERM 4
00355 !
00356 !       HERE SOME HORIZONTAL TERMS HAVE BEEN REMOVED BECAUSE
00357 !       HORIZONTAL FLUXES THROUGH BOTTOM AND TOP OF PRISM
00358 !       ARE NEGLECTED FOR MONOTONICITY
00359 !       SEE THE WHOLE TERM IN THE OPTION WITHOUT MONOTONICITY
00360 !
00361         XM(IELEM,03)=-NH1*SURFAC(IELEM)*XS03
00362         XM(IELEM,08)=-NH2*SURFAC(IELEM)*XS03
00363         XM(IELEM,12)=-NH3*SURFAC(IELEM)*XS03
00364 !
00365 !-----------------------------------------------------------------------
00366 !
00367       ENDDO
00368 !
00369 !-----------------------------------------------------------------------
00370 !
00371       ELSEIF(FORMUL(10:13).EQ.'1234') THEN
00372 !
00373 !-----------------------------------------------------------------------
00374 !
00375 !     VERSION WITHOUT MONOTONICITY AND ALL TERMS
00376 !
00377 !-----------------------------------------------------------------------
00378 !
00379       DO IELEM=1,NELEM
00380 !
00381         I1=IKLE(IELEM,1)
00382         I2=IKLE(IELEM,2)
00383         I3=IKLE(IELEM,3)
00384         I4=IKLE(IELEM,4)
00385         I5=IKLE(IELEM,5)
00386         I6=IKLE(IELEM,6)
00387 !       SUIVANT NPOU0, II4 SERA I4 OU I1, ETC.
00388         II4=I1+NPOU0
00389         II5=I2+NPOU0
00390         II6=I3+NPOU0
00391 !
00392         H1=Z(I4)-Z(I1)
00393         H2=Z(I5)-Z(I2)
00394         H3=Z(I6)-Z(I3)
00395 !
00396         IF((INCHYD.AND.MAX(Z(I1),Z(I2),Z(I3)).GT.
00397      &                 MIN(Z(I4),Z(I5),Z(I6)))    .OR.
00398      &      H1.LT.EPSI.OR.H2.LT.EPSI.OR.H3.LT.EPSI ) THEN
00399           NF1=0.D0
00400           NF2=0.D0
00401           NF3=0.D0
00402           NG1=0.D0
00403           NG2=0.D0
00404           NG3=0.D0
00405           NH1=0.D0
00406           NH2=0.D0
00407           NH3=0.D0
00408           AUX=0.D0
00409           NUXMOY=0.D0
00410           NUYMOY=0.D0
00411         ELSE
00412 !         SUIVANT LES CAS (II4=I1 OU I4, ETC.)
00413 !         VARIANTE AVEC VISCOSITE VERTICALE LINEAIRE (II4=I4)
00414 !         VARIANTE AVEC VISCOSITE VERTICALE P0 SUR LA VERTICALE (II4=I1)
00415           NF1=0.5D0*(F(I1)+F(II4))/H1
00416           NF2=0.5D0*(F(I2)+F(II5))/H2
00417           NF3=0.5D0*(F(I3)+F(II6))/H3
00418           NG1=0.5D0*(G(I1)+G(II4))/H1
00419           NG2=0.5D0*(G(I2)+G(II5))/H2
00420           NG3=0.5D0*(G(I3)+G(II6))/H3
00421           NH1=0.5D0*(H(I1)+H(II4))/H1
00422           NH2=0.5D0*(H(I2)+H(II5))/H2
00423           NH3=0.5D0*(H(I3)+H(II6))/H3
00424           AUX = XS24 / SURFAC(IELEM)
00425 !         TAKING INTO ACCOUNT AN AVERAGE DIFFUSION ON THE PRISM
00426           NUXMOY=(F(I1)+F(I2)+F(I3)+F(I4)+F(I5)+F(I6))/6.D0
00427           NUYMOY=(G(I1)+G(I2)+G(I3)+G(I4)+G(I5)+G(I6))/6.D0
00428         ENDIF
00429 !
00430 !       X2 = X(I2)-X(I1)
00431 !       X3 = X(I3)-X(I1)
00432 !       Y2 = Y(I2)-Y(I1)
00433 !       Y3 = Y(I3)-Y(I1)
00434         X2 = X(IELEM,2)
00435         X3 = X(IELEM,3)
00436         Y2 = Y(IELEM,2)
00437         Y3 = Y(IELEM,3)
00438 !
00439         X2X3  = X2*X3
00440         Y2Y3  = Y2*Y3
00441         X2AUX = X2X3-X2**2
00442         X3AUX = X3**2-X2X3
00443         Y2AUX = Y2Y3-Y2**2
00444         Y3AUX = Y3**2-Y2Y3
00445 !
00446 !       2D DIFFUSION, LOWER LEVEL (SEE MT02AA)
00447 !       LIKE IN 2D BUT MULTIPLIED BY DELTA(Z)/2
00448         SOMVX = ( F(I1)*H1+F(I2)*H2+F(I3)*H3 ) * AUX
00449         SOMVY = ( G(I1)*H1+G(I2)*H2+G(I3)*H3 ) * AUX
00450 !
00451 !       OFF-DIAGONAL TERMS FOR POINTS OF LOWER LEVEL
00452 !       WITH MONOTONICITY ENSURED
00453 !
00454 !       TERM 1
00455 !
00456         XM(IELEM, 1) = - SOMVY * X3AUX - SOMVX * Y3AUX
00457         XM(IELEM, 2) =   SOMVY * X2AUX + SOMVX * Y2AUX
00458         XM(IELEM, 6) = - SOMVY * X2X3  - SOMVX * Y2Y3
00459 !
00460 !       2D DIFFUSION, UPPER LEVEL
00461         SOMVX = ( F(I4)*H1+F(I5)*H2+F(I6)*H3 ) * AUX
00462         SOMVY = ( G(I4)*H1+G(I5)*H2+G(I6)*H3 ) * AUX
00463 !
00464 !       OFF-DIAGONAL TERMS FOR POINTS OF UPPER LEVEL
00465 !       WITH MONOTONICITY ENSURED
00466 !
00467 !       TERM 1
00468 !
00469         XM(IELEM,13) = - SOMVY * X3AUX - SOMVX * Y3AUX
00470         XM(IELEM,14) =   SOMVY * X2AUX + SOMVX * Y2AUX
00471         XM(IELEM,15) = - SOMVY * X2X3  - SOMVX * Y2Y3
00472 !
00473 !       AVERAGE OF NORMAL VECTOR TO PLANES (NOT NORMED)
00474 !       ONE CAN CHECK THAT WE GET -1 0 OR 0 -1 WITH Z=X OR Z=Y
00475 !       NPX=-DZ/DX   NPY=-DZ/DY
00476 !
00477         NPXL=-0.5D0*(Y2*(Z(I1)-Z(I3))+Y3*(Z(I2)-Z(I1)))
00478         NPXU=-0.5D0*(Y2*(Z(I4)-Z(I6))+Y3*(Z(I5)-Z(I4)))
00479         NPYL=-0.5D0*(X2*(Z(I3)-Z(I1))+X3*(Z(I1)-Z(I2)))
00480         NPYU=-0.5D0*(X2*(Z(I6)-Z(I4))+X3*(Z(I4)-Z(I5)))
00481         NPX=0.5D0*(NPXL+NPXU)/SURFAC(IELEM)
00482         NPY=0.5D0*(NPYL+NPYU)/SURFAC(IELEM)
00483         NPX2=NPX**2
00484         NPY2=NPY**2
00485 !
00486 !       TERMS 2 AND 3
00487 !
00488         NPX=NPX*NUXMOY
00489         NPY=NPY*NUYMOY
00490 !
00491 !       2D GRADIENT MATRIX  GX(3,3) AND GY(3,3)
00492 !       BUT GX(1,J)=GX(2,J)=GX(3,J)
00493 !       AND GY(1,J)=GY(2,J)=GY(3,J), HENCE ONLY GX(3) AND GY(3)
00494 !
00495         GX(2) =   Y3*XS06
00496         GX(3) = - Y2*XS06
00497         GX(1) = - GX(2) - GX(3)
00498 !
00499         GY(2) = - X3 * XS06
00500         GY(3) =   X2 * XS06
00501         GY(1) = - GY(2) - GY(3)
00502 !
00503 !       OFF-DIAGONAL TERMS, SOME (UP AND DOWN) ALREADY INITIALISED,
00504 !                           SOME (CROSSED TERMS) NOT
00505 !
00506 !       TERMS 2
00507 !
00508 !       TERM 1-2 (01)
00509         XM(IELEM,01)=XM(IELEM,01)+0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00510 !       TERM 1-3 (02)
00511         XM(IELEM,02)=XM(IELEM,02)+0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00512 !       TERM 1-4 (03)
00513         XM(IELEM,03)=            +0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00514 !       TERM 1-5 (04)
00515         XM(IELEM,04)=            +0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00516 !       TERM 1-6 (05)
00517         XM(IELEM,05)=            +0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00518 !       TERM 2-3 (06)
00519         XM(IELEM,06)=XM(IELEM,06)+0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00520 !       TERM 2-4 (07)
00521         XM(IELEM,07)=            +0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00522 !       TERM 2-5 (08)
00523         XM(IELEM,08)=            +0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00524 !       TERM 2-6 (09)
00525         XM(IELEM,09)=            +0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00526 !       TERM 3-4 (10)
00527         XM(IELEM,10)=            +0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00528 !       TERM 3-5 (11)
00529         XM(IELEM,11)=            +0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00530 !       TERM 3-6 (12)
00531         XM(IELEM,12)=            +0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00532 !       TERME 4-5 (13)
00533         XM(IELEM,13)=XM(IELEM,13)+0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00534 !       TERM 4-6 (14)
00535         XM(IELEM,14)=XM(IELEM,14)+0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00536 !       TERM 5-6 (15)
00537         XM(IELEM,15)=XM(IELEM,15)+0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00538 !
00539 !       TERMS 3
00540 !
00541 !       TERM 1-2 (01)
00542         XM(IELEM,01)=XM(IELEM,01)+0.5D0*( - ( NPX*GX(1)+NPY*GY(1)) )
00543 !       TERM 1-3 (02)
00544         XM(IELEM,02)=XM(IELEM,02)+0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00545 !       TERM 1-4 (03)
00546         XM(IELEM,03)=XM(IELEM,03)+0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00547 !       TERM 1-5 (04)
00548         XM(IELEM,04)=XM(IELEM,04)+0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00549 !       TERM 1-6 (05)
00550         XM(IELEM,05)=XM(IELEM,05)+0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00551 !       TERM 2-3 (06)
00552         XM(IELEM,06)=XM(IELEM,06)+0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00553 !       TERM 2-4 (07)
00554         XM(IELEM,07)=XM(IELEM,07)+0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00555 !       TERM 2-5 (08)
00556         XM(IELEM,08)=XM(IELEM,08)+0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00557 !       TERM 2-6 (09)
00558         XM(IELEM,09)=XM(IELEM,09)+0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00559 !       TERM 3-4 (10)
00560         XM(IELEM,10)=XM(IELEM,10)+0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00561 !       TERM 3-5 (11)
00562         XM(IELEM,11)=XM(IELEM,11)+0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00563 !       TERM 3-6 (12)
00564         XM(IELEM,12)=XM(IELEM,12)+0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00565 !       TERME 4-5 (13)
00566         XM(IELEM,13)=XM(IELEM,13)+0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00567 !       TERM 4-6 (14)
00568         XM(IELEM,14)=XM(IELEM,14)+0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00569 !       TERM 5-6 (15)
00570         XM(IELEM,15)=XM(IELEM,15)+0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00571 !
00572 !       TERM 4
00573 !
00574         XM(IELEM,03)=XM(IELEM,03)
00575      &              -(NPX2*NF1+NPY2*NG1+NH1)*SURFAC(IELEM)*XS03
00576         XM(IELEM,08)=XM(IELEM,08)
00577      &              -(NPX2*NF2+NPY2*NG2+NH2)*SURFAC(IELEM)*XS03
00578         XM(IELEM,12)=XM(IELEM,12)
00579      &              -(NPX2*NF3+NPY2*NG3+NH3)*SURFAC(IELEM)*XS03
00580 !
00581 !-----------------------------------------------------------------------
00582 !
00583       ENDDO
00584 !
00585 !-----------------------------------------------------------------------
00586 !
00587       ELSEIF(FORMUL(10:13).EQ.'1 3 ') THEN
00588 !
00589 !-----------------------------------------------------------------------
00590 !
00591 !     VERSION WITHOUT MONOTONICITY AND ONLY TERMS 1 AND 3
00592 !
00593 !-----------------------------------------------------------------------
00594 !
00595       DO IELEM=1,NELEM
00596 !
00597         I1=IKLE(IELEM,1)
00598         I2=IKLE(IELEM,2)
00599         I3=IKLE(IELEM,3)
00600         I4=IKLE(IELEM,4)
00601         I5=IKLE(IELEM,5)
00602         I6=IKLE(IELEM,6)
00603 !       SUIVANT NPOU0, II4 SERA I4 OU I1, ETC.
00604         II4=I1+NPOU0
00605         II5=I2+NPOU0
00606         II6=I3+NPOU0
00607 !
00608         H1=Z(I4)-Z(I1)
00609         H2=Z(I5)-Z(I2)
00610         H3=Z(I6)-Z(I3)
00611 !
00612         IF((INCHYD.AND.MAX(Z(I1),Z(I2),Z(I3)).GT.
00613      &                 MIN(Z(I4),Z(I5),Z(I6)))    .OR.
00614      &      H1.LT.EPSI.OR.H2.LT.EPSI.OR.H3.LT.EPSI ) THEN
00615           NF1=0.D0
00616           NF2=0.D0
00617           NF3=0.D0
00618           NG1=0.D0
00619           NG2=0.D0
00620           NG3=0.D0
00621           NH1=0.D0
00622           NH2=0.D0
00623           NH3=0.D0
00624           AUX=0.D0
00625           NUXMOY=0.D0
00626           NUYMOY=0.D0
00627         ELSE
00628 !         SUIVANT LES CAS (II4=I1 OU I4, ETC.)
00629 !         VARIANTE AVEC VISCOSITE VERTICALE LINEAIRE (II4=I4)
00630 !         VARIANTE AVEC VISCOSITE VERTICALE P0 SUR LA VERTICALE (II4=I1)
00631           NF1=0.5D0*(F(I1)+F(II4))/H1
00632           NF2=0.5D0*(F(I2)+F(II5))/H2
00633           NF3=0.5D0*(F(I3)+F(II6))/H3
00634           NG1=0.5D0*(G(I1)+G(II4))/H1
00635           NG2=0.5D0*(G(I2)+G(II5))/H2
00636           NG3=0.5D0*(G(I3)+G(II6))/H3
00637           NH1=0.5D0*(H(I1)+H(II4))/H1
00638           NH2=0.5D0*(H(I2)+H(II5))/H2
00639           NH3=0.5D0*(H(I3)+H(II6))/H3
00640           AUX = XS24 / SURFAC(IELEM)
00641 !         TAKING INTO ACCOUNT AN AVERAGE DIFFUSION ON THE PRISM
00642           NUXMOY=(F(I1)+F(I2)+F(I3)+F(I4)+F(I5)+F(I6))/6.D0
00643           NUYMOY=(G(I1)+G(I2)+G(I3)+G(I4)+G(I5)+G(I6))/6.D0
00644         ENDIF
00645 !
00646 !       X2 = X(I2)-X(I1)
00647 !       X3 = X(I3)-X(I1)
00648 !       Y2 = Y(I2)-Y(I1)
00649 !       Y3 = Y(I3)-Y(I1)
00650         X2 = X(IELEM,2)
00651         X3 = X(IELEM,3)
00652         Y2 = Y(IELEM,2)
00653         Y3 = Y(IELEM,3)
00654 !
00655         X2X3  = X2*X3
00656         Y2Y3  = Y2*Y3
00657         X2AUX = X2X3-X2**2
00658         X3AUX = X3**2-X2X3
00659         Y2AUX = Y2Y3-Y2**2
00660         Y3AUX = Y3**2-Y2Y3
00661 !
00662 !       2D DIFFUSION, LOWER LEVEL (SEE MT02AA)
00663 !       LIKE IN 2D BUT MULTIPLIED BY DELTA(Z)/2
00664         SOMVX = ( F(I1)*H1+F(I2)*H2+F(I3)*H3 ) * AUX
00665         SOMVY = ( G(I1)*H1+G(I2)*H2+G(I3)*H3 ) * AUX
00666 !
00667 !       OFF-DIAGONAL TERMS FOR POINTS OF LOWER LEVEL
00668 !       WITH MONOTONICITY ENSURED
00669 !
00670 !       TERM 1
00671 !
00672         XM(IELEM, 1) = - SOMVY * X3AUX - SOMVX * Y3AUX
00673         XM(IELEM, 2) =   SOMVY * X2AUX + SOMVX * Y2AUX
00674         XM(IELEM, 6) = - SOMVY * X2X3  - SOMVX * Y2Y3
00675         XM(IELEM,16) = XM(IELEM, 1)
00676         XM(IELEM,17) = XM(IELEM, 2)
00677         XM(IELEM,21) = XM(IELEM, 6)
00678 !
00679 !       2D DIFFUSION, UPPER LEVEL
00680         SOMVX = ( F(I4)*H1+F(I5)*H2+F(I6)*H3 ) * AUX
00681         SOMVY = ( G(I4)*H1+G(I5)*H2+G(I6)*H3 ) * AUX
00682 !
00683 !       OFF-DIAGONAL TERMS FOR POINTS OF UPPER LEVEL
00684 !       WITH MONOTONICITY ENSURED
00685 !
00686 !       TERM 1
00687 !
00688         XM(IELEM,13) = - SOMVY * X3AUX - SOMVX * Y3AUX
00689         XM(IELEM,14) =   SOMVY * X2AUX + SOMVX * Y2AUX
00690         XM(IELEM,15) = - SOMVY * X2X3  - SOMVX * Y2Y3
00691         XM(IELEM,28) = XM(IELEM,13)
00692         XM(IELEM,29) = XM(IELEM,14)
00693         XM(IELEM,30) = XM(IELEM,15)
00694 !
00695 !       AVERAGE OF NORMAL VECTOR TO PLANES (NOT NORMED)
00696 !       ONE CAN CHECK THAT WE GET -1 0 OR 0 -1 WITH Z=X OR Z=Y
00697 !       NPX=-DZ/DX   NPY=-DZ/DY
00698 !
00699         NPXL=-0.5D0*(Y2*(Z(I1)-Z(I3))+Y3*(Z(I2)-Z(I1)))
00700         NPXU=-0.5D0*(Y2*(Z(I4)-Z(I6))+Y3*(Z(I5)-Z(I4)))
00701         NPYL=-0.5D0*(X2*(Z(I3)-Z(I1))+X3*(Z(I1)-Z(I2)))
00702         NPYU=-0.5D0*(X2*(Z(I6)-Z(I4))+X3*(Z(I4)-Z(I5)))
00703         NPX=0.5D0*(NPXL+NPXU)/SURFAC(IELEM)
00704         NPY=0.5D0*(NPYL+NPYU)/SURFAC(IELEM)
00705         NPX2=NPX**2
00706         NPY2=NPY**2
00707 !
00708 !       TERM 3
00709 !
00710         NPX=NPX*NUXMOY
00711         NPY=NPY*NUYMOY
00712 !
00713 !       2D GRADIENT MATRIX  GX(3,3) AND GY(3,3)
00714 !       BUT GX(1,J)=GX(2,J)=GX(3,J)
00715 !       AND GY(1,J)=GY(2,J)=GY(3,J), HENCE ONLY GX(3) AND GY(3)
00716 !
00717         GX(2) =   Y3*XS06
00718         GX(3) = - Y2*XS06
00719         GX(1) = - GX(2) - GX(3)
00720 !
00721         GY(2) = - X3 * XS06
00722         GY(3) =   X2 * XS06
00723         GY(1) = - GY(2) - GY(3)
00724 !
00725 !       OFF-DIAGONAL TERMS, SOME (UP AND DOWN) ALREADY INITIALISED,
00726 !                           SOME (CROSSED TERMS) NOT
00727 !
00728 !       TERMS 3
00729 !
00730 !       TERM 1-2 (01)
00731         XM(IELEM,01)=XM(IELEM,01)+0.5D0*( - ( NPX*GX(1)+NPY*GY(1)) )
00732 !       TERM 1-3 (02)
00733         XM(IELEM,02)=XM(IELEM,02)+0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00734 !       TERM 1-4 (03)
00735         XM(IELEM,03)=            +0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00736 !       TERM 1-5 (04)
00737         XM(IELEM,04)=            +0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00738 !       TERM 1-6 (05)
00739         XM(IELEM,05)=            +0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00740 !       TERM 2-3 (06)
00741         XM(IELEM,06)=XM(IELEM,06)+0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00742 !       TERM 2-4 (07)
00743         XM(IELEM,07)=            +0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00744 !       TERM 2-5 (08)
00745         XM(IELEM,08)=            +0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00746 !       TERM 2-6 (09)
00747         XM(IELEM,09)=            +0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00748 !       TERM 3-4 (10)
00749         XM(IELEM,10)=            +0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00750 !       TERM 3-5 (11)
00751         XM(IELEM,11)=            +0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00752 !       TERM 3-6 (12)
00753         XM(IELEM,12)=            +0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00754 !       TERME 4-5 (13)
00755         XM(IELEM,13)=XM(IELEM,13)+0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00756 !       TERM 4-6 (14)
00757         XM(IELEM,14)=XM(IELEM,14)+0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00758 !       TERM 5-6 (15)
00759         XM(IELEM,15)=XM(IELEM,15)+0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00760 !
00761 !       TERM 2-1 (16)
00762         XM(IELEM,16)=XM(IELEM,16)+0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00763 !       TERM 3-1 (17)
00764         XM(IELEM,17)=XM(IELEM,17)+0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00765 !       TERM 4-1 (18)
00766         XM(IELEM,18)=             0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00767 !       TERM 5-1 (19)
00768         XM(IELEM,19)=             0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00769 !       TERM 6-1 (20)
00770         XM(IELEM,20)=             0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00771 !       TERM 3-2 (21)
00772         XM(IELEM,21)=XM(IELEM,21)+0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00773 !       TERM 4-2 (22)
00774         XM(IELEM,22)=             0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00775 !       TERM 5-2 (23)
00776         XM(IELEM,23)=             0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00777 !       TERM 6-2 (24)
00778         XM(IELEM,24)=             0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00779 !       TERM 4-3 (25)
00780         XM(IELEM,25)=             0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00781 !       TERM 5-3 (26)
00782         XM(IELEM,26)=             0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00783 !       TERM 5-4 (28)
00784         XM(IELEM,28)=XM(IELEM,28)+0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00785 !       TERM 6-3 (27)
00786         XM(IELEM,27)=             0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00787 !       TERM 6-4 (29)
00788         XM(IELEM,29)=XM(IELEM,29)+0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00789 !       TERM 6-5 (30)
00790         XM(IELEM,30)=XM(IELEM,30)+0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00791 !
00792 !-----------------------------------------------------------------------
00793 !
00794       ENDDO
00795 !
00796 !-----------------------------------------------------------------------
00797 !
00798       ELSEIF(FORMUL(10:13).EQ.' 2 4') THEN
00799 !
00800 !-----------------------------------------------------------------------
00801 !
00802 !     VERSION WITHOUT MONOTONICITY AND ONLY TERMS 2 AND 4
00803 !
00804 !-----------------------------------------------------------------------
00805 !
00806       DO IELEM=1,NELEM
00807 !
00808         I1=IKLE(IELEM,1)
00809         I2=IKLE(IELEM,2)
00810         I3=IKLE(IELEM,3)
00811         I4=IKLE(IELEM,4)
00812         I5=IKLE(IELEM,5)
00813         I6=IKLE(IELEM,6)
00814 !       SUIVANT NPOU0, II4 SERA I4 OU I1, ETC.
00815         II4=I1+NPOU0
00816         II5=I2+NPOU0
00817         II6=I3+NPOU0
00818 !
00819         H1=Z(I4)-Z(I1)
00820         H2=Z(I5)-Z(I2)
00821         H3=Z(I6)-Z(I3)
00822 !
00823         IF((INCHYD.AND.MAX(Z(I1),Z(I2),Z(I3)).GT.
00824      &                 MIN(Z(I4),Z(I5),Z(I6)))    .OR.
00825      &      H1.LT.EPSI.OR.H2.LT.EPSI.OR.H3.LT.EPSI ) THEN
00826           NF1=0.D0
00827           NF2=0.D0
00828           NF3=0.D0
00829           NG1=0.D0
00830           NG2=0.D0
00831           NG3=0.D0
00832           NH1=0.D0
00833           NH2=0.D0
00834           NH3=0.D0
00835           NUXMOY=0.D0
00836           NUYMOY=0.D0
00837         ELSE
00838 !         SUIVANT LES CAS (II4=I1 OU I4, ETC.)
00839 !         VARIANTE AVEC VISCOSITE VERTICALE LINEAIRE (II4=I4)
00840 !         VARIANTE AVEC VISCOSITE VERTICALE P0 SUR LA VERTICALE (II4=I1)
00841           NF1=0.5D0*(F(I1)+F(II4))/H1
00842           NF2=0.5D0*(F(I2)+F(II5))/H2
00843           NF3=0.5D0*(F(I3)+F(II6))/H3
00844           NG1=0.5D0*(G(I1)+G(II4))/H1
00845           NG2=0.5D0*(G(I2)+G(II5))/H2
00846           NG3=0.5D0*(G(I3)+G(II6))/H3
00847           NH1=0.5D0*(H(I1)+H(II4))/H1
00848           NH2=0.5D0*(H(I2)+H(II5))/H2
00849           NH3=0.5D0*(H(I3)+H(II6))/H3
00850 !         TAKING INTO ACCOUNT AN AVERAGE DIFFUSION ON THE PRISM
00851           NUXMOY=(F(I1)+F(I2)+F(I3)+F(I4)+F(I5)+F(I6))/6.D0
00852           NUYMOY=(G(I1)+G(I2)+G(I3)+G(I4)+G(I5)+G(I6))/6.D0
00853         ENDIF
00854 !
00855 !       X2 = X(I2)-X(I1)
00856 !       X3 = X(I3)-X(I1)
00857 !       Y2 = Y(I2)-Y(I1)
00858 !       Y3 = Y(I3)-Y(I1)
00859         X2 = X(IELEM,2)
00860         X3 = X(IELEM,3)
00861         Y2 = Y(IELEM,2)
00862         Y3 = Y(IELEM,3)
00863 !
00864 !       AVERAGE OF NORMAL VECTOR TO PLANES (NOT NORMED)
00865 !       ONE CAN CHECK THAT WE GET -1 0 OR 0 -1 WITH Z=X OR Z=Y
00866 !       NPX=-DZ/DX   NPY=-DZ/DY
00867 !
00868         NPXL=-0.5D0*(Y2*(Z(I1)-Z(I3))+Y3*(Z(I2)-Z(I1)))
00869         NPXU=-0.5D0*(Y2*(Z(I4)-Z(I6))+Y3*(Z(I5)-Z(I4)))
00870         NPYL=-0.5D0*(X2*(Z(I3)-Z(I1))+X3*(Z(I1)-Z(I2)))
00871         NPYU=-0.5D0*(X2*(Z(I6)-Z(I4))+X3*(Z(I4)-Z(I5)))
00872         NPX=0.5D0*(NPXL+NPXU)/SURFAC(IELEM)
00873         NPY=0.5D0*(NPYL+NPYU)/SURFAC(IELEM)
00874         NPX2=NPX**2
00875         NPY2=NPY**2
00876 !
00877 !       TERMS 2
00878 !
00879         NPX=NPX*NUXMOY
00880         NPY=NPY*NUYMOY
00881 !
00882 !       2D GRADIENT MATRIX  GX(3,3) AND GY(3,3)
00883 !       BUT GX(1,J)=GX(2,J)=GX(3,J)
00884 !       AND GY(1,J)=GY(2,J)=GY(3,J), HENCE ONLY GX(3) AND GY(3)
00885 !
00886         GX(2) =   Y3*XS06
00887         GX(3) = - Y2*XS06
00888         GX(1) = - GX(2) - GX(3)
00889 !
00890         GY(2) = - X3 * XS06
00891         GY(3) =   X2 * XS06
00892         GY(1) = - GY(2) - GY(3)
00893 !
00894 !       OFF-DIAGONAL TERMS
00895 !
00896 !       TERMS 2
00897 !
00898 !       TERM 1-2 (01)
00899         XM(IELEM,01)=0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00900 !       TERM 1-3 (02)
00901         XM(IELEM,02)=0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00902 !       TERM 1-4 (03)
00903         XM(IELEM,03)=0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00904 !       TERM 1-5 (04)
00905         XM(IELEM,04)=0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00906 !       TERM 1-6 (05)
00907         XM(IELEM,05)=0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00908 !       TERM 2-3 (06)
00909         XM(IELEM,06)=0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00910 !       TERM 2-4 (07)
00911         XM(IELEM,07)=0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00912 !       TERM 2-5 (08)
00913         XM(IELEM,08)=0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00914 !       TERM 2-6 (09)
00915         XM(IELEM,09)=0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00916 !       TERM 3-4 (10)
00917         XM(IELEM,10)=0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00918 !       TERM 3-5 (11)
00919         XM(IELEM,11)=0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00920 !       TERM 3-6 (12)
00921         XM(IELEM,12)=0.5D0*( - ( NPX*GX(3)+NPY*GY(3) ) )
00922 !       TERM 4-5 (13)
00923         XM(IELEM,13)=0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00924 !       TERM 4-6 (14)
00925         XM(IELEM,14)=0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00926 !       TERM 5-6 (15)
00927         XM(IELEM,15)=0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00928 !
00929 !       TERM 2-1 (16)
00930         XM(IELEM,16)=0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00931 !       TERM 3-1 (17)
00932         XM(IELEM,17)=0.5D0*( - ( NPX*GX(1)+NPY*GY(1) ) )
00933 !       TERM 4-1 (18)
00934         XM(IELEM,18)=0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00935 !       TERM 5-1 (19)
00936         XM(IELEM,19)=0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00937 !       TERM 6-1 (20)
00938         XM(IELEM,20)=0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00939 !       TERM 3-2 (21)
00940         XM(IELEM,21)=0.5D0*( - ( NPX*GX(2)+NPY*GY(2) ) )
00941 !       TERM 4-2 (22)
00942         XM(IELEM,22)=0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00943 !       TERM 5-2 (23)
00944         XM(IELEM,23)=0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00945 !       TERM 6-2 (24)
00946         XM(IELEM,24)=0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00947 !       TERM 4-3 (25)
00948         XM(IELEM,25)=0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00949 !       TERM 5-3 (26)
00950         XM(IELEM,26)=0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00951 !       TERM 6-3 (27)
00952         XM(IELEM,27)=0.5D0*( + ( NPX*GX(3)+NPY*GY(3) ) )
00953 !       TERM 5-4 (28)
00954         XM(IELEM,28)=0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00955 !       TERM 6-4 (29)
00956         XM(IELEM,29)=0.5D0*( + ( NPX*GX(1)+NPY*GY(1) ) )
00957 !       TERM 6-5 (30)
00958         XM(IELEM,30)=0.5D0*( + ( NPX*GX(2)+NPY*GY(2) ) )
00959 !
00960 !       TERM 4
00961 !
00962         XM(IELEM,03)=XM(IELEM,03)
00963      &              -(NPX2*NF1+NPY2*NG1+NH1)*SURFAC(IELEM)*XS03
00964         XM(IELEM,08)=XM(IELEM,08)
00965      &              -(NPX2*NF2+NPY2*NG2+NH2)*SURFAC(IELEM)*XS03
00966         XM(IELEM,12)=XM(IELEM,12)
00967      &              -(NPX2*NF3+NPY2*NG3+NH3)*SURFAC(IELEM)*XS03
00968         XM(IELEM,18)=XM(IELEM,18)
00969      &              -(NPX2*NF1+NPY2*NG1+NH1)*SURFAC(IELEM)*XS03
00970         XM(IELEM,23)=XM(IELEM,23)
00971      &              -(NPX2*NF2+NPY2*NG2+NH2)*SURFAC(IELEM)*XS03
00972         XM(IELEM,27)=XM(IELEM,27)
00973      &              -(NPX2*NF3+NPY2*NG3+NH3)*SURFAC(IELEM)*XS03
00974 !
00975 !-----------------------------------------------------------------------
00976 !
00977       ENDDO
00978 !
00979 !-----------------------------------------------------------------------
00980 !
00981       ELSE
00982         IF(LNG.EQ.1) WRITE(LU,*) 'MT02PP_STAR (BIEF): FORMULE INCONNUE'
00983         IF(LNG.EQ.2) WRITE(LU,*) 'MT02PP_STAR (BIEF): UNKNOWN FORMULA'
00984         WRITE(LU,*) FORMUL
00985         CALL PLANTE(1)
00986         STOP
00987       ENDIF
00988 !
00989 !-----------------------------------------------------------------------
00990 !
00991 !     DIAGONAL TERMS OBTAINED BY THE FACT THAT THE SUM OF TERMS IN A
00992 !     LINE IS 0
00993 !
00994 !-----------------------------------------------------------------------
00995 !
00996 !     IN SYMMETRIC MODE
00997 !
00998       IF(FORMUL(11:13).EQ.'MON'.OR.FORMUL(10:13).EQ.'1234') THEN
00999 !
01000       DO IELEM=1,NELEM
01001 !
01002         T(IELEM,1)= -XM(IELEM,01)
01003      &              -XM(IELEM,02)
01004      &              -XM(IELEM,03)
01005      &              -XM(IELEM,04)
01006      &              -XM(IELEM,05)
01007         T(IELEM,2)= -XM(IELEM,01)
01008      &              -XM(IELEM,06)
01009      &              -XM(IELEM,07)
01010      &              -XM(IELEM,08)
01011      &              -XM(IELEM,09)
01012         T(IELEM,3)= -XM(IELEM,02)
01013      &              -XM(IELEM,06)
01014      &              -XM(IELEM,10)
01015      &              -XM(IELEM,11)
01016      &              -XM(IELEM,12)
01017         T(IELEM,4)= -XM(IELEM,03)
01018      &              -XM(IELEM,07)
01019      &              -XM(IELEM,10)
01020      &              -XM(IELEM,13)
01021      &              -XM(IELEM,14)
01022         T(IELEM,5)= -XM(IELEM,04)
01023      &              -XM(IELEM,08)
01024      &              -XM(IELEM,11)
01025      &              -XM(IELEM,13)
01026      &              -XM(IELEM,15)
01027         T(IELEM,6)= -XM(IELEM,05)
01028      &              -XM(IELEM,09)
01029      &              -XM(IELEM,12)
01030      &              -XM(IELEM,14)
01031      &              -XM(IELEM,15)
01032 !
01033       ENDDO
01034 !
01035 !     IN NON SYMMETRIC MODE
01036 !
01037       ELSEIF(FORMUL(10:13).EQ.'1 3 '.OR.FORMUL(10:13).EQ.' 2 4') THEN
01038 !
01039       DO IELEM=1,NELEM
01040 !
01041         T(IELEM,1)= -XM(IELEM,01)
01042      &              -XM(IELEM,02)
01043      &              -XM(IELEM,03)
01044      &              -XM(IELEM,04)
01045      &              -XM(IELEM,05)
01046         T(IELEM,2)= -XM(IELEM,16)
01047      &              -XM(IELEM,06)
01048      &              -XM(IELEM,07)
01049      &              -XM(IELEM,08)
01050      &              -XM(IELEM,09)
01051         T(IELEM,3)= -XM(IELEM,17)
01052      &              -XM(IELEM,21)
01053      &              -XM(IELEM,10)
01054      &              -XM(IELEM,11)
01055      &              -XM(IELEM,12)
01056         T(IELEM,4)= -XM(IELEM,18)
01057      &              -XM(IELEM,22)
01058      &              -XM(IELEM,25)
01059      &              -XM(IELEM,13)
01060      &              -XM(IELEM,14)
01061         T(IELEM,5)= -XM(IELEM,19)
01062      &              -XM(IELEM,23)
01063      &              -XM(IELEM,26)
01064      &              -XM(IELEM,28)
01065      &              -XM(IELEM,15)
01066         T(IELEM,6)= -XM(IELEM,20)
01067      &              -XM(IELEM,24)
01068      &              -XM(IELEM,27)
01069      &              -XM(IELEM,29)
01070      &              -XM(IELEM,30)
01071 !
01072       ENDDO
01073 !
01074       ELSE
01075         IF(LNG.EQ.1) WRITE(LU,*) 'MT02PP_STAR (BIEF): FORMULE INCONNUE'
01076         IF(LNG.EQ.2) WRITE(LU,*) 'MT02PP_STAR (BIEF): UNKNOWN FORMULA'
01077         WRITE(LU,*) FORMUL
01078         CALL PLANTE(1)
01079         STOP
01080       ENDIF
01081 !
01082 !-----------------------------------------------------------------------
01083 !
01084       RETURN
01085       END

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