mt02pp.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\mt02pp.f
00002 !
00110                      SUBROUTINE MT02PP
00111 !                    *****************
00112 !
00113      &(T,XM,XMUL,SF,SG,SH,F,G,H,X,Y,Z,SURFAC,IKLE,NELEM,NELMAX,INCHYD,
00114      & FORMUL,NPLAN)
00115 !
00116 !***********************************************************************
00117 ! BIEF   V6P2                                   21/08/2010
00118 !***********************************************************************
00119 !
00120 !
00121 !
00122 !
00123 !
00124 !
00125 !
00126 !
00127 !
00128 !
00129 !
00130 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00131 !| F              |-->| FUNCTION USED IN THE FORMULA
00132 !| FORMUL         |-->| FORMULA DESCRIBING THE RESULTING MATRIX
00133 !| G              |-->| FUNCTION USED IN THE FORMULA
00134 !| H              |-->| FUNCTION USED IN THE FORMULA
00135 !| IKLE           |-->| CONNECTIVITY TABLE.
00136 !| INCHYD         |---| IF YES, TREATS HYDROSTATIC INCONSISTENCIES
00137 !| NELEM          |-->| NUMBER OF ELEMENTS
00138 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00139 !| NPLAN          |-->| NUMBER OF PLANES IN THE MESH OF PRISMS
00140 !| SF             |-->| STRUCTURE OF FUNCTIONS F
00141 !| SG             |-->| STRUCTURE OF FUNCTIONS G
00142 !| SH             |-->| STRUCTURE OF FUNCTIONS H
00143 !| SURFAC         |-->| AREA OF 2D ELEMENTS
00144 !| T              |<->| WORK ARRAY FOR ELEMENT BY ELEMENT DIAGONAL
00145 !| X              |-->| ABSCISSAE OF POINTS, PER ELEMENT
00146 !| Y              |-->| ORDINATES OF POINTS, PER ELEMENT
00147 !| Z              |-->| ELEVATIONS OF POINTS, PER POINT
00148 !| XM             |<->| OFF-DIAGONAL TERMS
00149 !| XMUL           |-->| COEFFICIENT FOR MULTIPLICATION
00150 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00151 !
00152       USE BIEF, EX_MT02PP => MT02PP
00153 !
00154       IMPLICIT NONE
00155       INTEGER LNG,LU
00156       COMMON/INFO/LNG,LU
00157 !
00158 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00159 !
00160       INTEGER, INTENT(IN) :: NELEM,NELMAX,NPLAN
00161       INTEGER, INTENT(IN) :: IKLE(NELMAX,6)
00162 !
00163       DOUBLE PRECISION, INTENT(INOUT) :: T(NELMAX,6),XM(NELMAX,15)
00164       DOUBLE PRECISION, INTENT(IN)    :: SURFAC(NELMAX)
00165 !
00166       DOUBLE PRECISION, INTENT(IN)    :: XMUL
00167       DOUBLE PRECISION, INTENT(IN)    :: F(*),G(*),H(*)
00168 !
00169 !     STRUCTURES OF F, G, H
00170 !
00171       TYPE(BIEF_OBJ), INTENT(IN)      :: SF,SG,SH
00172 !
00173       DOUBLE PRECISION, INTENT(IN)    :: X(NELMAX,6),Y(NELMAX,6),Z(*)
00174 !
00175       LOGICAL, INTENT(IN)             :: INCHYD
00176       CHARACTER(LEN=16), INTENT(IN)   :: FORMUL
00177 !
00178 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00179 !
00180 !     DECLARATIONS SPECIFIC TO THIS SUBROUTINE
00181 !
00182       DOUBLE PRECISION H1,H2,H3,RI,RS,R,RRI,RRS,RRRI,RRRS
00183       DOUBLE PRECISION D1,D2,D3,D12,D13,D23,D
00184       DOUBLE PRECISION SH1,SHH,SNI,SNS,SNHI,SNHS,SNHHI,SNHHS,SNHH
00185       DOUBLE PRECISION SNHI1,SNHS1,SNHI2,SNHS2,SNHI3,SNHS3
00186       DOUBLE PRECISION HHI12,HHS12,HH12,HHI13,HHS13,HH13
00187       DOUBLE PRECISION HHI23,HHS23,HH23
00188       DOUBLE PRECISION HRI1,HRS1,HRI2,HRS2,HRI3,HRS3
00189       DOUBLE PRECISION RR11,RR22,RR33,RR12,RR13,RR23
00190       DOUBLE PRECISION NF1,NF2,NF3,NF4,NF5,NF6
00191       DOUBLE PRECISION NG1,NG2,NG3,NG4,NG5,NG6
00192       DOUBLE PRECISION NH1,NH2,NH3,XS06,XS2880
00193 !
00194       INTEGER I1,I2,I3,I4,I5,I6,IELEM,II4,II5,II6,NPOU0
00195 !
00196       DOUBLE PRECISION EPSI
00197       DATA EPSI/1.D-4/
00198 !
00199 !-----------------------------------------------------------------------
00200 !
00201       XS06=XMUL/6.D0
00202       XS2880=XMUL/2880.D0
00203 !
00204 !-----------------------------------------------------------------------
00205 !
00206 !     VERY IMPORTANT !!!!!!!!!
00207 !
00208 !     SH%TYPR CHECKED TO CANCEL DIFFUSION ALONG Z
00209 !
00210 !     NOTE: XS06 ONLY USED WITH DIFFUSION ALONG Z !!!!!!!!!
00211 !
00212       IF(SH%TYPR.EQ.'0') THEN
00213         XS06=0.D0
00214       ENDIF
00215 !
00216 !     COULD BE OPTIMISED MORE BUT WOULD REQUEST SPLITTING LOOPS
00217 !     INTO HORIZONTAL AND VERTICAL
00218 !
00219 !-----------------------------------------------------------------------
00220 !
00221       IF(SF%ELM.NE.41) THEN
00222         IF (LNG.EQ.1) WRITE(LU,1000) SF%ELM
00223         IF (LNG.EQ.2) WRITE(LU,1001) SF%ELM
00224 1000    FORMAT(1X,'MT02PP (BIEF) : TYPE DE F NON PREVU : ',I6)
00225 1001    FORMAT(1X,'MT02PP (BIEF) : TYPE OF F NOT IMPLEMENTED: ',I6)
00226         CALL PLANTE(1)
00227         STOP
00228       ENDIF
00229       IF(SG%ELM.NE.41) THEN
00230         IF (LNG.EQ.1) WRITE(LU,2000) SG%ELM
00231         IF (LNG.EQ.2) WRITE(LU,2001) SG%ELM
00232 2000    FORMAT(1X,'MT02PP (BIEF) : TYPE DE G NON PREVU : ',I6)
00233 2001    FORMAT(1X,'MT02PP (BIEF) : TYPE OF G NOT IMPLEMENTED: ',I6)
00234         CALL PLANTE(1)
00235         STOP
00236       ENDIF
00237       IF(SH%ELM.NE.41) THEN
00238         IF (LNG.EQ.1) WRITE(LU,3000) SH%ELM
00239         IF (LNG.EQ.2) WRITE(LU,3001) SH%ELM
00240 3000    FORMAT(1X,'MT02PP (BIEF) : TYPE DE H NON PREVU : ',I6)
00241 3001    FORMAT(1X,'MT02PP (BIEF) : TYPE OF H NOT IMPLEMENTED: ',I6)
00242         CALL PLANTE(1)
00243         STOP
00244       ENDIF
00245 !
00246 !     SEE VISCLM OF TELEMAC-3D
00247 !     FOR THE TREATMENT OF P0 VERTICAL VISCOSITY ON THE VERTICAL
00248 !
00249       IF(SH%DIMDISC.EQ.0) THEN
00250 !       P1 VERTICAL VISCOSITY
00251         NPOU0=SH%DIM1/NPLAN
00252       ELSEIF(SH%DIMDISC.EQ.4111) THEN
00253 !       P0 VERTICAL VISCOSITY ON THE VERTICAL (SEE II4,5,6)
00254         NPOU0=0
00255       ELSE
00256         IF (LNG.EQ.1) WRITE(LU,4000) SH%DIMDISC
00257         IF (LNG.EQ.2) WRITE(LU,4001) SH%DIMDISC
00258 4000    FORMAT(1X,'MT02PP (BIEF) : DIMDISC DE H NON PREVU : ',I6)
00259 4001    FORMAT(1X,'MT02PP (BIEF): DIMDISC OF H NOT IMPLEMENTED: ',I6)
00260         CALL PLANTE(1)
00261         STOP
00262       ENDIF
00263 !
00264 !     VERSION WITH TREATMENT OF HYDROSTATIC INCONSISTENCIES
00265 !
00266       IF(INCHYD) THEN
00267 !
00268 !-----------------------------------------------------------------------
00269 !
00270 !     DIFFUSION ALONG X
00271 !
00272 !-----------------------------------------------------------------------
00273 !
00274 !   LOOP ON THE ELEMENTS
00275 !
00276       DO IELEM=1,NELEM
00277 !
00278         I1=IKLE(IELEM,1)
00279         I2=IKLE(IELEM,2)
00280         I3=IKLE(IELEM,3)
00281         I4=IKLE(IELEM,4)
00282         I5=IKLE(IELEM,5)
00283         I6=IKLE(IELEM,6)
00284 !       DEPENDING ON NPOU0, II4 WILL BE I4 OR I1, ETC
00285         II4=I1+NPOU0
00286         II5=I2+NPOU0
00287         II6=I3+NPOU0
00288 !
00289         H1=Z(I4)-Z(I1)
00290         H2=Z(I5)-Z(I2)
00291         H3=Z(I6)-Z(I3)
00292 !
00293         SH1=H1+H2+H3
00294         SHH=H1*H1+H2*H2+H3*H3
00295 !
00296         D1=Y(IELEM,2)-Y(IELEM,3)
00297 !       D2=Y(IELEM,3)-Y(IELEM,1)
00298         D2=Y(IELEM,3)
00299 !       D3=Y(IELEM,1)-Y(IELEM,2)
00300         D3=          -Y(IELEM,2)
00301 !
00302         D=XS2880/SURFAC(IELEM)
00303 !
00304         RI=-(Z(I1)*D1+Z(I2)*D2+Z(I3)*D3)
00305         RS=-(Z(I4)*D1+Z(I5)*D2+Z(I6)*D3)
00306         R=RI+RS
00307         RRI=R+RI+RI
00308         RRS=R+RS+RS
00309         RRRI=R*R+2*RI*RI
00310         RRRS=R*R+2*RS*RS
00311 !
00312         D12=D1*D2
00313         D13=D1*D3
00314         D23=D2*D3
00315 !
00316         IF(MAX(Z(I1),Z(I2),Z(I3)).GT.MIN(Z(I4),Z(I5),Z(I6)).OR.
00317      &     H1.LT.EPSI.OR.H2.LT.EPSI.OR.H3.LT.EPSI ) THEN
00318           NF1=0.D0
00319           NF2=0.D0
00320           NF3=0.D0
00321           NF4=0.D0
00322           NF5=0.D0
00323           NF6=0.D0
00324           NG1=0.D0
00325           NG2=0.D0
00326           NG3=0.D0
00327           NG4=0.D0
00328           NG5=0.D0
00329           NG6=0.D0
00330           NH1=0.D0
00331           NH2=0.D0
00332           NH3=0.D0
00333         ELSE
00334           NF1=F(I1)/H1
00335           NF2=F(I2)/H2
00336           NF3=F(I3)/H3
00337           NF4=F(I4)/H1
00338           NF5=F(I5)/H2
00339           NF6=F(I6)/H3
00340           NG1=G(I1)/H1
00341           NG2=G(I2)/H2
00342           NG3=G(I3)/H3
00343           NG4=G(I4)/H1
00344           NG5=G(I5)/H2
00345           NG6=G(I6)/H3
00346 !         DEPENDING ON THE CASE (II4=I1 OR I4, ETC.)
00347 !         ALTERNATIVE WITH VERTICAL LINEAR VISCOSITY (II4=I4)
00348 !         ALTERNATIVE WITH P0 VERTICAL VISCOSITY ON THE VERTICAL (II4=I1)
00349           NH1=(H(I1)+H(II4))/H1
00350           NH2=(H(I2)+H(II5))/H2
00351           NH3=(H(I3)+H(II6))/H3
00352         ENDIF
00353 !
00354         SNI=NF1+NF2+NF3
00355         SNS=NF4+NF5+NF6
00356         SNHI=NF1*H1+NF2*H2+NF3*H3
00357         SNHS=NF4*H1+NF5*H2+NF6*H3
00358         SNHHI=(SNI*SH1+SNHI+SNHI)*SH1+SNI*SHH
00359      &       +2*(NF1*H1*H1+NF2*H2*H2+NF3*H3*H3)
00360         SNHHS=(SNS*SH1+SNHS+SNHS)*SH1+SNS*SHH
00361      &       +2*(NF4*H1*H1+NF5*H2*H2+NF6*H3*H3)
00362         SNHH=SNHHI+SNHHS
00363         SNHHI=SNHH+SNHHI+SNHHI
00364         SNHHS=SNHH+SNHHS+SNHHS
00365 !
00366         HHI12=SNHHI*D12
00367         HHI13=SNHHI*D13
00368         HHI23=SNHHI*D23
00369         HHS12=SNHHS*D12
00370         HHS13=SNHHS*D13
00371         HHS23=SNHHS*D23
00372         HH12=SNHH*D12
00373         HH13=SNHH*D13
00374         HH23=SNHH*D23
00375 !
00376         SNHI1=(SNI+NF1)*(SH1+H1)+SNHI+NF1*H1
00377         SNHS1=(SNS+NF4)*(SH1+H1)+SNHS+NF4*H1
00378         SNHI2=(SNI+NF2)*(SH1+H2)+SNHI+NF2*H2
00379         SNHS2=(SNS+NF5)*(SH1+H2)+SNHS+NF5*H2
00380         SNHI3=(SNI+NF3)*(SH1+H3)+SNHI+NF3*H3
00381         SNHS3=(SNS+NF6)*(SH1+H3)+SNHS+NF6*H3
00382 !
00383         HRI1=RRI*SNHI1+R*SNHS1
00384         HRS1=RRS*SNHS1+R*SNHI1
00385         HRI2=RRI*SNHI2+R*SNHS2
00386         HRS2=RRS*SNHS2+R*SNHI2
00387         HRI3=RRI*SNHI3+R*SNHS3
00388         HRS3=RRS*SNHS3+R*SNHI3
00389 !
00390         RR11=2*(RRRI*(SNI+NF1+NF1)+RRRS*(SNS+NF4+NF4))
00391         RR22=2*(RRRI*(SNI+NF2+NF2)+RRRS*(SNS+NF5+NF5))
00392         RR33=2*(RRRI*(SNI+NF3+NF3)+RRRS*(SNS+NF6+NF6))
00393         RR12=   RRRI*(SNI+NF1+NF2)+RRRS*(SNS+NF4+NF5)
00394         RR13=   RRRI*(SNI+NF1+NF3)+RRRS*(SNS+NF4+NF6)
00395         RR23=   RRRI*(SNI+NF2+NF3)+RRRS*(SNS+NF5+NF6)
00396 !
00397 !       EXTRA-DIAGONAL TERMS
00398 !
00399         XM(IELEM, 1)=D*( HHI12  -D1*HRI2-D2*HRI1 +RR12)
00400         XM(IELEM, 2)=D*( HHI13  -D1*HRI3-D3*HRI1 +RR13)
00401         XM(IELEM, 3)=D*(-HH12-HH13+D1*(HRI1-HRS1)-RR11)
00402         XM(IELEM, 4)=D*( HH12   +D1*HRI2-D2*HRS1 -RR12)
00403         XM(IELEM, 5)=D*( HH13   +D1*HRI3-D3*HRS1 -RR13)
00404         XM(IELEM, 6)=D*( HHI23  -D2*HRI3-D3*HRI2 +RR23)
00405         XM(IELEM, 7)=D*( HH12   +D2*HRI1-D1*HRS2 -RR12)
00406         XM(IELEM, 8)=D*(-HH12-HH23+D2*(HRI2-HRS2)-RR22)
00407         XM(IELEM, 9)=D*( HH23   +D2*HRI3-D3*HRS2 -RR23)
00408         XM(IELEM,10)=D*( HH13   +D3*HRI1-D1*HRS3 -RR13)
00409         XM(IELEM,11)=D*( HH23   +D3*HRI2-D2*HRS3 -RR23)
00410         XM(IELEM,12)=D*(-HH13-HH23+D3*(HRI3-HRS3)-RR33)
00411         XM(IELEM,13)=D*( HHS12  +D1*HRS2+D2*HRS1 +RR12)
00412         XM(IELEM,14)=D*( HHS13  +D1*HRS3+D3*HRS1 +RR13)
00413         XM(IELEM,15)=D*( HHS23  +D2*HRS3+D3*HRS2 +RR23)
00414 !
00415 !-----------------------------------------------------------------------
00416 !
00417 !        DIFFUSION ALONG Y
00418 !
00419 !-----------------------------------------------------------------------
00420 !
00421         D1=X(IELEM,3)-X(IELEM,2)
00422 !       D2=X(IELEM,1)-X(IELEM,3)
00423         D2=          -X(IELEM,3)
00424 !       D3=X(IELEM,2)-X(IELEM,1)
00425         D3=X(IELEM,2)
00426 !
00427         RI=-(Z(I1)*D1+Z(I2)*D2+Z(I3)*D3)
00428         RS=-(Z(I4)*D1+Z(I5)*D2+Z(I6)*D3)
00429         R=RI+RS
00430         RRI=R+RI+RI
00431         RRS=R+RS+RS
00432         RRRI=R*R+2*RI*RI
00433         RRRS=R*R+2*RS*RS
00434 !
00435         D12=D1*D2
00436         D13=D1*D3
00437         D23=D2*D3
00438 !
00439         SNI=NG1+NG2+NG3
00440         SNS=NG4+NG5+NG6
00441         SNHI=NG1*H1+NG2*H2+NG3*H3
00442         SNHS=NG4*H1+NG5*H2+NG6*H3
00443         SNHHI=(SNI*SH1+SNHI+SNHI)*SH1+SNI*SHH
00444      &       +2*(NG1*H1*H1+NG2*H2*H2+NG3*H3*H3)
00445         SNHHS=(SNS*SH1+SNHS+SNHS)*SH1+SNS*SHH
00446      &       +2*(NG4*H1*H1+NG5*H2*H2+NG6*H3*H3)
00447         SNHH=SNHHI+SNHHS
00448         SNHHI=SNHH+SNHHI+SNHHI
00449         SNHHS=SNHH+SNHHS+SNHHS
00450 !
00451         HHI12=SNHHI*D12
00452         HHI13=SNHHI*D13
00453         HHI23=SNHHI*D23
00454         HHS12=SNHHS*D12
00455         HHS13=SNHHS*D13
00456         HHS23=SNHHS*D23
00457         HH12=SNHH*D12
00458         HH13=SNHH*D13
00459         HH23=SNHH*D23
00460 !
00461         SNHI1=(SNI+NG1)*(SH1+H1)+SNHI+NG1*H1
00462         SNHS1=(SNS+NG4)*(SH1+H1)+SNHS+NG4*H1
00463         SNHI2=(SNI+NG2)*(SH1+H2)+SNHI+NG2*H2
00464         SNHS2=(SNS+NG5)*(SH1+H2)+SNHS+NG5*H2
00465         SNHI3=(SNI+NG3)*(SH1+H3)+SNHI+NG3*H3
00466         SNHS3=(SNS+NG6)*(SH1+H3)+SNHS+NG6*H3
00467 !
00468         HRI1=RRI*SNHI1+R*SNHS1
00469         HRS1=RRS*SNHS1+R*SNHI1
00470         HRI2=RRI*SNHI2+R*SNHS2
00471         HRS2=RRS*SNHS2+R*SNHI2
00472         HRI3=RRI*SNHI3+R*SNHS3
00473         HRS3=RRS*SNHS3+R*SNHI3
00474 !
00475         RR11=2*(RRRI*(SNI+NG1+NG1)+RRRS*(SNS+NG4+NG4))
00476         RR22=2*(RRRI*(SNI+NG2+NG2)+RRRS*(SNS+NG5+NG5))
00477         RR33=2*(RRRI*(SNI+NG3+NG3)+RRRS*(SNS+NG6+NG6))
00478         RR12=   RRRI*(SNI+NG1+NG2)+RRRS*(SNS+NG4+NG5)
00479         RR13=   RRRI*(SNI+NG1+NG3)+RRRS*(SNS+NG4+NG6)
00480         RR23=   RRRI*(SNI+NG2+NG3)+RRRS*(SNS+NG5+NG6)
00481 !
00482 !       EXTRA-DIAGONAL TERMS
00483 !
00484         XM(IELEM, 1)=XM(IELEM, 1)+D*( HHI12  -D1*HRI2-D2*HRI1 +RR12)
00485         XM(IELEM, 2)=XM(IELEM, 2)+D*( HHI13  -D1*HRI3-D3*HRI1 +RR13)
00486         XM(IELEM, 3)=XM(IELEM, 3)+D*(-HH12-HH13+D1*(HRI1-HRS1)-RR11)
00487         XM(IELEM, 4)=XM(IELEM, 4)+D*( HH12   +D1*HRI2-D2*HRS1 -RR12)
00488         XM(IELEM, 5)=XM(IELEM, 5)+D*( HH13   +D1*HRI3-D3*HRS1 -RR13)
00489         XM(IELEM, 6)=XM(IELEM, 6)+D*( HHI23  -D2*HRI3-D3*HRI2 +RR23)
00490         XM(IELEM, 7)=XM(IELEM, 7)+D*( HH12   +D2*HRI1-D1*HRS2 -RR12)
00491         XM(IELEM, 8)=XM(IELEM, 8)+D*(-HH12-HH23+D2*(HRI2-HRS2)-RR22)
00492         XM(IELEM, 9)=XM(IELEM, 9)+D*( HH23   +D2*HRI3-D3*HRS2 -RR23)
00493         XM(IELEM,10)=XM(IELEM,10)+D*( HH13   +D3*HRI1-D1*HRS3 -RR13)
00494         XM(IELEM,11)=XM(IELEM,11)+D*( HH23   +D3*HRI2-D2*HRS3 -RR23)
00495         XM(IELEM,12)=XM(IELEM,12)+D*(-HH13-HH23+D3*(HRI3-HRS3)-RR33)
00496         XM(IELEM,13)=XM(IELEM,13)+D*( HHS12  +D1*HRS2+D2*HRS1 +RR12)
00497         XM(IELEM,14)=XM(IELEM,14)+D*( HHS13  +D1*HRS3+D3*HRS1 +RR13)
00498         XM(IELEM,15)=XM(IELEM,15)+D*( HHS23  +D2*HRS3+D3*HRS2 +RR23)
00499 !
00500 !-----------------------------------------------------------------------
00501 !
00502 !        DIFFUSION ALONG Z
00503 !
00504 !-----------------------------------------------------------------------
00505 !
00506 !       VERSION WITH SIMPLIFICATIONS TO ACHIEVE MONOTONY OF THE MATRIX
00507 !
00508         D=SURFAC(IELEM)*XS06
00509 !
00510 !       EXTRA-DIAGONAL TERMS
00511 !
00512         XM(IELEM, 3)=XM(IELEM, 3)-D*NH1
00513         XM(IELEM, 8)=XM(IELEM, 8)-D*NH2
00514         XM(IELEM,12)=XM(IELEM,12)-D*NH3
00515 !
00516 !-----------------------------------------------------------------------
00517 !
00518 !       OLD VERSION OF DIFFUSION ALONG Z
00519 !
00520 !       R=NH1+NH2+NH3+NH4+NH5+NH6
00521 !       D=((X(I2)-X(I1))*(Y(I3)-Y(I1))-(X(I3)-X(I1))*(Y(I2)-Y(I1)))
00522 !    *   *XMUL/240.D0
00523 !
00524 !       RR11=(R+NH1+NH1+NH4+NH4)*(D+D)
00525 !       RR22=(R+NH2+NH2+NH5+NH5)*(D+D)
00526 !       RR33=(R+NH3+NH3+NH6+NH6)*(D+D)
00527 !       RR12=(R+NH1+NH2+NH4+NH5)*D
00528 !       RR13=(R+NH1+NH3+NH4+NH6)*D
00529 !       RR23=(R+NH2+NH3+NH5+NH6)*D
00530 !
00531 !       EXTRA-DIAGONAL TERMS
00532 !
00533 !       XM(IELEM, 1)=XM(IELEM, 1)+RR12
00534 !       XM(IELEM, 2)=XM(IELEM, 2)+RR13
00535 !       XM(IELEM, 3)=XM(IELEM, 3)-RR11
00536 !       XM(IELEM, 4)=XM(IELEM, 4)-RR12
00537 !       XM(IELEM, 5)=XM(IELEM, 5)-RR13
00538 !       XM(IELEM, 6)=XM(IELEM, 6)+RR23
00539 !       XM(IELEM, 7)=XM(IELEM, 7)-RR12
00540 !       XM(IELEM, 8)=XM(IELEM, 8)-RR22
00541 !       XM(IELEM, 9)=XM(IELEM, 9)-RR23
00542 !       XM(IELEM,10)=XM(IELEM,10)-RR13
00543 !       XM(IELEM,11)=XM(IELEM,11)-RR23
00544 !       XM(IELEM,12)=XM(IELEM,12)-RR33
00545 !       XM(IELEM,13)=XM(IELEM,13)+RR12
00546 !       XM(IELEM,14)=XM(IELEM,14)+RR13
00547 !       XM(IELEM,15)=XM(IELEM,15)+RR23
00548 !
00549 !
00550 !-----------------------------------------------------------------------
00551 !
00552       ENDDO
00553 !
00554 !-----------------------------------------------------------------------
00555 !
00556       ELSE
00557 !
00558 !     VERSION WITHOUT TREATMENT OF HYDROSTATIC INCONSISTENCIES
00559 !
00560 !-----------------------------------------------------------------------
00561 !
00562 !     DIFFUSION ALONG X
00563 !
00564 !-----------------------------------------------------------------------
00565 !
00566 !   LOOP ON THE ELEMENTS
00567 !
00568       DO IELEM=1,NELEM
00569 !
00570         I1=IKLE(IELEM,1)
00571         I2=IKLE(IELEM,2)
00572         I3=IKLE(IELEM,3)
00573         I4=IKLE(IELEM,4)
00574         I5=IKLE(IELEM,5)
00575         I6=IKLE(IELEM,6)
00576 !       DEPENDING ON NPOU0, II4 WILL BE I4 OR I1, ETC
00577         II4=I1+NPOU0
00578         II5=I2+NPOU0
00579         II6=I3+NPOU0
00580 !
00581         H1=Z(I4)-Z(I1)
00582         H2=Z(I5)-Z(I2)
00583         H3=Z(I6)-Z(I3)
00584 !
00585         SH1=H1+H2+H3
00586         SHH=H1*H1+H2*H2+H3*H3
00587 !
00588         D1=Y(IELEM,2)-Y(IELEM,3)
00589 !       D2=Y(IELEM,3)-Y(IELEM,1)
00590         D2=Y(IELEM,3)
00591 !       D3=Y(IELEM,1)-Y(IELEM,2)
00592         D3=          -Y(IELEM,2)
00593 !
00594         D=XS2880/SURFAC(IELEM)
00595 !
00596         RI=-(Z(I1)*D1+Z(I2)*D2+Z(I3)*D3)
00597         RS=-(Z(I4)*D1+Z(I5)*D2+Z(I6)*D3)
00598         R=RI+RS
00599         RRI=R+RI+RI
00600         RRS=R+RS+RS
00601         RRRI=R*R+2*RI*RI
00602         RRRS=R*R+2*RS*RS
00603 !
00604         D12=D1*D2
00605         D13=D1*D3
00606         D23=D2*D3
00607 !
00608         IF(H1.LT.EPSI.OR.H2.LT.EPSI.OR.H3.LT.EPSI) THEN
00609           NF1=0.D0
00610           NF2=0.D0
00611           NF3=0.D0
00612           NF4=0.D0
00613           NF5=0.D0
00614           NF6=0.D0
00615           NG1=0.D0
00616           NG2=0.D0
00617           NG3=0.D0
00618           NG4=0.D0
00619           NG5=0.D0
00620           NG6=0.D0
00621           NH1=0.D0
00622           NH2=0.D0
00623           NH3=0.D0
00624         ELSE
00625           NF1=F(I1)/H1
00626           NF2=F(I2)/H2
00627           NF3=F(I3)/H3
00628           NF4=F(I4)/H1
00629           NF5=F(I5)/H2
00630           NF6=F(I6)/H3
00631           NG1=G(I1)/H1
00632           NG2=G(I2)/H2
00633           NG3=G(I3)/H3
00634           NG4=G(I4)/H1
00635           NG5=G(I5)/H2
00636           NG6=G(I6)/H3
00637 !         DEPENDING ON THE CASE (II4=I1 OR I4, ETC.)
00638 !         ALTERNATIVE WITH VERTICAL LINEAR VISCOSITY (II4=I4)
00639 !         ALTERNATIVE WITH P0 VERTICAL VISCOSITY ON THE VERTICAL (II4=I1)
00640           NH1=(H(I1)+H(II4))/H1
00641           NH2=(H(I2)+H(II5))/H2
00642           NH3=(H(I3)+H(II6))/H3
00643         ENDIF
00644 !
00645         SNI=NF1+NF2+NF3
00646         SNS=NF4+NF5+NF6
00647         SNHI=NF1*H1+NF2*H2+NF3*H3
00648         SNHS=NF4*H1+NF5*H2+NF6*H3
00649         SNHHI=(SNI*SH1+SNHI+SNHI)*SH1+SNI*SHH
00650      &       +2*(NF1*H1*H1+NF2*H2*H2+NF3*H3*H3)
00651         SNHHS=(SNS*SH1+SNHS+SNHS)*SH1+SNS*SHH
00652      &       +2*(NF4*H1*H1+NF5*H2*H2+NF6*H3*H3)
00653         SNHH=SNHHI+SNHHS
00654         SNHHI=SNHH+SNHHI+SNHHI
00655         SNHHS=SNHH+SNHHS+SNHHS
00656 !
00657         HHI12=SNHHI*D12
00658         HHI13=SNHHI*D13
00659         HHI23=SNHHI*D23
00660         HHS12=SNHHS*D12
00661         HHS13=SNHHS*D13
00662         HHS23=SNHHS*D23
00663         HH12=SNHH*D12
00664         HH13=SNHH*D13
00665         HH23=SNHH*D23
00666 !
00667         SNHI1=(SNI+NF1)*(SH1+H1)+SNHI+NF1*H1
00668         SNHS1=(SNS+NF4)*(SH1+H1)+SNHS+NF4*H1
00669         SNHI2=(SNI+NF2)*(SH1+H2)+SNHI+NF2*H2
00670         SNHS2=(SNS+NF5)*(SH1+H2)+SNHS+NF5*H2
00671         SNHI3=(SNI+NF3)*(SH1+H3)+SNHI+NF3*H3
00672         SNHS3=(SNS+NF6)*(SH1+H3)+SNHS+NF6*H3
00673 !
00674         HRI1=RRI*SNHI1+R*SNHS1
00675         HRS1=RRS*SNHS1+R*SNHI1
00676         HRI2=RRI*SNHI2+R*SNHS2
00677         HRS2=RRS*SNHS2+R*SNHI2
00678         HRI3=RRI*SNHI3+R*SNHS3
00679         HRS3=RRS*SNHS3+R*SNHI3
00680 !
00681         RR11=2*(RRRI*(SNI+NF1+NF1)+RRRS*(SNS+NF4+NF4))
00682         RR22=2*(RRRI*(SNI+NF2+NF2)+RRRS*(SNS+NF5+NF5))
00683         RR33=2*(RRRI*(SNI+NF3+NF3)+RRRS*(SNS+NF6+NF6))
00684         RR12=   RRRI*(SNI+NF1+NF2)+RRRS*(SNS+NF4+NF5)
00685         RR13=   RRRI*(SNI+NF1+NF3)+RRRS*(SNS+NF4+NF6)
00686         RR23=   RRRI*(SNI+NF2+NF3)+RRRS*(SNS+NF5+NF6)
00687 !
00688 !       EXTRA-DIAGONAL TERMS
00689 !
00690         XM(IELEM, 1)=D*( HHI12  -D1*HRI2-D2*HRI1 +RR12)
00691         XM(IELEM, 2)=D*( HHI13  -D1*HRI3-D3*HRI1 +RR13)
00692         XM(IELEM, 3)=D*(-HH12-HH13+D1*(HRI1-HRS1)-RR11)
00693         XM(IELEM, 4)=D*( HH12   +D1*HRI2-D2*HRS1 -RR12)
00694         XM(IELEM, 5)=D*( HH13   +D1*HRI3-D3*HRS1 -RR13)
00695         XM(IELEM, 6)=D*( HHI23  -D2*HRI3-D3*HRI2 +RR23)
00696         XM(IELEM, 7)=D*( HH12   +D2*HRI1-D1*HRS2 -RR12)
00697         XM(IELEM, 8)=D*(-HH12-HH23+D2*(HRI2-HRS2)-RR22)
00698         XM(IELEM, 9)=D*( HH23   +D2*HRI3-D3*HRS2 -RR23)
00699         XM(IELEM,10)=D*( HH13   +D3*HRI1-D1*HRS3 -RR13)
00700         XM(IELEM,11)=D*( HH23   +D3*HRI2-D2*HRS3 -RR23)
00701         XM(IELEM,12)=D*(-HH13-HH23+D3*(HRI3-HRS3)-RR33)
00702         XM(IELEM,13)=D*( HHS12  +D1*HRS2+D2*HRS1 +RR12)
00703         XM(IELEM,14)=D*( HHS13  +D1*HRS3+D3*HRS1 +RR13)
00704         XM(IELEM,15)=D*( HHS23  +D2*HRS3+D3*HRS2 +RR23)
00705 !
00706 !-----------------------------------------------------------------------
00707 !
00708 !        DIFFUSION ALONG Y
00709 !
00710 !-----------------------------------------------------------------------
00711 !
00712         D1=X(IELEM,3)-X(IELEM,2)
00713 !       D2=X(IELEM,1)-X(IELEM,3)
00714         D2=          -X(IELEM,3)
00715 !       D3=X(IELEM,2)-X(IELEM,1)
00716         D3=X(IELEM,2)
00717 !
00718         RI=-(Z(I1)*D1+Z(I2)*D2+Z(I3)*D3)
00719         RS=-(Z(I4)*D1+Z(I5)*D2+Z(I6)*D3)
00720         R=RI+RS
00721         RRI=R+RI+RI
00722         RRS=R+RS+RS
00723         RRRI=R*R+2*RI*RI
00724         RRRS=R*R+2*RS*RS
00725 !
00726         D12=D1*D2
00727         D13=D1*D3
00728         D23=D2*D3
00729 !
00730         SNI=NG1+NG2+NG3
00731         SNS=NG4+NG5+NG6
00732         SNHI=NG1*H1+NG2*H2+NG3*H3
00733         SNHS=NG4*H1+NG5*H2+NG6*H3
00734         SNHHI=(SNI*SH1+SNHI+SNHI)*SH1+SNI*SHH
00735      &       +2*(NG1*H1*H1+NG2*H2*H2+NG3*H3*H3)
00736         SNHHS=(SNS*SH1+SNHS+SNHS)*SH1+SNS*SHH
00737      &       +2*(NG4*H1*H1+NG5*H2*H2+NG6*H3*H3)
00738         SNHH=SNHHI+SNHHS
00739         SNHHI=SNHH+SNHHI+SNHHI
00740         SNHHS=SNHH+SNHHS+SNHHS
00741 !
00742         HHI12=SNHHI*D12
00743         HHI13=SNHHI*D13
00744         HHI23=SNHHI*D23
00745         HHS12=SNHHS*D12
00746         HHS13=SNHHS*D13
00747         HHS23=SNHHS*D23
00748         HH12=SNHH*D12
00749         HH13=SNHH*D13
00750         HH23=SNHH*D23
00751 !
00752         SNHI1=(SNI+NG1)*(SH1+H1)+SNHI+NG1*H1
00753         SNHS1=(SNS+NG4)*(SH1+H1)+SNHS+NG4*H1
00754         SNHI2=(SNI+NG2)*(SH1+H2)+SNHI+NG2*H2
00755         SNHS2=(SNS+NG5)*(SH1+H2)+SNHS+NG5*H2
00756         SNHI3=(SNI+NG3)*(SH1+H3)+SNHI+NG3*H3
00757         SNHS3=(SNS+NG6)*(SH1+H3)+SNHS+NG6*H3
00758 !
00759         HRI1=RRI*SNHI1+R*SNHS1
00760         HRS1=RRS*SNHS1+R*SNHI1
00761         HRI2=RRI*SNHI2+R*SNHS2
00762         HRS2=RRS*SNHS2+R*SNHI2
00763         HRI3=RRI*SNHI3+R*SNHS3
00764         HRS3=RRS*SNHS3+R*SNHI3
00765 !
00766         RR11=2*(RRRI*(SNI+NG1+NG1)+RRRS*(SNS+NG4+NG4))
00767         RR22=2*(RRRI*(SNI+NG2+NG2)+RRRS*(SNS+NG5+NG5))
00768         RR33=2*(RRRI*(SNI+NG3+NG3)+RRRS*(SNS+NG6+NG6))
00769         RR12=   RRRI*(SNI+NG1+NG2)+RRRS*(SNS+NG4+NG5)
00770         RR13=   RRRI*(SNI+NG1+NG3)+RRRS*(SNS+NG4+NG6)
00771         RR23=   RRRI*(SNI+NG2+NG3)+RRRS*(SNS+NG5+NG6)
00772 !
00773 !       EXTRA-DIAGONAL TERMS
00774 !
00775         XM(IELEM, 1)=XM(IELEM, 1)+D*( HHI12  -D1*HRI2-D2*HRI1 +RR12)
00776         XM(IELEM, 2)=XM(IELEM, 2)+D*( HHI13  -D1*HRI3-D3*HRI1 +RR13)
00777         XM(IELEM, 3)=XM(IELEM, 3)+D*(-HH12-HH13+D1*(HRI1-HRS1)-RR11)
00778         XM(IELEM, 4)=XM(IELEM, 4)+D*( HH12   +D1*HRI2-D2*HRS1 -RR12)
00779         XM(IELEM, 5)=XM(IELEM, 5)+D*( HH13   +D1*HRI3-D3*HRS1 -RR13)
00780         XM(IELEM, 6)=XM(IELEM, 6)+D*( HHI23  -D2*HRI3-D3*HRI2 +RR23)
00781         XM(IELEM, 7)=XM(IELEM, 7)+D*( HH12   +D2*HRI1-D1*HRS2 -RR12)
00782         XM(IELEM, 8)=XM(IELEM, 8)+D*(-HH12-HH23+D2*(HRI2-HRS2)-RR22)
00783         XM(IELEM, 9)=XM(IELEM, 9)+D*( HH23   +D2*HRI3-D3*HRS2 -RR23)
00784         XM(IELEM,10)=XM(IELEM,10)+D*( HH13   +D3*HRI1-D1*HRS3 -RR13)
00785         XM(IELEM,11)=XM(IELEM,11)+D*( HH23   +D3*HRI2-D2*HRS3 -RR23)
00786         XM(IELEM,12)=XM(IELEM,12)+D*(-HH13-HH23+D3*(HRI3-HRS3)-RR33)
00787         XM(IELEM,13)=XM(IELEM,13)+D*( HHS12  +D1*HRS2+D2*HRS1 +RR12)
00788         XM(IELEM,14)=XM(IELEM,14)+D*( HHS13  +D1*HRS3+D3*HRS1 +RR13)
00789         XM(IELEM,15)=XM(IELEM,15)+D*( HHS23  +D2*HRS3+D3*HRS2 +RR23)
00790 !
00791 !-----------------------------------------------------------------------
00792 !
00793 !       DIFFUSION ALONG Z
00794 !
00795 !-----------------------------------------------------------------------
00796 !
00797 !       VERSION WITH SIMPLIFICATIONS TO ACHIEVE MONOTONY OF THE MATRIX
00798 !
00799         D=SURFAC(IELEM)*XS06
00800 !
00801 !       EXTRA-DIAGONAL TERMS
00802 !
00803         XM(IELEM, 3)=XM(IELEM, 3)-D*NH1
00804         XM(IELEM, 8)=XM(IELEM, 8)-D*NH2
00805         XM(IELEM,12)=XM(IELEM,12)-D*NH3
00806 !
00807 !-----------------------------------------------------------------------
00808 !
00809       ENDDO
00810 !
00811 !     IF(INCHYD) THEN
00812       ENDIF
00813 !
00814 !-----------------------------------------------------------------------
00815 !
00816 !     POSITIVE EXTRA-DIAGONAL TERMS REMOVED
00817 !     TO ENSURE MONOTONY
00818 !
00819       IF(FORMUL(14:16).EQ.'MON') THEN
00820 !
00821         IF(XMUL.GT.0.D0) THEN
00822           DO IELEM=1,NELEM
00823             XM(IELEM, 1)=MIN(XM(IELEM, 1),0.D0)
00824             XM(IELEM, 2)=MIN(XM(IELEM, 2),0.D0)
00825             XM(IELEM, 3)=MIN(XM(IELEM, 3),0.D0)
00826             XM(IELEM, 4)=MIN(XM(IELEM, 4),0.D0)
00827             XM(IELEM, 5)=MIN(XM(IELEM, 5),0.D0)
00828             XM(IELEM, 6)=MIN(XM(IELEM, 6),0.D0)
00829             XM(IELEM, 7)=MIN(XM(IELEM, 7),0.D0)
00830             XM(IELEM, 8)=MIN(XM(IELEM, 8),0.D0)
00831             XM(IELEM, 9)=MIN(XM(IELEM, 9),0.D0)
00832             XM(IELEM,10)=MIN(XM(IELEM,10),0.D0)
00833             XM(IELEM,11)=MIN(XM(IELEM,11),0.D0)
00834             XM(IELEM,12)=MIN(XM(IELEM,12),0.D0)
00835             XM(IELEM,13)=MIN(XM(IELEM,13),0.D0)
00836             XM(IELEM,14)=MIN(XM(IELEM,14),0.D0)
00837             XM(IELEM,15)=MIN(XM(IELEM,15),0.D0)
00838           ENDDO
00839         ELSE
00840           DO IELEM=1,NELEM
00841             XM(IELEM, 1)=MAX(XM(IELEM, 1),0.D0)
00842             XM(IELEM, 2)=MAX(XM(IELEM, 2),0.D0)
00843             XM(IELEM, 3)=MAX(XM(IELEM, 3),0.D0)
00844             XM(IELEM, 4)=MAX(XM(IELEM, 4),0.D0)
00845             XM(IELEM, 5)=MAX(XM(IELEM, 5),0.D0)
00846             XM(IELEM, 6)=MAX(XM(IELEM, 6),0.D0)
00847             XM(IELEM, 7)=MAX(XM(IELEM, 7),0.D0)
00848             XM(IELEM, 8)=MAX(XM(IELEM, 8),0.D0)
00849             XM(IELEM, 9)=MAX(XM(IELEM, 9),0.D0)
00850             XM(IELEM,10)=MAX(XM(IELEM,10),0.D0)
00851             XM(IELEM,11)=MAX(XM(IELEM,11),0.D0)
00852             XM(IELEM,12)=MAX(XM(IELEM,12),0.D0)
00853             XM(IELEM,13)=MAX(XM(IELEM,13),0.D0)
00854             XM(IELEM,14)=MAX(XM(IELEM,14),0.D0)
00855             XM(IELEM,15)=MAX(XM(IELEM,15),0.D0)
00856           ENDDO
00857         ENDIF
00858       ENDIF
00859 !
00860 !-----------------------------------------------------------------------
00861 !
00862 !     DIAGONAL TERMS OBTAINED GIVEN THAT :
00863 !     SUM OF THE TERMS IN A ROW =0
00864 !
00865 !-----------------------------------------------------------------------
00866 !
00867       DO IELEM=1,NELEM
00868 !
00869         T(IELEM,1)= -XM(IELEM,01)
00870      &              -XM(IELEM,02)
00871      &              -XM(IELEM,03)
00872      &              -XM(IELEM,04)
00873      &              -XM(IELEM,05)
00874         T(IELEM,2)= -XM(IELEM,01)
00875      &              -XM(IELEM,06)
00876      &              -XM(IELEM,07)
00877      &              -XM(IELEM,08)
00878      &              -XM(IELEM,09)
00879         T(IELEM,3)= -XM(IELEM,02)
00880      &              -XM(IELEM,06)
00881      &              -XM(IELEM,10)
00882      &              -XM(IELEM,11)
00883      &              -XM(IELEM,12)
00884         T(IELEM,4)= -XM(IELEM,03)
00885      &              -XM(IELEM,07)
00886      &              -XM(IELEM,10)
00887      &              -XM(IELEM,13)
00888      &              -XM(IELEM,14)
00889         T(IELEM,5)= -XM(IELEM,04)
00890      &              -XM(IELEM,08)
00891      &              -XM(IELEM,11)
00892      &              -XM(IELEM,13)
00893      &              -XM(IELEM,15)
00894         T(IELEM,6)= -XM(IELEM,05)
00895      &              -XM(IELEM,09)
00896      &              -XM(IELEM,12)
00897      &              -XM(IELEM,14)
00898      &              -XM(IELEM,15)
00899 !
00900       ENDDO
00901 !
00902 !-----------------------------------------------------------------------
00903 !
00904       RETURN
00905       END

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