mt14pp.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\mt14pp.f
00002 !
00086                      SUBROUTINE MT14PP
00087 !                    *****************
00088 !
00089      &( T,XM,PPQ,LEGO,XMUL,SU,SV,SW,U,V,W,SF,SG,SH,F,G,H,
00090      &  SURFAC,IKLE,NELEM,NELMAX,SIGMAG,SPECAD)
00091 !
00092 !***********************************************************************
00093 ! BIEF   V6P3                                  21/08/2010
00094 !***********************************************************************
00095 !
00096 !
00097 !
00098 !reference  J-M HERVOUET THESIS OR BOOK: MURD SCHEME IN 3 DIMENSIONS
00099 !+             CHAPTER 6.
00100 !
00101 !
00102 !
00103 !
00104 !
00105 !
00106 !
00107 !
00108 !
00109 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00110 !| F              |-->| FUNCTION USED IN THE FORMULA
00111 !| G              |-->| FUNCTION USED IN THE FORMULA
00112 !| H              |-->| FUNCTION USED IN THE FORMULA
00113 !| IKLE           |-->| CONNECTIVITY TABLE
00114 !| LEGO           |-->| LOGICAL. IF YES RESULT ASSEMBLED.
00115 !| NELEM          |-->| NUMBER OF ELEMENTS
00116 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00117 !| PPQ            |-->| STORAGE IN XM OF OFF-DIAGONAL ELEMENTS
00118 !| SF             |-->| BIEF_OBJ STRUCTURE OF F
00119 !| SG             |-->| BIEF_OBJ STRUCTURE OF G
00120 !| SH             |-->| BIEF_OBJ STRUCTURE OF H
00121 !| SIGMAG         |-->| IF YES, GENERALISED SIGMA TRANSFORMATION
00122 !| SPECAD         |-->| IF YES, SPECIAL FORM OF ADVECTION FIELD
00123 !| SU             |-->| BIEF_OBJ STRUCTURE OF U
00124 !| SV             |-->| BIEF_OBJ STRUCTURE OF V
00125 !| SW             |-->| BIEF_OBJ STRUCTURE OF W
00126 !| SURFAC         |-->| AREA OF 2D ELEMENTS
00127 !| T              |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
00128 !| U              |-->| FUNCTION USED IN THE FORMULA
00129 !| V              |-->| FUNCTION USED IN THE FORMULA
00130 !| W              |-->| FUNCTION USED IN THE FORMULA
00131 !| XM             |<->| OFF-DIAGONAL TERMS
00132 !| XMUL           |-->| COEFFICIENT FOR MULTIPLICATION
00133 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00134 !
00135       USE BIEF, EX_MT14PP => MT14PP
00136 !
00137       IMPLICIT NONE
00138       INTEGER LNG,LU
00139       COMMON/INFO/LNG,LU
00140 !
00141 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00142 !
00143       INTEGER, INTENT(IN) :: NELEM,NELMAX
00144       INTEGER, INTENT(IN) :: IKLE(NELMAX,6),PPQ(6,6)
00145 !
00146       DOUBLE PRECISION, INTENT(INOUT) :: T(NELMAX,6),XM(30,NELMAX)
00147 !
00148       DOUBLE PRECISION, INTENT(IN) :: XMUL
00149       DOUBLE PRECISION, INTENT(IN) :: U(*),V(*),W(*),F(*),G(*)
00150       DOUBLE PRECISION, INTENT(IN) :: H(NELMAX,6)
00151 !
00152       LOGICAL, INTENT(IN) :: LEGO,SIGMAG,SPECAD
00153 !
00154 !     STRUCTURES DE U,V,W
00155 !
00156       TYPE(BIEF_OBJ), INTENT(IN) :: SU,SV,SW,SF,SG,SH
00157 !
00158       DOUBLE PRECISION, INTENT(IN) :: SURFAC(NELMAX)
00159 !
00160 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00161 !
00162 !     DECLARATIONS SPECIFIC TO THIS SUBROUTINE
00163 !
00164       INTEGER IELEM,IELMW
00165       DOUBLE PRECISION XSUR3
00166       DOUBLE PRECISION LI0J0,LI0K0,LJ0I0,LJ0K0,LK0I0,LK0J0
00167       DOUBLE PRECISION LI3J3,LI3K3,LJ3I3,LJ3K3,LK3I3,LK3J3
00168       DOUBLE PRECISION LI3J0,LI3K0,LJ3I0,LJ3K0,LK3I0,LK3J0
00169       DOUBLE PRECISION LI3I0,LJ3J0,LK3K0,LI0I3
00170       DOUBLE PRECISION ALFA,ALFAJ,ALFAK,ALFA1,ALFA2,BETA,SOM0,SOM3
00171       DOUBLE PRECISION ALFAI0,ALFAJ0,ALFAK0,ALFAI3,ALFAJ3,ALFAK3
00172       DOUBLE PRECISION ALFAII,ALFAIJ,ALFAIK,ALFAJI,ALFAJJ,ALFAJK
00173       DOUBLE PRECISION ALFAKI,ALFAKJ,ALFAKK,EPS
00174       DATA EPS /1.D-10/
00175 !
00176       INTEGER IPLUS1(6),IPLUS3(6),INDIC(0:7)
00177       INTEGER IXM,I0,J0,K0,I3,J3,K3
00178 !
00179 !-----------------------------------------------------------------------
00180 !
00181 !     DONNEES
00182 !
00183       DATA IPLUS1 / 2 , 3 , 1 , 5 , 6 , 4 /
00184       DATA IPLUS3 / 4 , 5 , 6 , 1 , 2 , 3 /
00185       DATA INDIC  / 4 , 4 , 5 , 3 , 6 , 2 , 1 , 1 /
00186 !
00187 !=======================================================================
00188 !
00189       IELMW = SW%ELM
00190 !                                                         *   *
00191 !     COMPUTES THE COEFFICIENTS A(I) (INTEGRAL OF H U.GRAD(PSI), ONLY
00192 !     CONSIDERING THE HORIZONTAL PART OF U). SEE JMH THESIS
00193 !
00194 !     NOTE: THIS VC04PP CALL IS ALREADY DONE BY FLUX3D IN TELEMAC-3D
00195 !           THROUGH A CALL TO VECTOR. THE EQUIVALENT OF T HAS BEEN
00196 !           SAVED IN WHAT IS HERE H
00197 !
00198 !     FORMUL='             HOR'
00199 !     CALL VC04PP(XMUL,SU,SV,SW,U,V,W,SF,SG,SH,F,G,H,X,Y,Z,
00200 !    * IKLE(1,1),IKLE(1,2),IKLE(1,3),IKLE(1,4),IKLE(1,5),IKLE(1,6),
00201 !    * NELEM,NELMAX,T(1,1),T(1,2),T(1,3),T(1,4),T(1,5),T(1,6),
00202 !    * SPECAD,FORMUL,NPLAN)
00203 !
00204 !     NOW T IS RETRIEVED BY STORAGE DONE INTO FLUX3D
00205 !     FUNCTION H MUST CONTAIN THE EQUIVALENT OF T ABOVE, WHICH IS
00206 !     THE NON ASSEMBLED VALUE OF FLUINT
00207 !
00208 !     DO IELEM=1,NELEM
00209 !       T(IELEM,1)=H(         IELEM)
00210 !       T(IELEM,2)=H(  NELMAX+IELEM)
00211 !       T(IELEM,3)=H(2*NELMAX+IELEM)
00212 !       T(IELEM,4)=H(3*NELMAX+IELEM)
00213 !       T(IELEM,5)=H(4*NELMAX+IELEM)
00214 !       T(IELEM,6)=H(5*NELMAX+IELEM)
00215 !     ENDDO
00216 !
00217 !     FORMUL NORMALLY STARTS WITH VGRADP OR VGRADP2, OR VGRADP22 TO KNOW
00218 !     SIGMAG AND SPECAD, USELESS HERE.
00219       XSUR3 = XMUL/3.D0
00220 !
00221 !-----------------------------------------------------------------------
00222 !
00223       IF(IELMW.NE.41) THEN
00224         IF(LNG.EQ.1) WRITE(LU,*) 'MT14PP DISCRETISATION NON PREVUE'
00225         IF(LNG.EQ.2) WRITE(LU,*) 'MT14PP DISCRETISATION NOT HANDLED'
00226         CALL PLANTE(1)
00227         STOP
00228       ENDIF
00229 !
00230 !     COMPUTES THE COEFFICIENTS B(I) BY ELEMENTS
00231 !     AND COMPUTES LAMBDA(I,J)
00232 !
00233       DO IELEM = 1 , NELEM
00234 !
00235         IXM = 0
00236         IF (W(IKLE(IELEM,1)).GT.0.D0) IXM = 1
00237         IF (W(IKLE(IELEM,2)).GT.0.D0) IXM = IXM + 2
00238         IF (W(IKLE(IELEM,3)).GT.0.D0) IXM = IXM + 4
00239 !
00240 !       NOTE JMH:
00241 !       IT SEEMS THAT INDIC IS THE POINT THAT RECEIVES THE VERTICAL
00242 !       VELOCITY WHICH IS THE ONLY ONE OF ITS SIGN
00243 !
00244 !       ALL W < 0         INDIC = 4 (RANDOM, WHY NOT 1 ?)
00245 !       ONLY W1 > 0       INDIC = 4
00246 !       ONLY W2 > 0       INCIC = 5
00247 !       ONLY W3 > 0       INDIC = 6
00248 !       ALL W > 0         INDIC = 1 (RANDOM, WHY NOT 4 ?)
00249 !       ONLY W1 < 0       INDIC = 1
00250 !       ONLY W2 < 0       INDIC = 2
00251 !       ONLY W3 < 0       INDIC = 3
00252 !
00253         I0 = INDIC(IXM)
00254         J0 = IPLUS1(I0)
00255         K0 = IPLUS1(J0)
00256         I3 = IPLUS3(I0)
00257         J3 = IPLUS1(I3)
00258         K3 = IPLUS1(J3)
00259 !
00260         ALFA = XSUR3 * SURFAC(IELEM)
00261 !
00262         LI3I0 = ALFA * ABS(W(IKLE(IELEM,MIN(I0,I3))))
00263         LI0I3 = 0.D0
00264         IF (IXM.GE.1.AND.IXM.LE.6) THEN
00265            LI0I3 = LI3I0
00266            LI3I0 = 0.D0
00267         ENDIF
00268         LJ3J0 = ALFA * ABS(W(IKLE(IELEM,MIN(J0,J3))))
00269         XM(PPQ(J0,J3),IELEM) = 0.D0
00270         LK3K0 = ALFA * ABS(W(IKLE(IELEM,MIN(K0,K3))))
00271         XM(PPQ(K0,K3),IELEM) = 0.D0
00272 !
00273         ALFAI0 = -H(IELEM,I0)
00274         ALFAJ0 = -H(IELEM,J0)
00275         ALFAK0 = -H(IELEM,K0)
00276         ALFAI3 = -H(IELEM,I3)
00277         ALFAJ3 = -H(IELEM,J3)
00278         ALFAK3 = -H(IELEM,K3)
00279 !
00280         LJ0I0 = MAX(0.D0,MIN(ALFAI0,-ALFAJ0))
00281         LK0I0 = MAX(0.D0,MIN(ALFAI0,-ALFAK0))
00282         LI3J3 = MAX(0.D0,MIN(ALFAJ3,-ALFAI3))
00283         LI3K3 = MAX(0.D0,MIN(ALFAK3,-ALFAI3))
00284 !
00285         ALFAJ = MIN(LJ3J0,MAX(LI3J3,LJ0I0))
00286         ALFAK = MIN(LK3K0,MAX(LI3K3,LK0I0))
00287         ALFA  = MIN(LI0I3,ALFAJ+ALFAK)
00288         IF (ALFA.GT.EPS) THEN
00289 !
00290 !          JMH CODE
00291 !
00292 !          IF(ALFA.LT.ALFAJ+ALFAK) THEN
00293 !            LIMITATION OF ALFAJ AND ALFAK WITH THE IDEA THAT THE
00294 !            TWO CORRECTED FLUXES WILL BE AS EQUAL AS POSSIBLE
00295 !            MOREOVER WE MUST HAVE ALFAJ+ALFAK=ALFA
00296 !            IF(ALFAK.GT.ALFAJ) THEN
00297 !              DELTA=MIN(ALFA,ALFAK-ALFAJ)
00298 !              ALFAK=0.5D0*(ALFA+DELTA)
00299 !              ALFAJ=0.5D0*(ALFA-DELTA)
00300 !            ELSE
00301 !              DELTA=MIN(ALFA,ALFAJ-ALFAK)
00302 !              ALFAK=0.5D0*(ALFA-DELTA)
00303 !              ALFAJ=0.5D0*(ALFA+DELTA)
00304 !            ENDIF
00305 !          ENDIF
00306 !
00307 !          JMJ CODE
00308            BETA = ALFA / (ALFAJ+ALFAK)
00309            ALFAJ = ALFAJ * BETA
00310            ALFAK = ALFAK * BETA
00311 !          END OF JMJ CODE
00312 !
00313            LI0I3 = LI0I3-ALFA
00314            LJ3J0 = LJ3J0-ALFAJ
00315            LK3K0 = LK3K0-ALFAK
00316 !
00317            ALFAI3 = ALFAI3+ALFA
00318            ALFAI0 = ALFAI0-ALFA
00319            ALFAJ0 = ALFAJ0+ALFAJ
00320            ALFAJ3 = ALFAJ3-ALFAJ
00321            ALFAK0 = ALFAK0+ALFAK
00322            ALFAK3 = ALFAK3-ALFAK
00323 !
00324            LJ0I0 = MAX(0.D0,MIN(ALFAI0,-ALFAJ0))
00325            LK0I0 = MAX(0.D0,MIN(ALFAI0,-ALFAK0))
00326            LI3J3 = MAX(0.D0,MIN(ALFAJ3,-ALFAI3))
00327            LI3K3 = MAX(0.D0,MIN(ALFAK3,-ALFAI3))
00328 !
00329         ENDIF
00330 !
00331         LI0J0 = MAX(0.D0,MIN(ALFAJ0,-ALFAI0))
00332         LK0J0 = MAX(0.D0,MIN(ALFAJ0,-ALFAK0))
00333         LI0K0 = MAX(0.D0,MIN(ALFAK0,-ALFAI0))
00334         LJ0K0 = MAX(0.D0,MIN(ALFAK0,-ALFAJ0))
00335         LJ3I3 = MAX(0.D0,MIN(ALFAI3,-ALFAJ3))
00336         LJ3K3 = MAX(0.D0,MIN(ALFAK3,-ALFAJ3))
00337         LK3I3 = MAX(0.D0,MIN(ALFAI3,-ALFAK3))
00338         LK3J3 = MAX(0.D0,MIN(ALFAJ3,-ALFAK3))
00339 !
00340         XM(PPQ(J0,I3),IELEM) = 0.D0
00341         XM(PPQ(K0,I3),IELEM) = 0.D0
00342         XM(PPQ(I0,J3),IELEM) = 0.D0
00343         XM(PPQ(K0,J3),IELEM) = 0.D0
00344         XM(PPQ(I0,K3),IELEM) = 0.D0
00345         XM(PPQ(J0,K3),IELEM) = 0.D0
00346 !
00347         ALFAI0 = 0.D0
00348         ALFAJ0 = 0.D0
00349         ALFAK0 = 0.D0
00350         ALFAI3 = 0.D0
00351         ALFAJ3 = 0.D0
00352         ALFAK3 = 0.D0
00353 !
00354         SOM0 = LJ0I0+LK0I0
00355         SOM3 = LI3J3+LI3K3
00356 !
00357         ALFA1 = MIN(LI0I3,SOM0)
00358         IF (ALFA1.GT.EPS) THEN
00359 !
00360            BETA = 1.D0 / SOM0
00361            ALFAJ0 = LJ0I0 * BETA
00362            ALFAK0 = LK0I0 * BETA
00363            ALFA = MAX(0.D0,ALFA1-SOM3)
00364 !
00365            LJ0I0 = LJ0I0 - ALFAJ0*ALFA1
00366            LK0I0 = LK0I0 - ALFAK0*ALFA1
00367            XM(PPQ(J0,I3),IELEM) = ALFAJ0*ALFA
00368            XM(PPQ(K0,I3),IELEM) = ALFAK0*ALFA
00369 !
00370         ENDIF
00371 !
00372         ALFA2 = MIN(LI0I3,SOM3)
00373         IF (ALFA2.GT.EPS) THEN
00374 !
00375            BETA = 1.D0 / SOM3
00376            ALFAJ3 = LI3J3 * BETA
00377            ALFAK3 = LI3K3 * BETA
00378            ALFA = MAX(0.D0,ALFA2-SOM0)
00379 !
00380            LI3J3 = LI3J3 - ALFAJ3*ALFA2
00381            LI3K3 = LI3K3 - ALFAK3*ALFA2
00382            XM(PPQ(I0,J3),IELEM) = ALFAJ3*ALFA
00383            XM(PPQ(I0,K3),IELEM) = ALFAK3*ALFA
00384 !
00385         ENDIF
00386 !
00387         ALFA = MIN(ALFA1,ALFA2)
00388         IF (ALFA.GT.0.D0) THEN
00389 !
00390            ALFAJ3 = ALFAJ3 * ALFA
00391            ALFAK3 = ALFAK3 * ALFA
00392            ALFAJJ = ALFAJ3 * ALFAJ0
00393            ALFAKK = ALFAK3 * ALFAK0
00394            ALFAJK = ALFAJ3 * ALFAK0
00395            ALFAKJ = ALFAK3 * ALFAJ0
00396            ALFA = MIN(ALFAJK,ALFAKJ)
00397 !
00398            XM(PPQ(J0,J3),IELEM) = ALFAJJ + ALFA
00399            XM(PPQ(K0,K3),IELEM) = ALFAKK + ALFA
00400            XM(PPQ(K0,J3),IELEM) = ALFAJK - ALFA
00401            XM(PPQ(J0,K3),IELEM) = ALFAKJ - ALFA
00402 !
00403         ENDIF
00404 !
00405         XM(PPQ(I0,I3),IELEM) = MAX(0.D0,LI0I3-MAX(SOM3,SOM0))
00406 !
00407         LJ3I0 = 0.D0
00408         LK3I0 = 0.D0
00409         LI3J0 = 0.D0
00410         LK3J0 = 0.D0
00411         LI3K0 = 0.D0
00412         LJ3K0 = 0.D0
00413 !
00414         SOM0 = LI0J0+LI0K0
00415         SOM3 = LJ3I3+LK3I3
00416 !
00417         ALFA1 = MIN(LI3I0,SOM0)
00418         IF (ALFA1.GT.EPS) THEN
00419 !
00420            BETA = 1.D0 / SOM0
00421            ALFAJ0 = LI0J0 * BETA
00422            ALFAK0 = LI0K0 * BETA
00423            ALFA = MAX(0.D0,ALFA1-SOM3)
00424 !
00425            LI0J0 = LI0J0 - ALFAJ0*ALFA1
00426            LI0K0 = LI0K0 - ALFAK0*ALFA1
00427            LI3J0 = LI3J0 + ALFAJ0*ALFA
00428            LI3K0 = LI3K0 + ALFAK0*ALFA
00429 !
00430         ENDIF
00431 !
00432         XM(PPQ(I0,J0),IELEM) = LI0J0
00433         XM(PPQ(I0,K0),IELEM) = LI0K0
00434 !
00435         ALFA2 = MIN(LI3I0,SOM3)
00436         IF (ALFA2.GT.EPS) THEN
00437 !
00438            BETA = 1.D0 / SOM3
00439            ALFAJ3 = LJ3I3 * BETA
00440            ALFAK3 = LK3I3 * BETA
00441            ALFA = MAX(0.D0,ALFA2-SOM0)
00442 !
00443            LJ3I3 = LJ3I3 - ALFAJ3*ALFA2
00444            LK3I3 = LK3I3 - ALFAK3*ALFA2
00445            LJ3I0 = LJ3I0 + ALFAJ3*ALFA
00446            LK3I0 = LK3I0 + ALFAK3*ALFA
00447 !
00448         ENDIF
00449 !
00450         XM(PPQ(J3,I3),IELEM) = LJ3I3
00451         XM(PPQ(K3,I3),IELEM) = LK3I3
00452 !
00453         ALFA = MIN(ALFA1,ALFA2)
00454         IF (ALFA.GT.0.D0) THEN
00455 !
00456            ALFAJ0 = ALFAJ0 * ALFA
00457            ALFAK0 = ALFAK0 * ALFA
00458            ALFAJJ = ALFAJ0 * ALFAJ3
00459            ALFAKK = ALFAK0 * ALFAK3
00460            ALFAJK = ALFAJ0 * ALFAK3
00461            ALFAKJ = ALFAK0 * ALFAJ3
00462            ALFA = MIN(ALFAJK,ALFAKJ)
00463 !
00464            LJ3J0 = LJ3J0 + ALFAJJ + ALFA
00465            LK3K0 = LK3K0 + ALFAKK + ALFA
00466            LK3J0 = LK3J0 + ALFAJK - ALFA
00467            LJ3K0 = LJ3K0 + ALFAKJ - ALFA
00468 !
00469         ENDIF
00470 !
00471         LI3I0 = MAX(0.D0,LI3I0-MAX(SOM0,SOM3))
00472 !
00473         SOM0 = LJ0I0+LJ0K0
00474         SOM3 = LI3J3+LK3J3
00475 !
00476         ALFA1 = MIN(LJ3J0,SOM0)
00477         IF (ALFA1.GT.EPS) THEN
00478 !
00479            BETA = 1.D0 / SOM0
00480            ALFAI0 = LJ0I0 * BETA
00481            ALFAK0 = LJ0K0 * BETA
00482            ALFA = MAX(0.D0,ALFA1-SOM3)
00483 !
00484            LJ0I0 = LJ0I0 - ALFAI0*ALFA1
00485            LJ0K0 = LJ0K0 - ALFAK0*ALFA1
00486            LJ3I0 = LJ3I0 + ALFAI0*ALFA
00487            LJ3K0 = LJ3K0 + ALFAK0*ALFA
00488 !
00489         ENDIF
00490 !
00491         XM(PPQ(J0,I0),IELEM) = LJ0I0
00492         XM(PPQ(J0,K0),IELEM) = LJ0K0
00493 !
00494         ALFA2 = MIN(LJ3J0,SOM3)
00495         IF (ALFA2.GT.EPS) THEN
00496 !
00497            BETA = 1.D0 / SOM3
00498            ALFAI3 = LI3J3 * BETA
00499            ALFAK3 = LK3J3 * BETA
00500            ALFA = MAX(0.D0,ALFA2-SOM0)
00501 !
00502            LI3J3 = LI3J3 - ALFAI3*ALFA2
00503            LK3J3 = LK3J3 - ALFAK3*ALFA2
00504            LI3J0 = LI3J0 + ALFAI3*ALFA
00505            LK3J0 = LK3J0 + ALFAK3*ALFA
00506 !
00507         ENDIF
00508 !
00509         XM(PPQ(I3,J3),IELEM) = LI3J3
00510         XM(PPQ(K3,J3),IELEM) = LK3J3
00511 !
00512         ALFA = MIN(ALFA1,ALFA2)
00513         IF (ALFA.GT.0.D0) THEN
00514 !
00515            ALFAI0 = ALFAI0 * ALFA
00516            ALFAK0 = ALFAK0 * ALFA
00517            ALFAII = ALFAI0 * ALFAI3
00518            ALFAKK = ALFAK0 * ALFAK3
00519            ALFAIK = ALFAI0 * ALFAK3
00520            ALFAKI = ALFAK0 * ALFAI3
00521            ALFA = MIN(ALFAIK,ALFAKI)
00522 !
00523            LI3I0 = LI3I0 + ALFAII + ALFA
00524            LK3K0 = LK3K0 + ALFAKK + ALFA
00525            LK3I0 = LK3I0 + ALFAIK - ALFA
00526            LI3K0 = LI3K0 + ALFAKI - ALFA
00527 !
00528         ENDIF
00529 !
00530         LJ3J0 = MAX(0.D0,LJ3J0-MAX(SOM0,SOM3))
00531 !
00532         SOM0 = LK0I0+LK0J0
00533         SOM3 = LI3K3+LJ3K3
00534 !
00535         ALFA1 = MIN(LK3K0,SOM0)
00536         IF (ALFA1.GT.EPS) THEN
00537 !
00538            BETA = 1.D0 / SOM0
00539            ALFAI0 = LK0I0 * BETA
00540            ALFAJ0 = LK0J0 * BETA
00541            ALFA = MAX(0.D0,ALFA1-SOM3)
00542 !
00543            LK0I0 = LK0I0 - ALFAI0*ALFA1
00544            LK0J0 = LK0J0 - ALFAJ0*ALFA1
00545            LK3I0 = LK3I0 + ALFAI0*ALFA
00546            LK3J0 = LK3J0 + ALFAJ0*ALFA
00547 !
00548         ENDIF
00549 !
00550         XM(PPQ(K0,I0),IELEM) = LK0I0
00551         XM(PPQ(K0,J0),IELEM) = LK0J0
00552 !
00553         ALFA2 = MIN(LK3K0,SOM3)
00554         IF (ALFA2.GT.EPS) THEN
00555 !
00556            BETA = 1.D0 / SOM3
00557            ALFAI3 = LI3K3 * BETA
00558            ALFAJ3 = LJ3K3 * BETA
00559            ALFA = MAX(0.D0,ALFA2-SOM0)
00560 !
00561            LI3K3 = LI3K3 - ALFAI3*ALFA2
00562            LJ3K3 = LJ3K3 - ALFAJ3*ALFA2
00563            LI3K0 = LI3K0 + ALFAI3*ALFA
00564            LJ3K0 = LJ3K0 + ALFAJ3*ALFA
00565 !
00566         ENDIF
00567 !
00568         XM(PPQ(I3,K3),IELEM) = LI3K3
00569         XM(PPQ(J3,K3),IELEM) = LJ3K3
00570 !
00571         ALFA = MIN(ALFA1,ALFA2)
00572         IF (ALFA.GT.0.D0) THEN
00573 !
00574            ALFAI0 = ALFAI0 * ALFA
00575            ALFAJ0 = ALFAJ0 * ALFA
00576            ALFAII = ALFAI0 * ALFAI3
00577            ALFAJJ = ALFAJ0 * ALFAJ3
00578            ALFAIJ = ALFAI0 * ALFAJ3
00579            ALFAJI = ALFAJ0 * ALFAI3
00580            ALFA = MIN(ALFAIJ,ALFAJI)
00581 !
00582            LI3I0 = LI3I0 + ALFAII + ALFA
00583            LJ3J0 = LJ3J0 + ALFAJJ + ALFA
00584            LJ3I0 = LJ3I0 + ALFAIJ - ALFA
00585            LI3J0 = LI3J0 + ALFAJI - ALFA
00586 !
00587         ENDIF
00588 !
00589         LK3K0 = MAX(0.D0,LK3K0-MAX(SOM0,SOM3))
00590 !
00591         XM(PPQ(I3,I0),IELEM) = LI3I0
00592         XM(PPQ(J3,I0),IELEM) = LJ3I0
00593         XM(PPQ(K3,I0),IELEM) = LK3I0
00594         XM(PPQ(I3,J0),IELEM) = LI3J0
00595         XM(PPQ(J3,J0),IELEM) = LJ3J0
00596         XM(PPQ(K3,J0),IELEM) = LK3J0
00597         XM(PPQ(I3,K0),IELEM) = LI3K0
00598         XM(PPQ(J3,K0),IELEM) = LJ3K0
00599         XM(PPQ(K3,K0),IELEM) = LK3K0
00600 !
00601       ENDDO ! IELEM
00602 !
00603 !-----------------------------------------------------------------------
00604 !
00605 !     COMPUTES THE SUM OF EACH ROW (WITH A - SIGN)
00606 !     THE DIAGONAL TERMS ARE 0
00607 !
00608       IF(LEGO) THEN
00609 !
00610         DO IELEM = 1,NELEM
00611 !
00612           T(IELEM,1) = -XM(01,IELEM)-XM(02,IELEM)
00613      &                 -XM(03,IELEM)-XM(04,IELEM)-XM(05,IELEM)
00614           T(IELEM,2) = -XM(16,IELEM)-XM(06,IELEM)
00615      &                 -XM(07,IELEM)-XM(08,IELEM)-XM(09,IELEM)
00616           T(IELEM,3) = -XM(17,IELEM)-XM(21,IELEM)
00617      &                 -XM(10,IELEM)-XM(11,IELEM)-XM(12,IELEM)
00618           T(IELEM,4) = -XM(18,IELEM)-XM(22,IELEM)
00619      &                 -XM(25,IELEM)-XM(13,IELEM)-XM(14,IELEM)
00620           T(IELEM,5) = -XM(19,IELEM)-XM(23,IELEM)
00621      &                 -XM(26,IELEM)-XM(28,IELEM)-XM(15,IELEM)
00622           T(IELEM,6) = -XM(20,IELEM)-XM(24,IELEM)
00623      &                 -XM(27,IELEM)-XM(29,IELEM)-XM(30,IELEM)
00624 !
00625         ENDDO
00626 !
00627       ENDIF
00628 !
00629 !-----------------------------------------------------------------------
00630 !
00631       RETURN
00632       END

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