flux_waf.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\flux_waf.f
00002 !
00056                         SUBROUTINE FLUX_WAF
00057 !                       *******************
00058 
00059      &(XI,H1,H2,U1,U2,V1,V2,PSI1,PSI2,
00060      & HL_UP,HR_UP,VL_UP,VR_UP,PSIL_UP,PSIR_UP,
00061      & XNN,YNN,DT,DX,WAFFLX)
00062 !
00063 !***********************************************************************
00064 ! TELEMAC 2D VERSION 6.3                                         R. ATA
00065 !
00066 !***********************************************************************
00067 !
00068 !     FUNCTION  : SUBROUTINE COMPUTES WAF FLUX: THREE HYDRODYNAMICAL
00069 !                 COMPENENTS + TRACER TRANSPORT
00070 !      SEE TORO: SHOCK CAPTURING METHODS FOR FREE
00071 !            SURFACE FLOWS (WILEY 2005)
00072 !
00073 !
00074 !
00075 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00076 ! |  DT            | -->|  TIME STEP                                   |
00077 ! |  DX            | -->|  CHAACTERISTIC LENGTH (FOR COURANT NUMBER)   |
00078 ! |  H1,H2         | -->|  LEFT AND RIGHT WATER DEPTHS                 |
00079 ! |  HL_UP,HR_UP   | -->|  LEFT AND RIGHT NEGHBORS WATER DEPTHS        |
00080 ! |  PSI1,PSI2     | -->|  LEFT AND RIGHT TRACER DENSITIES             |
00081 ! |  PSIL_UP,PSIR_UP| -->|  LEFT AND RIGHT NEIGHBORS TRACER DENSITIES  |
00082 ! |  ROT           | -->|  EXECUTE FINAL ROTATION OR NO                |
00083 ! |  U1,U2         | -->|  LEFT AND RIGHT VELOCITY X-COMPONENTS        |
00084 ! |  V1,V2         | -->|  LEFT AND RIGHT VELOCITY Y-COMPONENTS        |
00085 ! |  VL_UP,VR_UP   | -->|  LEFT AND RIGHT NEIGHBOR VELOCITY COMPONENTS |
00086 ! |  XI            | -->|  LOCATION OF THE INTERFACE IN RIEMANN DIAGRAM|
00087 ! |  XNN,YNN       | -->|  X AND Y COMPONENT OF THE OUTWARD NORMAL     |
00088 ! |  WAFFLX        | <--|  FLUX COMPONENTS AT THE INTERFACE            |
00089 ! ______________________________________________________________________
00090 !
00091 !  MODE: -->(UNCHANGEABLE INPUT),<--(OUTPUT),<-->(CHANGEABLE INPUT)
00092 !-----------------------------------------------------------------------
00093 !  CALLING SUBROUTINE FLUX_WAF OR FLUX_HLLC OR FLUXZZ
00094 !
00095 !***********************************************************************
00096 !
00097       USE BIEF
00098       USE INTERFACE_TELEMAC2D,ONLY: LIMITER,FLUX_HLLC
00099 !
00100       IMPLICIT NONE
00101       INTEGER LNG,LU
00102       COMMON/INFO/LNG,LU
00103 !
00104 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00105 !
00106       DOUBLE PRECISION, INTENT(IN)    :: H1,H2,U1,U2,PSI1,PSI2
00107       DOUBLE PRECISION, INTENT(IN)    :: XI,V1,V2,XNN,YNN,DT
00108       DOUBLE PRECISION, INTENT(IN)    :: HL_UP,HR_UP,VL_UP,VR_UP
00109       DOUBLE PRECISION, INTENT(IN)    :: PSIL_UP,PSIR_UP
00110       DOUBLE PRECISION, INTENT(IN)    :: DX
00111       DOUBLE PRECISION, INTENT(INOUT) :: WAFFLX(4)
00112 !
00113 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00114 !
00115       INTEGER                         :: I,ILIM
00116       LOGICAL                         :: ROT,TVD
00117 !
00118       DOUBLE PRECISION, PARAMETER     :: G = 9.81D0
00119       DOUBLE PRECISION                :: HL,UL,VL,PSI_L
00120       DOUBLE PRECISION                :: HR,UR,VR,PSI_R
00121       DOUBLE PRECISION                :: AL,AR,HSTAR,USTAR
00122       DOUBLE PRECISION                :: PQL,PQR,SL,SR
00123       DOUBLE PRECISION                :: FL(4),FR(4)
00124 !
00125       DOUBLE PRECISION                :: GSUR2,EPS,DTDX
00126       DOUBLE PRECISION                :: CL,CR,CSTAR,WL,WR
00127       DOUBLE PRECISION                :: WLR,WLSTAR,WRSTAR
00128       DOUBLE PRECISION                :: FLU2X,FLU2Y
00129       DOUBLE PRECISION                :: U0,SSTAR
00130       DOUBLE PRECISION                :: FLX(4),HLLCFLX(4)
00131 !
00132       DOUBLE PRECISION                :: LIM_RL,LIM_RR,LIM_RSTAR
00133       DOUBLE PRECISION                :: RL,RR,RSTAR,DELTA
00134 !
00135       INTRINSIC SIGN
00136 !
00137 !-----------------------------------------------------------------------
00138 !
00139       EPS   = 1.E-6
00140       GSUR2 = G/2.0D0
00141       ROT   = .FALSE.
00142       TVD   = .TRUE.
00143       ILIM  = 4
00144       PQL   = 0.D0
00145       PQR   = 0.D0
00146       USTAR = 0.D0
00147       HSTAR = 0.D0
00148       AL    = 0.D0
00149       AR    = 0.D0
00150 !
00151 !***********************************************************************
00152 !     INITIALIZATION OF FLX, HLLCFLX AND WAFFLX
00153       DO I=1,4
00154         FLX(I)     = 0.D0
00155         HLLCFLX(I) = 0.D0
00156         WAFFLX(I)  = 0.D0
00157       ENDDO
00158 !
00159 !-----------------------------------------------------------------------
00160 !     DEPTHS, VELOCITIES, TRACERS
00161       HL    = H1
00162       UL    = U1
00163       VL    = V1
00164       PSI_L = PSI1
00165 !
00166       HR    = H2
00167       UR    = U2
00168       VR    = V2
00169       PSI_R = PSI2
00170 !
00171 ! LET'S START BY COMPUT GNIHLLC FLUX (WITHOUT INVERSE ROTATION IN THE END)
00172 !
00173       CALL FLUX_HLLC(XI,HL,HR,UL,UR,VL,VR,PSI_L,PSI_R,
00174      &               XNN,YNN,ROT,HLLCFLX)
00175 !
00176 ! ROTATION
00177 !
00178       U0  = UL
00179       UL  = XNN*U0+YNN*VL
00180       VL  =-YNN*U0+XNN*VL
00181 !
00182       U0  = UR
00183       UR  = XNN*U0+YNN*VR
00184       VR  =-YNN*U0+XNN*VR
00185 !
00186 ! CASE WITH DRY LEFT AND RIGHT
00187       IF(HL.LT.EPS.AND.HR.LT.EPS) GOTO 20
00188 !
00189 ! CELERITIES
00190 !
00191       AL = SQRT(G*HL)
00192       AR = SQRT(G*HR)
00193 !
00194 ! STAR VARIABLES
00195 !
00196       HSTAR = 0.5D0*(HL+HR)-0.25D0*(UR-UL)*(HL+HR)/(AL+AR)
00197 !RA BUG FIXED WHEN COMPUTING U STAR, THANKS TO L.STADLER (BAW)
00198 !     USTAR = 0.5D0*(UL+UR)-0.25D0*(HR-HL)*(AL+AR)/(HL+HR)
00199       USTAR = 0.5D0*(UL+UR)-       (HR-HL)*(AL+AR)/(HL+HR)
00200 
00201 ! COMPUTE PQL AND PQR:
00202 ! IT DEPENDS IF WE ARE IN PRESENCE OF SHOCK OR RAREFACTION WAVE
00203       IF(HSTAR.LT.HL)THEN
00204 !       RAREFACTION
00205         PQL = 1.0D0
00206       ELSE
00207 !       SHOCK
00208         IF(HL.GT.EPS)THEN
00209           PQL = SQRT(0.5D0*(HSTAR+HL)*HSTAR/HL**2)
00210         ELSE
00211           PQL = 0.0D0
00212         ENDIF
00213       ENDIF
00214       IF(HSTAR.LT.HR)THEN
00215 !       RAREFACTION
00216         PQR = 1.0D0
00217       ELSE
00218 !       SHOCK
00219         IF(HR.GT.EPS)THEN
00220           PQR = SQRT(0.5D0*(HSTAR+HR)*HSTAR/HR**2)
00221         ELSE
00222           PQR = 0.0D0
00223         ENDIF
00224       ENDIF
00225 !
00226 20    CONTINUE
00227 !
00228 ! FL AND FR
00229 !
00230       FL(1)   = HL*UL
00231       FL(2)   = HL*UL**2 +GSUR2*HL**2
00232       FL(3)   = HL*UL*VL
00233       FL(4)   = HL*UL*PSI_L
00234 !
00235       FR(1)   = HR*UR
00236       FR(2)   = HR*UR**2 +GSUR2*HR**2
00237       FR(3)   = HR*UR*VR
00238       FR(4)   = HR*UR*PSI_R
00239 !
00240 !     SL, SR AND SSTAR  (WE CONSIDER DRY CASES)
00241 !
00242       IF(HL.GT.EPS) THEN
00243         SL    = UL-AL*PQL
00244       ELSE
00245         SL    = UR - 2.D0*AR
00246         SR    = UR + AR
00247         SSTAR = SL
00248       ENDIF
00249 !
00250       IF(HR.GT.EPS)THEN
00251         SR    = UR + AR*PQR
00252       ELSE
00253         SL    = UL - AL
00254         SR    = UL + 2.D0*AL
00255         SSTAR = SR
00256         GOTO 35
00257       ENDIF
00258       SSTAR   = USTAR
00259 
00260 35    CONTINUE
00261 !
00262 ! WEIGHTING COEFFICIENTS WL,WLR, WR WLSTAR AND WRSTAR
00263 !
00264 !     COURANT NUMBERS FOR ALL WAVES
00265       DTDX  = DT/DX
00266       CL    = SL*DTDX
00267       CR    = SR*DTDX
00268       CSTAR = SSTAR*DTDX
00269 !
00270 !===================================================
00271 !   NON TVD WAF SCHEME
00272 !===================================================
00273 !
00274       IF(.NOT.TVD) THEN
00275 !
00276 !     COEFFICIENTS
00277       WL     = 0.5D0*(1.D0 + CL)
00278       WR     = 0.5D0*(1.D0 - CR)
00279       WLR    = 0.5D0*(CR - CL)
00280       WLSTAR = 0.5D0*(1.D0 + CSTAR)
00281       WRSTAR = 0.5D0*(1.D0 - CSTAR)
00282 !
00283 !     FINAL FLUX (BEFORE ROTATION)
00284 !
00285       FLX(1) = WL*FL(1) + WLR*HLLCFLX(1) + WR*FR(1)
00286       FLX(2) = WL*FL(2) + WLR*HLLCFLX(2) + WR*FR(2)
00287       FLX(3) = (WLSTAR*VL + WRSTAR*VR)*FLX(1)
00288       FLX(4) = (WLSTAR*PSI_L + WRSTAR*PSI_R)*FLX(1)
00289 !
00290 !===================================================
00291 !    TVD WAF SCHEME
00292 !===================================================
00293 !
00294       ELSE
00295 !
00296 !     LIMITERS
00297 !     PREPARE rK BEFORE CALLING LIMITER
00298 !     COMPUTE ALL rK (SEE LOUKILI ET AL. PAGE 4)
00299 !       RL
00300         IF(SL.GT.0.D0)THEN
00301           DELTA = HL-HL_UP
00302         ELSE
00303           DELTA = HR_UP-HR
00304         ENDIF
00305         RL = DELTA/(HR-HL + EPS)
00306 !       RR
00307         IF(SR.GT.0.0D0)THEN
00308           DELTA = HL-HL_UP
00309         ELSE
00310           DELTA = HR_UP-HR
00311         ENDIF
00312         RR = DELTA/(HR-HL + EPS)
00313 !       r*
00314         IF(SSTAR.GT.0.D0)THEN
00315           DELTA = VL-VL_UP
00316         ELSE
00317           DELTA = VR_UP-VR
00318         ENDIF
00319         RSTAR = DELTA/(VR-VL+EPS)
00320 !
00321         LIM_RL    = LIMITER(ILIM,RL,CL)
00322         LIM_RR    = LIMITER(ILIM,RR,CR)
00323         LIM_RSTAR = LIMITER(ILIM,RSTAR,CSTAR)
00324 !
00325 !   TVD COEFFICIENTS
00326 !
00327       WL     = 0.5D0*(1.D0 + SIGN(1.D0,CL)*LIM_RL) !DSIGN(A,B)=|A|*SIGN(B)
00328       WR     = 0.5D0*(1.D0 - SIGN(1.D0,CR)*LIM_RR)
00329       WLR    = 0.5D0*(SIGN(1.D0,CR)*LIM_RR - SIGN(1.D0,CL)*LIM_RL)
00330       WLSTAR = 0.5D0*(1.D0 + SIGN(1.D0,CSTAR)*LIM_RSTAR)
00331       WRSTAR = 0.5D0*(1.D0 - SIGN(1.D0,CSTAR)*LIM_RSTAR)
00332 !
00333 ! FINAL FLUX (BEFORE ROTATION)
00334 !
00335       FLX(1) = WL*FL(1) + WLR*HLLCFLX(1) + WR*FR(1)
00336       FLX(2) = WL*FL(2) + WLR*HLLCFLX(2) + WR*FR(2)
00337       FLX(3) = (WLSTAR*VL    + WRSTAR*VR   )*FLX(1)
00338       FLX(4) = (WLSTAR*PSI_L + WRSTAR*PSI_R)*FLX(1)
00339 !
00340       ENDIF
00341 !
00342 !
00343 ! INVERSE ROTATION
00344 !
00345       FLU2X  = XNN*FLX(2) - YNN*FLX(3)
00346       FLU2Y  = YNN*FLX(2) + XNN*FLX(3)
00347 !
00348 ! FINAL WAF FLUX
00349 !
00350       WAFFLX(1) = FLX(1)
00351       WAFFLX(2) = FLU2X
00352       WAFFLX(3) = FLU2Y
00353       WAFFLX(4) = FLX(4)
00354 !
00355 !
00356 !-----------------------------------------------------------------------
00357 !
00358       RETURN
00359       END

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