flusrc.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\flusrc.f
00002 !
00046                      SUBROUTINE FLUSRC
00047 !                    *****************
00048 !
00049      &(IEL1,IEL2,ISEGIN,VNOIN,W,FLUSCE,X,Y,AIRS,NPOIN,NSEG,ZF,EPS,G)
00050 !
00051 !***********************************************************************
00052 ! TELEMAC2D   V6P1                                   21/08/2010
00053 !***********************************************************************
00054 !
00055 !
00056 !
00057 !
00058 !
00059 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00060 !| AIRS           |-->| AREA OF THE CELLS.
00061 !| EPS            |-->| TOLERANCE
00062 !| FLUSCE         |<->| SOURCE FLUXES
00063 !| G              |-->| GRAVITY
00064 !| IEL1           |-->| FIRST ELEMENT NUMBER
00065 !| IEL2           |-->| SECOND ELEMENT NUMBER
00066 !| ISEGIN         |-->| SEGMENT NUMBER
00067 !| NPOIN          |-->| TOTAL NUMBER OF NODES
00068 !| NSEG           |-->| TOTAL NUMBER OF SEGMENTS IN THE MESH
00069 !| VNOIN          |-->| NORMAL VECTOR TO THE INTERFACE
00070 !|                |   | (2 FIRST COMPONENTS) AND
00071 !|                |   | LENGTH OF THE SEGMENT (3RD COMPONENT)
00072 !| W              |-->| CONSERVATIVE VARIABLE OF THE PROBLEM AT TIME TN
00073 !| X              |-->| X COORDINATES
00074 !| Y              |-->| Y COORDINATES
00075 !| ZF             |-->| BATHYMETRY
00076 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00077 !
00078       IMPLICIT NONE
00079       INTEGER LNG,LU
00080       COMMON/INFO/LNG,LU
00081 !
00082 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00083 !
00084       INTEGER, INTENT(IN)             :: NPOIN,NSEG,ISEGIN,IEL1,IEL2
00085       DOUBLE PRECISION, INTENT(IN)    :: G,EPS,VNOIN(3,NSEG),ZF(NPOIN)
00086       DOUBLE PRECISION, INTENT(IN)    :: AIRS(NPOIN),X(NPOIN),Y(NPOIN)
00087       DOUBLE PRECISION, INTENT(IN)    :: W(3,NPOIN)
00088       DOUBLE PRECISION, INTENT(INOUT) :: FLUSCE(3,NPOIN)
00089 !
00090 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00091 !
00092       INTEGER INDIC(2)
00093 !
00094       DOUBLE PRECISION XGI, YGI, XGJ, YGJ, DIJ, A1, A2
00095       DOUBLE PRECISION HI,UI,VI,HJ,VJ,UJ,XN,YN
00096       DOUBLE PRECISION CT2,CT,RLAMB0,RLAMBM,ALPHA,CI2
00097       DOUBLE PRECISION CJ,CJ2,RLAMBJ,RLAMBI,RLAMBP
00098       DOUBLE PRECISION RI,RJ,UT,VT,CI,UN
00099       DOUBLE PRECISION T11(3),T21(3),TS11(3),TS21(3)
00100       DOUBLE PRECISION T12(3),T22(3),TS12(3),TS22(3)
00101       DOUBLE PRECISION GE(3),FE(3),ZF1,ZF2,PSA1,PSA2,PSA
00102 !
00103       INTRINSIC MIN,MAX
00104 !
00105 !-----------------------------------------------------------------------
00106 !
00107 !------
00108 ! 1. COMPUTES SOURCES TERMS AT THE INTERFACE IEL1 , IEL2
00109 !------
00110 !
00111 !
00112 !
00113 !   --->    SOME INTERMEDIATE CALCULATIONS
00114 !           ------------------------------
00115 !
00116       HI = W(1,IEL1)
00117       IF (HI.GT.EPS) THEN
00118         UI = W(2,IEL1) / HI
00119         VI = W(3,IEL1) / HI
00120         INDIC(1) = 0
00121       ELSE
00122         UI = 0.D0
00123         VI = 0.D0
00124         INDIC(1) = 1
00125       ENDIF
00126 !
00127       HJ = W(1,IEL2)
00128       IF ( HJ.GT.EPS) THEN
00129         UJ = W(2,IEL2) / HJ
00130         VJ = W(3,IEL2) / HJ
00131         INDIC(2) = 0
00132       ELSE
00133         UJ = 0.D0
00134         VJ = 0.D0
00135         INDIC(2) = 1
00136       ENDIF
00137 !
00138       XN = VNOIN (1,ISEGIN)
00139       YN = VNOIN (2,ISEGIN)
00140 !
00141 !    BOTTOM MODIFICATION: AT REST WITH A DRY ELEMENT
00142 !
00143       UN = UJ*XN+VJ*YN
00144       ZF1 = ZF(IEL1)
00145       IF(INDIC(1).EQ.1) THEN
00146         IF((ZF1+EPS.GT.ZF(IEL2)+HJ).AND.(UN.GE.-EPS)) THEN
00147           ZF1 = ZF(IEL2)+HJ-EPS
00148         ENDIF
00149       ENDIF
00150 !
00151       UN = UI*XN+VI*YN
00152       ZF2 = ZF(IEL2)
00153       IF(INDIC(2).EQ.1) THEN
00154         IF((ZF2+EPS.GT.ZF(IEL1)+HI).AND.(UN.LE.EPS)) THEN
00155           ZF2 = ZF(IEL1)+HI-EPS
00156         ENDIF
00157       ENDIF
00158 !
00159 !   --->    COMPUTES THE AVERAGES OF ROE OF U,V,H,C**2 AND C
00160 !           ---------------------------------------------
00161 !
00162       IF(HI.LE.0.D0) THEN
00163         HI = 0.D0
00164 ! FOLLOWING LINE INSERTED BY JMH
00165         RI = 0.D0
00166       ELSE
00167         RI = SQRT ( HI )
00168       ENDIF
00169       IF(HJ.LE.0.D0) THEN
00170         HJ = 0.D0
00171 ! FOLLOWING LINE INSERTED BY JMH
00172         RJ = 0.D0
00173       ELSE
00174         RJ = SQRT ( HJ )
00175       ENDIF
00176 !     MAX OF THE TWO FOLLOWING LINES INSERTED BY JMH
00177       UT = ( RI * UI + RJ * UJ ) / MAX(RI+RJ,1.D-8)
00178       VT = ( RI * VI + RJ * VJ ) / MAX(RI+RJ,1.D-8)
00179       CT2 = G*(HI+HJ)/2.D0
00180       CT = SQRT ( CT2 )
00181 !
00182 !   --->  TEST ON THE SIGN OF THE EIGENVALUE LAMB0 =
00183 !           ----------------------------------------------------------
00184 !
00185       RLAMB0 = UT * XN + VT * YN
00186 !
00187 !     COMPUTES EIGENVALUES MATRICES
00188 !--------------------------------------------
00189 !
00190       T11(1) = 1.D0
00191       T11(2) = UT - CT * XN
00192       T11(3) = VT - CT * YN
00193       T21(1) = 0.D0
00194       T21(2) = CT * YN
00195       T21(3) = -CT * XN
00196 !
00197       T12(1) = 1.D0
00198       T12(2) = UT + CT * XN
00199       T12(3) = VT + CT * YN
00200       T22(1) = 0.D0
00201       T22(2) = -CT * YN
00202       T22(3) = +CT * XN
00203 !
00204       TS11(1) = (UT * XN + VT * YN) * CT + CT2
00205       TS21(1) = (2.D0 * VT * XN - 2.D0 * UT * YN) * CT
00206       TS11(2) = -XN * CT
00207       TS21(2) = 2.D0 * YN * CT
00208       TS11(3) = -YN * CT
00209       TS21(3) = -2.D0 * XN * CT
00210 !
00211       TS12(1) = -(UT * XN + VT * YN) * CT + CT2
00212       TS22(1) = -(2.D0 * VT * XN - 2.D0 * UT * YN) * CT
00213       TS12(2) = +XN * CT
00214       TS22(2) = -2.D0 * YN * CT
00215       TS12(3) = +YN * CT
00216       TS22(3) = +2.D0* XN * CT
00217 !
00218 !----------CALCULS POUR LES TERMES SOURCES--------------------
00219 !
00220 !
00221       XGI = X(IEL1)
00222       YGI = Y(IEL1)
00223       XGJ = X(IEL2)
00224       YGJ = Y(IEL2)
00225 !
00226       DIJ = SQRT ((XGJ -XGI)**2 + (YGJ - YGI)**2)
00227       A1  = VNOIN(3,ISEGIN)*DIJ/2.D0
00228       A2  = VNOIN(3,ISEGIN)*DIJ/2.D0
00229 !
00230 !  BOTTOM GRADIENTS
00231 !
00232       GE(1)=0.D0
00233       GE(2)=G*((HI+HJ)/2.D0)*(ZF2-ZF1)*XN/DIJ
00234       GE(3)=G*((HI+HJ)/2.D0)*(ZF2-ZF1)*YN/DIJ
00235 !
00236 !  FRICTION TERMS
00237 !
00238 !     CH = 900.D0
00239 !     H = (CT2/G)**(1./3.)
00240       FE(1)= 0.D0
00241 !     FE(2)= G*UT*SQRT(UT**2 + VT**2)/((CH**2)*H)
00242 !     FE(3)= G*VT*SQRT(UT**2 + VT**2)/((CH**2)*H)
00243       FE(2) = 0.D0
00244       FE(3) = 0.D0
00245 !
00246       GE(1)=GE(1)+FE(1)
00247       GE(2)=GE(2)+FE(2)
00248       GE(3)=GE(3)+FE(3)
00249 !
00250 !---------------------------------------------------------------------
00251       IF(RLAMB0.GE.-.000001D0) THEN
00252 !     ---- EXIT SEGMENT ---------
00253 !
00254 !   --->   SMALL CALCULATIONS
00255 !
00256         RLAMBM = RLAMB0 - CT
00257 !
00258 !
00259         ALPHA = UI * XN + VI * YN
00260 !
00261 !TBTB BEGINNING: MODIFICATION OF RLAMBM IF RLAMBM
00262 !
00263         CI2 = G*HI
00264         CI = SQRT (CI2)
00265         CJ2 =  G*HJ
00266         CJ = SQRT (CJ2)
00267         RLAMBI = ALPHA - CI
00268         RLAMBJ = UJ * XN + VJ * YN - CJ
00269 !
00270         IF ( RLAMBI .LT. 0.D0 .AND. RLAMBJ .GT. 0.D0
00271      &                                                ) THEN
00272           RLAMBM = MIN(0.D0,RLAMBM) - ABS(RLAMBI - RLAMBJ) / 4.D0
00273         ENDIF
00274 !
00275 !------------CALCUL DES TERMES SOURCES ------------------------
00276 !
00277         FLUSCE (1,IEL1) = 0.D0
00278         FLUSCE (2,IEL1) = 0.D0
00279         FLUSCE (3,IEL1) = 0.D0
00280 !
00281         FLUSCE (1,IEL2) = 0.D0
00282         FLUSCE (2,IEL2) = 0.D0
00283         FLUSCE (3,IEL2) = 0.D0
00284 !
00285 !
00286 !
00287 !   --->TEST ON THE SIGN OF LAMBDAM
00288 !       ----------------------------
00289 !
00290         IF ( RLAMBM . LT . 0.D0 ) THEN
00291 !       - - - - - - - - - - - - - -
00292 !
00293 !----------CALCUL DES TERMES SOURCES --------------------------
00294 !
00295           PSA = TS11(1)*GE(1)+TS11(2)*GE(2)+TS11(3)*GE(3)
00296 !
00297           FLUSCE(1,IEL1) = PSA*T11(1)
00298           FLUSCE(2,IEL1) = PSA*T11(2)
00299           FLUSCE(3,IEL1) = PSA*T11(3)
00300 !
00301 !
00302           PSA1= TS12(1)*GE(1)+TS12(2)*GE(2)+TS12(3)*GE(3)
00303           PSA2= TS22(1)*GE(1)+TS22(2)*GE(2)+TS22(3)*GE(3)
00304 !
00305 !
00306           FLUSCE(1,IEL2) = (PSA1*T12(1)+PSA2*T22(1))
00307           FLUSCE(2,IEL2) = (PSA1*T12(2)+PSA2*T22(2))
00308           FLUSCE(3,IEL2) = (PSA1*T12(3)+PSA2*T22(3))
00309 !
00310         ELSE
00311 !           -----
00312 !
00313 !
00314           FLUSCE(1,IEL1) = 0.D0
00315           FLUSCE(2,IEL1) = 0.D0
00316           FLUSCE(3,IEL1) = 0.D0
00317 !
00318           FLUSCE(1,IEL2) = GE(1)*CT2*2.D0
00319           FLUSCE(2,IEL2) = GE(2)*CT2*2.D0
00320           FLUSCE(3,IEL2) = GE(3)*CT2*2.D0
00321 !
00322         ENDIF
00323 !          -----
00324 !      TESTEST
00325       ELSE
00326 !      TESTEST
00327 !
00328 !   --->   SMALL CALCULATIONS
00329 !          --------------
00330 !
00331         RLAMBP = RLAMB0 + CT
00332         ALPHA = UI * XN + VI * YN
00333 !
00334         CI2 = G*HI
00335         CI = SQRT (CI2)
00336         CJ2 =  G*HJ
00337         CJ = SQRT (CJ2)
00338         RLAMBI = ALPHA - CI
00339         RLAMBJ = UJ * XN + VJ * YN - CJ
00340 !
00341         IF(RLAMBI .LT. 0.D0 .AND. RLAMBJ .GT. 0.D0) THEN
00342           RLAMBP = MAX(0.D0,RLAMBP) + ABS(RLAMBI - RLAMBJ)/4.D0
00343         ENDIF
00344 !
00345 !-----------COMPUTATION OF SOURCE TERMS --------------------------
00346 !
00347         FLUSCE(1,IEL1) = GE(1)*CT2*2.D0
00348         FLUSCE(2,IEL1) = GE(2)*CT2*2.D0
00349         FLUSCE(3,IEL1) = GE(3)*CT2*2.D0
00350 !
00351         FLUSCE(1,IEL2) = 0.D0
00352         FLUSCE(2,IEL2) = 0.D0
00353         FLUSCE(3,IEL2) = 0.D0
00354 !
00355 !
00356 !   --->    TEST ON THE SIGN OF LAMBDAP
00357 !           ----------------------------
00358 !
00359         IF ( RLAMBP . GT . 0.D0 ) THEN
00360 !       - - - - - - - - - - - - - -
00361 !
00362 !-----------COMPUTATION OF SOURCE TERMS------------------
00363 !
00364 !
00365           FLUSCE(1,IEL1) = 0.D0
00366           FLUSCE(2,IEL1) = 0.D0
00367           FLUSCE(3,IEL1) = 0.D0
00368 !
00369           FLUSCE(1,IEL2) = 0.D0
00370           FLUSCE(2,IEL2) = 0.D0
00371           FLUSCE(3,IEL2) = 0.D0
00372           PSA1= TS11(1)*GE(1)+TS11(2)*GE(2)+TS11(3)*GE(3)
00373           PSA2= TS21(1)*GE(1)+TS21(2)*GE(2)+TS21(3)*GE(3)
00374 !
00375 !
00376           FLUSCE(1,IEL1) = (PSA1*T11(1)+PSA2*T21(1))
00377           FLUSCE(2,IEL1) = (PSA1*T11(2)+PSA2*T21(2))
00378           FLUSCE(3,IEL1) = (PSA1*T11(3)+PSA2*T21(3))
00379 !
00380           PSA = TS12(1)*GE(1)+TS12(2)*GE(2)+TS12(3)*GE(3)
00381 !
00382           FLUSCE(1,IEL2) = PSA*T12(1)
00383           FLUSCE(2,IEL2) = PSA*T12(2)
00384           FLUSCE(3,IEL2) = PSA*T12(3)
00385 !
00386         ENDIF
00387 !       -----
00388 !
00389 !       TESTEST
00390       ENDIF
00391 !       TESTEST
00392       FLUSCE(1,IEL1)=FLUSCE(1,IEL1)*A1/CT2
00393       FLUSCE(2,IEL1)=FLUSCE(2,IEL1)*A1/CT2
00394       FLUSCE(3,IEL1)=FLUSCE(3,IEL1)*A1/CT2
00395       FLUSCE(1,IEL2)=FLUSCE(1,IEL2)*A2/CT2
00396       FLUSCE(2,IEL2)=FLUSCE(2,IEL2)*A2/CT2
00397       FLUSCE(3,IEL2)=FLUSCE(3,IEL2)*A2/CT2
00398 !
00399 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00400 !
00401 !
00402       RETURN
00403       END

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