mt05pp.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\mt05pp.f
00002 !
00087                      SUBROUTINE MT05PP
00088 !                    *****************
00089 !
00090      &(T,XM,XMUL,SU,SV,SW,U,V,W,SF,SG,SH,F,G,H,
00091      & X,Y,Z,IKLE,NELEM,NELMAX,SURFAC,SIGMAG,SPECAD,NPLAN)
00092 !
00093 !***********************************************************************
00094 ! BIEF   V6P3                                   21/08/2010
00095 !***********************************************************************
00096 !
00097 !
00098 !
00099 !
00100 !
00101 !
00102 !
00103 !
00104 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00105 !| F              |-->| FUNCTION USED IN THE FORMULA
00106 !| G              |-->| FUNCTION USED IN THE FORMULA
00107 !| FORMUL         |-->| FORMULA DESCRIBING THE RESULTING MATRIX
00108 !| H              |-->| FUNCTION USED IN THE FORMULA
00109 !| IKLE           |-->| CONNECTIVITY TABLE.
00110 !| NELEM          |-->| NUMBER OF ELEMENTS
00111 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00112 !| SF             |-->| STRUCTURE OF FUNCTIONS F
00113 !| SG             |-->| STRUCTURE OF FUNCTIONS G
00114 !| SH             |-->| STRUCTURE OF FUNCTIONS H
00115 !| SIGMAG         |-->| IF YES, GENERALISED SIGMA TRANSFORMATION
00116 !| SPECAD         |-->| IF YES, SPECIAL ADVECTION FIELD (SEE IMPLEMENTATION)
00117 !| SU             |-->| STRUCTURE OF FUNCTIONS U
00118 !| SV             |-->| STRUCTURE OF FUNCTIONS V
00119 !| SW             |-->| STRUCTURE OF FUNCTIONS W
00120 !| SURFAC         |-->| AREA OF 2D ELEMENTS
00121 !| T              |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
00122 !| U              |-->| FUNCTION USED IN THE FORMULA
00123 !| V              |-->| FUNCTION USED IN THE FORMULA
00124 !| W              |-->| FUNCTION USED IN THE FORMULA
00125 !| X              |-->| ABSCISSAE OF POINTS, PER ELEMENT
00126 !| Y              |-->| ORDINATES OF POINTS, PER ELEMENT
00127 !| Z              |-->| ELEVATIONS OF POINTS, PER POINTS !!!!
00128 !| XM             |<->| OFF-DIAGONAL TERMS
00129 !| XMUL           |-->| COEFFICIENT FOR MULTIPLICATION
00130 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00131 !
00132       USE BIEF, EX_MT05PP => MT05PP
00133 !
00134       IMPLICIT NONE
00135       INTEGER LNG,LU
00136       COMMON/INFO/LNG,LU
00137 !
00138 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00139 !
00140       INTEGER, INTENT(IN)             :: NELEM,NELMAX,NPLAN
00141       INTEGER, INTENT(IN)             :: IKLE(NELMAX,6)
00142 !
00143       DOUBLE PRECISION, INTENT(INOUT) :: T(NELMAX,6),XM(NELMAX,30)
00144 !
00145       DOUBLE PRECISION, INTENT(IN)    :: XMUL
00146       DOUBLE PRECISION, INTENT(IN)    :: U(*),V(*),W(*),SURFAC(*)
00147       DOUBLE PRECISION, INTENT(IN)    :: F(*),G(*),H(*)
00148 !
00149 !     STRUCTURES OF U, V, W
00150 !
00151       TYPE(BIEF_OBJ), INTENT(IN)      :: SU,SV,SW,SF,SG,SH
00152 !
00153       DOUBLE PRECISION, INTENT(IN)    :: X(NELMAX,6),Y(NELMAX,6),Z(*)
00154       LOGICAL, INTENT(IN)             :: SIGMAG,SPECAD
00155 !
00156 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00157 !
00158 !     DECLARATIONS SPECIFIC TO THIS SUBROUTINE
00159 !
00160       DOUBLE PRECISION INT1,X2,X3,Y2,Y3,DEN,DET,G2,G3,GRADGX,GRADGY
00161       DOUBLE PRECISION PZ1,PX1,PX2,PX3,PY1,PY2,PY3
00162       DOUBLE PRECISION U1,U2,U3,U4,U5,U6,V1,V2,V3,V4,V5,V6
00163       DOUBLE PRECISION Q1,Q2,Q3,H1,H2,H3,HT
00164       DOUBLE PRECISION SUS,SUI,SHUI,SHUS
00165       DOUBLE PRECISION SHUI1,SHUS1,SHUI2,SHUS2,SHUI3,SHUS3
00166       DOUBLE PRECISION HU1S3I,HU1I3S,HU1SI,HU2S3I,HU2I3S,HU2SI
00167       DOUBLE PRECISION HU3S3I,HU3I3S,HU3SI
00168       DOUBLE PRECISION SVS,SVI,SHVI,SHVS
00169       DOUBLE PRECISION SHVI1,SHVS1,SHVI2,SHVS2,SHVI3,SHVS3
00170       DOUBLE PRECISION HV1S3I,HV1I3S,HV1SI,HV2S3I,HV2I3S,HV2SI
00171       DOUBLE PRECISION HV3S3I,HV3I3S,HV3SI,HH1,HH2,HH3,DZ1,DZ2,DZ3
00172 !
00173       INTEGER I1,I2,I3,I4,I5,I6,IELEM,IELEM2,IELEM3,NELEM2
00174 !
00175 !**********************************************************************
00176 !
00177       IF(SU%ELM.NE.41) THEN
00178         IF (LNG.EQ.1) WRITE(LU,1000) SU%ELM
00179         IF (LNG.EQ.2) WRITE(LU,1001) SU%ELM
00180 1000    FORMAT(1X,'MT05PP (BIEF) : TYPE DE U NON PREVU : ',I6)
00181 1001    FORMAT(1X,'MT05PP (BIEF) : TYPE OF U NOT IMPLEMENTED: ',I6)
00182         CALL PLANTE(1)
00183         STOP
00184       ENDIF
00185       IF(SV%ELM.NE.41) THEN
00186         IF (LNG.EQ.1) WRITE(LU,2000) SV%ELM
00187         IF (LNG.EQ.2) WRITE(LU,2001) SV%ELM
00188 2000    FORMAT(1X,'MT05PP (BIEF) : TYPE DE V NON PREVU : ',I6)
00189 2001    FORMAT(1X,'MT05PP (BIEF) : TYPE OF V NOT IMPLEMENTED: ',I6)
00190         CALL PLANTE(1)
00191         STOP
00192       ENDIF
00193 !     IN FACT W DEFINED PER LAYER BUT DECLARED 41
00194       IF(SW%ELM.NE.41) THEN
00195         IF (LNG.EQ.1) WRITE(LU,3000) SW%ELM
00196         IF (LNG.EQ.2) WRITE(LU,3001) SW%ELM
00197 3000    FORMAT(1X,'MT05PP (BIEF) : TYPE DE W NON PREVU : ',I6)
00198 3001    FORMAT(1X,'MT05PP (BIEF) : TYPE OF W NOT IMPLEMENTED: ',I6)
00199         CALL PLANTE(1)
00200         STOP
00201       ENDIF
00202 !
00203 !-----------------------------------------------------------------------
00204 !
00205 !     WITH NORMAL SIGMA TRANSFORMATION
00206 !
00207       IF(.NOT.SIGMAG) THEN
00208 !
00209 !     LOOP ON THE ELEMENTS
00210 !
00211 !     STANDARD CASE
00212 !
00213       IF(.NOT.SPECAD) THEN
00214 !
00215       DO IELEM=1,NELEM
00216 !
00217         I1 = IKLE(IELEM,1)
00218         I2 = IKLE(IELEM,2)
00219         I3 = IKLE(IELEM,3)
00220         I4 = IKLE(IELEM,4)
00221         I5 = IKLE(IELEM,5)
00222         I6 = IKLE(IELEM,6)
00223 !
00224 !       X2  =  X(I2) - X(I1)
00225 !       X3  =  X(I3) - X(I1)
00226 !       Y2  =  Y(I2) - Y(I1)
00227 !       Y3  =  Y(I3) - Y(I1)
00228         X2  =  X(IELEM,2)
00229         X3  =  X(IELEM,3)
00230         Y2  =  Y(IELEM,2)
00231         Y3  =  Y(IELEM,3)
00232 !
00233         U1  =  U(I1)
00234         U2  =  U(I2)
00235         U3  =  U(I3)
00236         U4  =  U(I4)
00237         U5  =  U(I5)
00238         U6  =  U(I6)
00239         V1  =  V(I1)
00240         V2  =  V(I2)
00241         V3  =  V(I3)
00242         V4  =  V(I4)
00243         V5  =  V(I5)
00244         V6  =  V(I6)
00245         Q1  =  W(I1)
00246         Q2  =  W(I2)
00247         Q3  =  W(I3)
00248 !
00249         H1  =  Z(I4)-Z(I1)
00250         H2  =  Z(I5)-Z(I2)
00251         H3  =  Z(I6)-Z(I3)
00252 !
00253 !    INTERMEDIATE COMPUTATIONS
00254 !
00255         DEN=XMUL/1440.D0
00256         PX2=Y3*DEN
00257         PX3=-Y2*DEN
00258         PX1=-PX2-PX3
00259         PY2=-X3*DEN
00260         PY3=X2*DEN
00261         PY1=-PY3-PY2
00262 !
00263         HT = H1+H2+H3
00264         SUS = U6+U5+U4
00265         SUI = U3+U2+U1
00266         SHUI = H3*U3+H2*U2+H1*U1
00267         SHUS = H3*U6+H2*U5+H1*U4
00268         SHUS1 = (HT+H1)*(SUS+U4)+SHUS+U4*H1
00269         SHUI1 = (HT+H1)*(SUI+U1)+SHUI+U1*H1
00270         SHUS2 = (HT+H2)*(SUS+U5)+SHUS+U5*H2
00271         SHUI2 = (HT+H2)*(SUI+U2)+SHUI+U2*H2
00272         SHUS3 = (HT+H3)*(SUS+U6)+SHUS+U6*H3
00273         SHUI3 = (HT+H3)*(SUI+U3)+SHUI+U3*H3
00274         HU1S3I = SHUS1+3*SHUI1
00275         HU1I3S = SHUI1+3*SHUS1
00276         HU1SI = SHUS1+SHUI1
00277         HU2S3I = SHUS2+3*SHUI2
00278         HU2I3S = SHUI2+3*SHUS2
00279         HU2SI = SHUS2+SHUI2
00280         HU3S3I = SHUS3+3*SHUI3
00281         HU3I3S = SHUI3+3*SHUS3
00282         HU3SI = SHUS3+SHUI3
00283 !
00284         SVS = V6+V5+V4
00285         SVI = V3+V2+V1
00286         SHVI = H3*V3+H2*V2+H1*V1
00287         SHVS = H3*V6+H2*V5+H1*V4
00288         SHVS1  = (HT+H1)*(SVS+V4)+SHVS+V4*H1
00289         SHVI1  = (HT+H1)*(SVI+V1)+SHVI+V1*H1
00290         SHVS2  = (HT+H2)*(SVS+V5)+SHVS+V5*H2
00291         SHVI2  = (HT+H2)*(SVI+V2)+SHVI+V2*H2
00292         SHVS3  = (HT+H3)*(SVS+V6)+SHVS+V6*H3
00293         SHVI3  = (HT+H3)*(SVI+V3)+SHVI+V3*H3
00294         HV1S3I = SHVS1+3*SHVI1
00295         HV1I3S = SHVI1+3*SHVS1
00296         HV1SI  = SHVS1+  SHVI1
00297         HV2S3I = SHVS2+3*SHVI2
00298         HV2I3S = SHVI2+3*SHVS2
00299         HV2SI  = SHVS2+  SHVI2
00300         HV3S3I = SHVS3+3*SHVI3
00301         HV3I3S = SHVI3+3*SHVS3
00302         HV3SI  = SHVS3+  SHVI3
00303 !
00304 !       OPTION WITH MASS-LUMPING ON VERTICAL VELOCITIES
00305 !       COMPATIBILITY WITH TRIDW2 (LUMPING TO COMPUTE DELTA(Z)*WSTAR)
00306 !
00307         PZ1=-DEN*(X2*Y3-Y2*X3)*30
00308 !
00309         T(IELEM,1)=PZ1*2*Q1
00310         XM(IELEM,3) =-T(IELEM,1)
00311         XM(IELEM,18)= T(IELEM,1)
00312         T(IELEM,4)  =-T(IELEM,1)
00313 !
00314         T(IELEM,2)=PZ1*2*Q2
00315         XM(IELEM,8) =-T(IELEM,2)
00316         XM(IELEM,23)= T(IELEM,2)
00317         T(IELEM,5)  =-T(IELEM,2)
00318 !
00319         T(IELEM,3)=PZ1*2*Q3
00320         XM(IELEM,12)=-T(IELEM,3)
00321         XM(IELEM,27)= T(IELEM,3)
00322         T(IELEM,6)  =-T(IELEM,3)
00323 !
00324         XM(IELEM,16)=PZ1*Q1
00325         XM(IELEM,7) =-XM(IELEM,16)
00326         XM(IELEM,17)= XM(IELEM,16)
00327         XM(IELEM,10)=-XM(IELEM,16)
00328         XM(IELEM,19)= XM(IELEM,16)
00329         XM(IELEM,28)=-XM(IELEM,16)
00330         XM(IELEM,20)= XM(IELEM,16)
00331         XM(IELEM,29)=-XM(IELEM,16)
00332 !
00333         XM(IELEM,1)=PZ1*Q2
00334         XM(IELEM,4) =-XM(IELEM,1)
00335         XM(IELEM,22)= XM(IELEM,1)
00336         XM(IELEM,13)=-XM(IELEM,1)
00337         XM(IELEM,21)= XM(IELEM,1)
00338         XM(IELEM,11)=-XM(IELEM,1)
00339         XM(IELEM,24)= XM(IELEM,1)
00340         XM(IELEM,30)=-XM(IELEM,1)
00341 !
00342         XM(IELEM,2)=PZ1*Q3
00343         XM(IELEM,5)=-XM(IELEM,2)
00344         XM(IELEM,25)=XM(IELEM,2)
00345         XM(IELEM,14)=-XM(IELEM,2)
00346         XM(IELEM,26)=XM(IELEM,2)
00347         XM(IELEM,15)=-XM(IELEM,2)
00348         XM(IELEM,6)=XM(IELEM,2)
00349         XM(IELEM,9)=-XM(IELEM,2)
00350 !
00351         T(IELEM,1)=T(IELEM,1)+PX1*HU1S3I+PY1*HV1S3I
00352         T(IELEM,2)=T(IELEM,2)+PX2*HU2S3I+PY2*HV2S3I
00353         T(IELEM,3)=T(IELEM,3)+PX3*HU3S3I+PY3*HV3S3I
00354         T(IELEM,4)=T(IELEM,4)+PX1*HU1I3S+PY1*HV1I3S
00355         T(IELEM,5)=T(IELEM,5)+PX2*HU2I3S+PY2*HV2I3S
00356         T(IELEM,6)=T(IELEM,6)+PX3*HU3I3S+PY3*HV3I3S
00357 !
00358         XM(IELEM,01)=XM(IELEM,01)+PX2*HU1S3I+PY2*HV1S3I
00359         XM(IELEM,02)=XM(IELEM,02)+PX3*HU1S3I+PY3*HV1S3I
00360         XM(IELEM,03)=XM(IELEM,03)+PX1*HU1SI +PY1*HV1SI
00361         XM(IELEM,04)=XM(IELEM,04)+PX2*HU1SI +PY2*HV1SI
00362         XM(IELEM,05)=XM(IELEM,05)+PX3*HU1SI +PY3*HV1SI
00363         XM(IELEM,06)=XM(IELEM,06)+PX3*HU2S3I+PY3*HV2S3I
00364         XM(IELEM,07)=XM(IELEM,07)+PX1*HU2SI +PY1*HV2SI
00365         XM(IELEM,08)=XM(IELEM,08)+PX2*HU2SI +PY2*HV2SI
00366         XM(IELEM,09)=XM(IELEM,09)+PX3*HU2SI +PY3*HV2SI
00367         XM(IELEM,10)=XM(IELEM,10)+PX1*HU3SI +PY1*HV3SI
00368         XM(IELEM,11)=XM(IELEM,11)+PX2*HU3SI +PY2*HV3SI
00369         XM(IELEM,12)=XM(IELEM,12)+PX3*HU3SI +PY3*HV3SI
00370         XM(IELEM,13)=XM(IELEM,13)+PX2*HU1I3S+PY2*HV1I3S
00371         XM(IELEM,14)=XM(IELEM,14)+PX3*HU1I3S+PY3*HV1I3S
00372         XM(IELEM,15)=XM(IELEM,15)+PX3*HU2I3S+PY3*HV2I3S
00373         XM(IELEM,16)=XM(IELEM,16)+PX1*HU2S3I+PY1*HV2S3I
00374         XM(IELEM,17)=XM(IELEM,17)+PX1*HU3S3I+PY1*HV3S3I
00375         XM(IELEM,18)=XM(IELEM,18)+PX1*HU1SI +PY1*HV1SI
00376         XM(IELEM,19)=XM(IELEM,19)+PX1*HU2SI +PY1*HV2SI
00377         XM(IELEM,20)=XM(IELEM,20)+PX1*HU3SI +PY1*HV3SI
00378         XM(IELEM,21)=XM(IELEM,21)+PX2*HU3S3I+PY2*HV3S3I
00379         XM(IELEM,22)=XM(IELEM,22)+PX2*HU1SI +PY2*HV1SI
00380         XM(IELEM,23)=XM(IELEM,23)+PX2*HU2SI +PY2*HV2SI
00381         XM(IELEM,24)=XM(IELEM,24)+PX2*HU3SI +PY2*HV3SI
00382         XM(IELEM,25)=XM(IELEM,25)+PX3*HU1SI +PY3*HV1SI
00383         XM(IELEM,26)=XM(IELEM,26)+PX3*HU2SI +PY3*HV2SI
00384         XM(IELEM,27)=XM(IELEM,27)+PX3*HU3SI +PY3*HV3SI
00385         XM(IELEM,28)=XM(IELEM,28)+PX1*HU2I3S+PY1*HV2I3S
00386         XM(IELEM,29)=XM(IELEM,29)+PX1*HU3I3S+PY1*HV3I3S
00387         XM(IELEM,30)=XM(IELEM,30)+PX2*HU3I3S+PY2*HV3I3S
00388 !
00389       ENDDO
00390 !
00391       ELSE
00392 !
00393 !     WITH SPECIFIC ADVECTING FIELD
00394 !
00395       NELEM2 = NELEM/(NPLAN-1)
00396 !
00397       DO IELEM=1,NELEM
00398 !
00399         I1 = IKLE(IELEM,1)
00400         I2 = IKLE(IELEM,2)
00401         I3 = IKLE(IELEM,3)
00402         I4 = IKLE(IELEM,4)
00403         I5 = IKLE(IELEM,5)
00404         I6 = IKLE(IELEM,6)
00405 !
00406 !       X2  =  X(I2) - X(I1)
00407 !       X3  =  X(I3) - X(I1)
00408 !       Y2  =  Y(I2) - Y(I1)
00409 !       Y3  =  Y(I3) - Y(I1)
00410         X2  =  X(IELEM,2)
00411         X3  =  X(IELEM,3)
00412         Y2  =  Y(IELEM,2)
00413         Y3  =  Y(IELEM,3)
00414 !
00415         DET= X2*Y3-X3*Y2
00416         IELEM2 = MOD(IELEM-1,NELEM2) + 1
00417 !       G IS PIECE-WISE LINEAR
00418 !       IT IS ZCONV IN TELEMAC-3D
00419         G2 = G(IELEM2+NELEM2)-G(IELEM2)
00420         G3 = G(IELEM2+2*NELEM2)-G(IELEM2)
00421         GRADGX=(G2*Y3-G3*Y2)/DET
00422         GRADGY=(X2*G3-X3*G2)/DET
00423 !
00424         U1=U(I1)+F(I1)*GRADGX
00425         U2=U(I2)+F(I2)*GRADGX
00426         U3=U(I3)+F(I3)*GRADGX
00427         U4=U(I4)+F(I4)*GRADGX
00428         U5=U(I5)+F(I5)*GRADGX
00429         U6=U(I6)+F(I6)*GRADGX
00430         V1=V(I1)+F(I1)*GRADGY
00431         V2=V(I2)+F(I2)*GRADGY
00432         V3=V(I3)+F(I3)*GRADGY
00433         V4=V(I4)+F(I4)*GRADGY
00434         V5=V(I5)+F(I5)*GRADGY
00435         V6=V(I6)+F(I6)*GRADGY
00436 !
00437         Q1  =  W(I1)
00438         Q2  =  W(I2)
00439         Q3  =  W(I3)
00440 !
00441         H1  =  Z(I4)-Z(I1)
00442         H2  =  Z(I5)-Z(I2)
00443         H3  =  Z(I6)-Z(I3)
00444 !
00445 !    INTERMEDIATE COMPUTATIONS
00446 !
00447         DEN=XMUL/1440.D0
00448         PX2=Y3*DEN
00449         PX3=-Y2*DEN
00450         PX1=-PX2-PX3
00451         PY2=-X3*DEN
00452         PY3=X2*DEN
00453         PY1=-PY3-PY2
00454 !
00455         HT = H1+H2+H3
00456         SUS = U6+U5+U4
00457         SUI = U3+U2+U1
00458         SHUI = H3*U3+H2*U2+H1*U1
00459         SHUS = H3*U6+H2*U5+H1*U4
00460         SHUS1 = (HT+H1)*(SUS+U4)+SHUS+U4*H1
00461         SHUI1 = (HT+H1)*(SUI+U1)+SHUI+U1*H1
00462         SHUS2 = (HT+H2)*(SUS+U5)+SHUS+U5*H2
00463         SHUI2 = (HT+H2)*(SUI+U2)+SHUI+U2*H2
00464         SHUS3 = (HT+H3)*(SUS+U6)+SHUS+U6*H3
00465         SHUI3 = (HT+H3)*(SUI+U3)+SHUI+U3*H3
00466         HU1S3I = SHUS1+3*SHUI1
00467         HU1I3S = SHUI1+3*SHUS1
00468         HU1SI = SHUS1+SHUI1
00469         HU2S3I = SHUS2+3*SHUI2
00470         HU2I3S = SHUI2+3*SHUS2
00471         HU2SI = SHUS2+SHUI2
00472         HU3S3I = SHUS3+3*SHUI3
00473         HU3I3S = SHUI3+3*SHUS3
00474         HU3SI = SHUS3+SHUI3
00475 !
00476         SVS = V6+V5+V4
00477         SVI = V3+V2+V1
00478         SHVI = H3*V3+H2*V2+H1*V1
00479         SHVS = H3*V6+H2*V5+H1*V4
00480         SHVS1  = (HT+H1)*(SVS+V4)+SHVS+V4*H1
00481         SHVI1  = (HT+H1)*(SVI+V1)+SHVI+V1*H1
00482         SHVS2  = (HT+H2)*(SVS+V5)+SHVS+V5*H2
00483         SHVI2  = (HT+H2)*(SVI+V2)+SHVI+V2*H2
00484         SHVS3  = (HT+H3)*(SVS+V6)+SHVS+V6*H3
00485         SHVI3  = (HT+H3)*(SVI+V3)+SHVI+V3*H3
00486         HV1S3I = SHVS1+3*SHVI1
00487         HV1I3S = SHVI1+3*SHVS1
00488         HV1SI  = SHVS1+  SHVI1
00489         HV2S3I = SHVS2+3*SHVI2
00490         HV2I3S = SHVI2+3*SHVS2
00491         HV2SI  = SHVS2+  SHVI2
00492         HV3S3I = SHVS3+3*SHVI3
00493         HV3I3S = SHVI3+3*SHVS3
00494         HV3SI  = SHVS3+  SHVI3
00495 !
00496 !       OPTION WITH MASS-LUMPING ON VERTICAL VELOCITIES
00497 !       COMPATIBILITY WITH TRIDW2 (LUMPING TO COMPUTE DELTA(Z)*WSTAR)
00498 !
00499         PZ1=-DEN*(X2*Y3-Y2*X3)*30
00500 !
00501         INT1=PZ1*2*Q1
00502         T(IELEM,1)=INT1
00503         XM(IELEM,3)=-INT1
00504         XM(IELEM,18)=INT1
00505         T(IELEM,4)=-INT1
00506 !
00507         INT1=PZ1*2*Q2
00508         T(IELEM,2)=INT1
00509         XM(IELEM,8)= -INT1
00510         XM(IELEM,23)=INT1
00511         T(IELEM,5)=-INT1
00512 !
00513         INT1=PZ1*2*Q3
00514         T(IELEM,3)=INT1
00515         XM(IELEM,12)=-INT1
00516         XM(IELEM,27)=INT1
00517         T(IELEM,6)=-INT1
00518 !
00519         INT1=PZ1*Q1
00520         XM(IELEM,16)=INT1
00521         XM(IELEM,7)=-INT1
00522         XM(IELEM,17)=INT1
00523         XM(IELEM,10)=-INT1
00524         XM(IELEM,19)=INT1
00525         XM(IELEM,28)=-INT1
00526         XM(IELEM,20)=INT1
00527         XM(IELEM,29)=-INT1
00528 !
00529         INT1=PZ1*Q2
00530         XM(IELEM,1)=INT1
00531         XM(IELEM,4)=-INT1
00532         XM(IELEM,22)=INT1
00533         XM(IELEM,13)=-INT1
00534         XM(IELEM,21)=INT1
00535         XM(IELEM,11)=-INT1
00536         XM(IELEM,24)=INT1
00537         XM(IELEM,30)=-INT1
00538 !
00539         INT1=PZ1*Q3
00540         XM(IELEM,2)=INT1
00541         XM(IELEM,5)=-INT1
00542         XM(IELEM,25)=INT1
00543         XM(IELEM,14)=-INT1
00544         XM(IELEM,26)=INT1
00545         XM(IELEM,15)=-INT1
00546         XM(IELEM,6)=INT1
00547         XM(IELEM,9)=-INT1
00548 !
00549         T(IELEM,1)=T(IELEM,1)+PX1*HU1S3I+PY1*HV1S3I
00550         T(IELEM,2)=T(IELEM,2)+PX2*HU2S3I+PY2*HV2S3I
00551         T(IELEM,3)=T(IELEM,3)+PX3*HU3S3I+PY3*HV3S3I
00552         T(IELEM,4)=T(IELEM,4)+PX1*HU1I3S+PY1*HV1I3S
00553         T(IELEM,5)=T(IELEM,5)+PX2*HU2I3S+PY2*HV2I3S
00554         T(IELEM,6)=T(IELEM,6)+PX3*HU3I3S+PY3*HV3I3S
00555 !
00556         XM(IELEM,01)=XM(IELEM,01)+PX2*HU1S3I+PY2*HV1S3I
00557         XM(IELEM,02)=XM(IELEM,02)+PX3*HU1S3I+PY3*HV1S3I
00558         XM(IELEM,03)=XM(IELEM,03)+PX1*HU1SI +PY1*HV1SI
00559         XM(IELEM,04)=XM(IELEM,04)+PX2*HU1SI +PY2*HV1SI
00560         XM(IELEM,05)=XM(IELEM,05)+PX3*HU1SI +PY3*HV1SI
00561         XM(IELEM,06)=XM(IELEM,06)+PX3*HU2S3I+PY3*HV2S3I
00562         XM(IELEM,07)=XM(IELEM,07)+PX1*HU2SI +PY1*HV2SI
00563         XM(IELEM,08)=XM(IELEM,08)+PX2*HU2SI +PY2*HV2SI
00564         XM(IELEM,09)=XM(IELEM,09)+PX3*HU2SI +PY3*HV2SI
00565         XM(IELEM,10)=XM(IELEM,10)+PX1*HU3SI +PY1*HV3SI
00566         XM(IELEM,11)=XM(IELEM,11)+PX2*HU3SI +PY2*HV3SI
00567         XM(IELEM,12)=XM(IELEM,12)+PX3*HU3SI +PY3*HV3SI
00568         XM(IELEM,13)=XM(IELEM,13)+PX2*HU1I3S+PY2*HV1I3S
00569         XM(IELEM,14)=XM(IELEM,14)+PX3*HU1I3S+PY3*HV1I3S
00570         XM(IELEM,15)=XM(IELEM,15)+PX3*HU2I3S+PY3*HV2I3S
00571         XM(IELEM,16)=XM(IELEM,16)+PX1*HU2S3I+PY1*HV2S3I
00572         XM(IELEM,17)=XM(IELEM,17)+PX1*HU3S3I+PY1*HV3S3I
00573         XM(IELEM,18)=XM(IELEM,18)+PX1*HU1SI +PY1*HV1SI
00574         XM(IELEM,19)=XM(IELEM,19)+PX1*HU2SI +PY1*HV2SI
00575         XM(IELEM,20)=XM(IELEM,20)+PX1*HU3SI +PY1*HV3SI
00576         XM(IELEM,21)=XM(IELEM,21)+PX2*HU3S3I+PY2*HV3S3I
00577         XM(IELEM,22)=XM(IELEM,22)+PX2*HU1SI +PY2*HV1SI
00578         XM(IELEM,23)=XM(IELEM,23)+PX2*HU2SI +PY2*HV2SI
00579         XM(IELEM,24)=XM(IELEM,24)+PX2*HU3SI +PY2*HV3SI
00580         XM(IELEM,25)=XM(IELEM,25)+PX3*HU1SI +PY3*HV1SI
00581         XM(IELEM,26)=XM(IELEM,26)+PX3*HU2SI +PY3*HV2SI
00582         XM(IELEM,27)=XM(IELEM,27)+PX3*HU3SI +PY3*HV3SI
00583         XM(IELEM,28)=XM(IELEM,28)+PX1*HU2I3S+PY1*HV2I3S
00584         XM(IELEM,29)=XM(IELEM,29)+PX1*HU3I3S+PY1*HV3I3S
00585         XM(IELEM,30)=XM(IELEM,30)+PX2*HU3I3S+PY2*HV3I3S
00586 !
00587       ENDDO
00588 !
00589       ENDIF
00590 !
00591 !-----------------------------------------------------------------------
00592 !
00593 !     WITH GENERALIZED SIGMA TRANSFORMATION
00594 !
00595       ELSE
00596 !
00597       IF(SPECAD) THEN
00598         IF (LNG.EQ.1) WRITE(LU,201)
00599         IF (LNG.EQ.2) WRITE(LU,202)
00600 201     FORMAT(1X,'MT05PP (BIEF) :',/,
00601      &         1X,'SIGMAG=T ET SPECAD=T : ',1I6,' CAS NON PREVU',/,
00602      &         1X,'NOM REEL DE U : ',A6)
00603 202     FORMAT(1X,'MT05PP (BIEF) :',/,
00604      &         1X,'SIGMAG=T AND SPECAD=T : ',1I6,' NOT IMPLEMENTED',/,
00605      &         1X,'REAL NAME OF U : ',A6)
00606         CALL PLANTE(1)
00607         STOP
00608       ENDIF
00609 !
00610       NELEM2 = NELEM/(NPLAN-1)
00611 !
00612       DO IELEM=1,NELEM
00613 !
00614         IELEM2 = MOD(IELEM-1,NELEM2) + 1
00615         IELEM3 = NELEM - NELEM2 + IELEM2
00616 !
00617         I1 = IKLE(IELEM,1)
00618         I2 = IKLE(IELEM,2)
00619         I3 = IKLE(IELEM,3)
00620         I4 = IKLE(IELEM,4)
00621         I5 = IKLE(IELEM,5)
00622         I6 = IKLE(IELEM,6)
00623 !
00624 !       COMPUTES DELTA Z
00625 !
00626         DZ1 = Z(I4) - Z(I1)
00627         DZ2 = Z(I5) - Z(I2)
00628         DZ3 = Z(I6) - Z(I3)
00629 !
00630 !       COMPUTES WATER DEPTH
00631 !
00632         H1 = Z(IKLE(IELEM3,4))-Z(IKLE(IELEM2,1))
00633         H2 = Z(IKLE(IELEM3,5))-Z(IKLE(IELEM2,2))
00634         H3 = Z(IKLE(IELEM3,6))-Z(IKLE(IELEM2,3))
00635 !
00636         HH1=MAX(H1,1.D-6)
00637         HH2=MAX(H2,1.D-6)
00638         HH3=MAX(H3,1.D-6)
00639 !
00640 !       X2  =  X(I2) - X(I1)
00641 !       X3  =  X(I3) - X(I1)
00642 !       Y2  =  Y(I2) - Y(I1)
00643 !       Y3  =  Y(I3) - Y(I1)
00644         X2  =  X(IELEM,2)
00645         X3  =  X(IELEM,3)
00646         Y2  =  Y(IELEM,2)
00647         Y3  =  Y(IELEM,3)
00648 !
00649         U1=DZ1*U(I1)/HH1
00650         U2=DZ2*U(I2)/HH2
00651         U3=DZ3*U(I3)/HH3
00652         U4=DZ1*U(I4)/HH1
00653         U5=DZ2*U(I5)/HH2
00654         U6=DZ3*U(I6)/HH3
00655         V1=DZ1*V(I1)/HH1
00656         V2=DZ2*V(I2)/HH2
00657         V3=DZ3*V(I3)/HH3
00658         V4=DZ1*V(I4)/HH1
00659         V5=DZ2*V(I5)/HH2
00660         V6=DZ3*V(I6)/HH3
00661 !
00662         Q1  =  W(I1)
00663         Q2  =  W(I2)
00664         Q3  =  W(I3)
00665 !
00666 !    INTERMEDIATE COMPUTATIONS
00667 !
00668         DEN=XMUL/1440.D0
00669         PX2=Y3*DEN
00670         PX3=-Y2*DEN
00671         PX1=-PX2-PX3
00672         PY2=-X3*DEN
00673         PY3=X2*DEN
00674         PY1=-PY3-PY2
00675 !
00676         HT = H1+H2+H3
00677         SUS = U6+U5+U4
00678         SUI = U3+U2+U1
00679         SHUI = H3*U3+H2*U2+H1*U1
00680         SHUS = H3*U6+H2*U5+H1*U4
00681         SHUS1 = (HT+H1)*(SUS+U4)+SHUS+U4*H1
00682         SHUI1 = (HT+H1)*(SUI+U1)+SHUI+U1*H1
00683         SHUS2 = (HT+H2)*(SUS+U5)+SHUS+U5*H2
00684         SHUI2 = (HT+H2)*(SUI+U2)+SHUI+U2*H2
00685         SHUS3 = (HT+H3)*(SUS+U6)+SHUS+U6*H3
00686         SHUI3 = (HT+H3)*(SUI+U3)+SHUI+U3*H3
00687         HU1S3I = SHUS1+3*SHUI1
00688         HU1I3S = SHUI1+3*SHUS1
00689         HU1SI = SHUS1+SHUI1
00690         HU2S3I = SHUS2+3*SHUI2
00691         HU2I3S = SHUI2+3*SHUS2
00692         HU2SI = SHUS2+SHUI2
00693         HU3S3I = SHUS3+3*SHUI3
00694         HU3I3S = SHUI3+3*SHUS3
00695         HU3SI = SHUS3+SHUI3
00696 !
00697         SVS = V6+V5+V4
00698         SVI = V3+V2+V1
00699         SHVI = H3*V3+H2*V2+H1*V1
00700         SHVS = H3*V6+H2*V5+H1*V4
00701         SHVS1  = (HT+H1)*(SVS+V4)+SHVS+V4*H1
00702         SHVI1  = (HT+H1)*(SVI+V1)+SHVI+V1*H1
00703         SHVS2  = (HT+H2)*(SVS+V5)+SHVS+V5*H2
00704         SHVI2  = (HT+H2)*(SVI+V2)+SHVI+V2*H2
00705         SHVS3  = (HT+H3)*(SVS+V6)+SHVS+V6*H3
00706         SHVI3  = (HT+H3)*(SVI+V3)+SHVI+V3*H3
00707         HV1S3I = SHVS1+3*SHVI1
00708         HV1I3S = SHVI1+3*SHVS1
00709         HV1SI  = SHVS1+  SHVI1
00710         HV2S3I = SHVS2+3*SHVI2
00711         HV2I3S = SHVI2+3*SHVS2
00712         HV2SI  = SHVS2+  SHVI2
00713         HV3S3I = SHVS3+3*SHVI3
00714         HV3I3S = SHVI3+3*SHVS3
00715         HV3SI  = SHVS3+  SHVI3
00716 !
00717 !       OPTION WITH MASS-LUMPING ON VERTICAL VELOCITIES
00718 !       COMPATIBILITY WITH TRIDW2 (LUMPING TO COMPUTE DELTA(Z)*WSTAR)
00719 !
00720         PZ1=-DEN*(X2*Y3-Y2*X3)*30
00721 !
00722         INT1=PZ1*2*Q1
00723         T(IELEM,1)=INT1
00724         XM(IELEM,3)=-INT1
00725         XM(IELEM,18)=INT1
00726         T(IELEM,4)=-INT1
00727 !
00728         INT1=PZ1*2*Q2
00729         T(IELEM,2)=INT1
00730         XM(IELEM,8)= -INT1
00731         XM(IELEM,23)=INT1
00732         T(IELEM,5)=-INT1
00733 !
00734         INT1=PZ1*2*Q3
00735         T(IELEM,3)=INT1
00736         XM(IELEM,12)=-INT1
00737         XM(IELEM,27)=INT1
00738         T(IELEM,6)=-INT1
00739 !
00740         INT1=PZ1*Q1
00741         XM(IELEM,16)=INT1
00742         XM(IELEM,7)=-INT1
00743         XM(IELEM,17)=INT1
00744         XM(IELEM,10)=-INT1
00745         XM(IELEM,19)=INT1
00746         XM(IELEM,28)=-INT1
00747         XM(IELEM,20)=INT1
00748         XM(IELEM,29)=-INT1
00749 !
00750         INT1=PZ1*Q2
00751         XM(IELEM,1)=INT1
00752         XM(IELEM,4)=-INT1
00753         XM(IELEM,22)=INT1
00754         XM(IELEM,13)=-INT1
00755         XM(IELEM,21)=INT1
00756         XM(IELEM,11)=-INT1
00757         XM(IELEM,24)=INT1
00758         XM(IELEM,30)=-INT1
00759 !
00760         INT1=PZ1*Q3
00761         XM(IELEM,2)=INT1
00762         XM(IELEM,5)=-INT1
00763         XM(IELEM,25)=INT1
00764         XM(IELEM,14)=-INT1
00765         XM(IELEM,26)=INT1
00766         XM(IELEM,15)=-INT1
00767         XM(IELEM,6)=INT1
00768         XM(IELEM,9)=-INT1
00769 !
00770         T(IELEM,1)=T(IELEM,1)+PX1*HU1S3I+PY1*HV1S3I
00771         T(IELEM,2)=T(IELEM,2)+PX2*HU2S3I+PY2*HV2S3I
00772         T(IELEM,3)=T(IELEM,3)+PX3*HU3S3I+PY3*HV3S3I
00773         T(IELEM,4)=T(IELEM,4)+PX1*HU1I3S+PY1*HV1I3S
00774         T(IELEM,5)=T(IELEM,5)+PX2*HU2I3S+PY2*HV2I3S
00775         T(IELEM,6)=T(IELEM,6)+PX3*HU3I3S+PY3*HV3I3S
00776 !
00777         XM(IELEM,01)=XM(IELEM,01)+PX2*HU1S3I+PY2*HV1S3I
00778         XM(IELEM,02)=XM(IELEM,02)+PX3*HU1S3I+PY3*HV1S3I
00779         XM(IELEM,03)=XM(IELEM,03)+PX1*HU1SI +PY1*HV1SI
00780         XM(IELEM,04)=XM(IELEM,04)+PX2*HU1SI +PY2*HV1SI
00781         XM(IELEM,05)=XM(IELEM,05)+PX3*HU1SI +PY3*HV1SI
00782         XM(IELEM,06)=XM(IELEM,06)+PX3*HU2S3I+PY3*HV2S3I
00783         XM(IELEM,07)=XM(IELEM,07)+PX1*HU2SI +PY1*HV2SI
00784         XM(IELEM,08)=XM(IELEM,08)+PX2*HU2SI +PY2*HV2SI
00785         XM(IELEM,09)=XM(IELEM,09)+PX3*HU2SI +PY3*HV2SI
00786         XM(IELEM,10)=XM(IELEM,10)+PX1*HU3SI +PY1*HV3SI
00787         XM(IELEM,11)=XM(IELEM,11)+PX2*HU3SI +PY2*HV3SI
00788         XM(IELEM,12)=XM(IELEM,12)+PX3*HU3SI +PY3*HV3SI
00789         XM(IELEM,13)=XM(IELEM,13)+PX2*HU1I3S+PY2*HV1I3S
00790         XM(IELEM,14)=XM(IELEM,14)+PX3*HU1I3S+PY3*HV1I3S
00791         XM(IELEM,15)=XM(IELEM,15)+PX3*HU2I3S+PY3*HV2I3S
00792         XM(IELEM,16)=XM(IELEM,16)+PX1*HU2S3I+PY1*HV2S3I
00793         XM(IELEM,17)=XM(IELEM,17)+PX1*HU3S3I+PY1*HV3S3I
00794         XM(IELEM,18)=XM(IELEM,18)+PX1*HU1SI +PY1*HV1SI
00795         XM(IELEM,19)=XM(IELEM,19)+PX1*HU2SI +PY1*HV2SI
00796         XM(IELEM,20)=XM(IELEM,20)+PX1*HU3SI +PY1*HV3SI
00797         XM(IELEM,21)=XM(IELEM,21)+PX2*HU3S3I+PY2*HV3S3I
00798         XM(IELEM,22)=XM(IELEM,22)+PX2*HU1SI +PY2*HV1SI
00799         XM(IELEM,23)=XM(IELEM,23)+PX2*HU2SI +PY2*HV2SI
00800         XM(IELEM,24)=XM(IELEM,24)+PX2*HU3SI +PY2*HV3SI
00801         XM(IELEM,25)=XM(IELEM,25)+PX3*HU1SI +PY3*HV1SI
00802         XM(IELEM,26)=XM(IELEM,26)+PX3*HU2SI +PY3*HV2SI
00803         XM(IELEM,27)=XM(IELEM,27)+PX3*HU3SI +PY3*HV3SI
00804         XM(IELEM,28)=XM(IELEM,28)+PX1*HU2I3S+PY1*HV2I3S
00805         XM(IELEM,29)=XM(IELEM,29)+PX1*HU3I3S+PY1*HV3I3S
00806         XM(IELEM,30)=XM(IELEM,30)+PX2*HU3I3S+PY2*HV3I3S
00807 !
00808       ENDDO
00809 !
00810       ENDIF
00811 !
00812 !-----------------------------------------------------------------------
00813 !
00814       RETURN
00815       END

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