fluroe.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\fluroe.f
00002 !
00069                      SUBROUTINE FLUROE
00070 !                    *****************
00071 !
00072      &(W,FLUSCE,NUBO,VNOIN,WINF,FLUX,FLUSORT,FLUENT,
00073      & NELEM,NSEG,NPTFR,NPOIN,X,Y,AIRS,ZF,EPS,DDMIN,G,
00074      & XNEBOR,YNEBOR,LIMPRO,NBOR,KDIR,KNEU,KDDL,FLBOR,
00075      & ELTSEG,IFABOR,MESH)
00076 !
00077 !***********************************************************************
00078 ! TELEMAC2D   V6P3                                          21/06/2013
00079 !***********************************************************************
00080 !
00081 !
00082 !
00083 !
00084 !
00085 !
00086 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00087 !| AIRS           |-->| CELL AREA
00088 !| DDMIN          |<--| MINIMUM DISTANCE
00089 !| EPS            |-->| TOLERANCE
00090 !| FLUENT         |<--| MASS FLUX MASSE INLET
00091 !| FLUSCE         |-->| SOURCE FLUXES
00092 !| FLUSORT        |<--| MASS FLUX MASSE OUTLET
00093 !| FLUX           |<--| ROE FLUX
00094 !| G              |-->| GRAVITY
00095 !| KDDL           |-->| CONVENTION FOR FREE POINTS (BC)
00096 !| KDIR           |-->| CONVENTION FOR DIRICHLET POINTS
00097 !| KFROT          |-->| BED FRICTION LAW
00098 !| KNEU           |-->| CONVENTION NEUMANN POINTS
00099 !| LIMPRO         |-->| TYPES OF BOUNDARY CONDITION!
00100 !| NBOR           |-->| NUMER OF BORD NODES
00101 !| NELEM          |-->| NUMBER OF ELEMENTS
00102 !| NPTFR          |-->| TOTAL NUMBER OF BOUNDARY NODES
00103 !| NREJET         |-->| NUMBER OF SOURCE/SINK
00104 !| NSEG           |-->| NUMBER OF EDGES
00105 !| NUBO           |-->| GLOBAL INDICES OF EDGE EXTREMITIES
00106 !| VNOIN          |-->| NORMAL TO THE INTERFACE
00107 !|                |   | (2 FIRS COMPOSANTES) AND
00108 !|                |   | SEGMENT LENGTH (3RD COMPONENT)
00109 !| W              |<->| WORKING TABLE
00110 !| WINF           |-->| PRESCRIBED VALUES AT THE INLET AND OUTLET
00111 !| X,Y            |-->| COORDINATES FOR MESH NODES
00112 !| XNEBOR,YNEBOR  |-->| NORMAL TO BOUNDARY POINTS
00113 !| ZF             |-->| BED TOPOGRAPHY (BATHYMETRY)
00114 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00115 !
00116       USE BIEF
00117       USE INTERFACE_TELEMAC2D, EX_FLUROE => FLUROE
00118 !
00119       IMPLICIT NONE
00120       INTEGER LNG,LU
00121       COMMON/INFO/LNG,LU
00122 !
00123 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00124 !
00125       INTEGER, INTENT(IN) :: NPOIN,NELEM,NSEG,NPTFR,KDIR,KNEU,KDDL
00126       INTEGER, INTENT(IN) :: NUBO(2,*),LIMPRO(NPTFR,6),NBOR(NPTFR)
00127       INTEGER, INTENT(IN) :: ELTSEG(NELEM,3)
00128       DOUBLE PRECISION, INTENT(IN) :: X(NPOIN),Y(NPOIN),W(3,NPOIN)
00129       DOUBLE PRECISION, INTENT(IN) :: AIRS(NPOIN),ZF(NPOIN),VNOIN(3,*)
00130       DOUBLE PRECISION, INTENT(IN) :: XNEBOR(2*NPTFR),YNEBOR(2*NPTFR)
00131       DOUBLE PRECISION, INTENT(IN) :: WINF(3,NPTFR),G,EPS
00132       DOUBLE PRECISION, INTENT(INOUT) :: FLUX(NPOIN,3),DDMIN
00133       DOUBLE PRECISION, INTENT(INOUT) :: FLUSCE(3,NPOIN),FLUSORT,FLUENT
00134 ! RA
00135       TYPE(BIEF_OBJ), INTENT(INOUT)  :: FLBOR
00136       INTEGER, INTENT(IN)            :: IFABOR(NELEM,3)
00137       TYPE(BIEF_MESH),INTENT(INOUT)  :: MESH
00138 !
00139 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00140 !
00141       INTEGER NB,NSG,NUBO1,NUBO2,INDIC,K,I,IEL,IER
00142 !
00143       DOUBLE PRECISION DIJ,FLULOC(3),INFLOW,OUTFLOW,PROD_SCAL,DEMI
00144       DOUBLE PRECISION HI,UI,VI,HJ,VJ,UJ,XN,YN,RNORM,CT,PRI,USN
00145       LOGICAL, ALLOCATABLE ::  YESNO(:)
00146 !
00147       INTRINSIC SQRT
00148 !-----------------------------------------------------------------------
00149 !------
00150 ! 1. INITIALISATIONS
00151 !
00152       DEMI=0.5D0
00153       FLULOC(1) = 0.D0
00154       FLULOC(2) = 0.D0
00155       FLULOC(3) = 0.D0
00156 !
00157       ALLOCATE(YESNO(NSEG),STAT=IER)
00158       IF(IER.NE.0)THEN
00159         IF(LNG.EQ.1)WRITE(LU,*)'FLUX_TCH: ERREUR D''ALLOCATION'
00160         IF(LNG.EQ.2)WRITE(LU,*)'FLUX_TCH: ALLOCATION ERROR '
00161         CALL PLANTE(1)
00162         STOP
00163       ENDIF
00164 !
00165       DO  NB =  1 , NPOIN
00166         FLUX(NB,1) = 0.D0
00167         FLUX(NB,2) = 0.D0
00168         FLUX(NB,3) = 0.D0
00169         FLUSCE(1,NB) = 0.D0
00170         FLUSCE(2,NB) = 0.D0
00171         FLUSCE(3,NB) = 0.D0
00172       ENDDO
00173 !    INITIALIZATION OF YESNO
00174       DO I=1,NSEG
00175         YESNO(I)=.FALSE.
00176       ENDDO
00177 !
00178 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00179 !
00180 !------
00181 ! 2. COMPUTE INTERNAL FLUXES
00182 !------
00183 !
00184 !  LOOP OVER INTERNAL SEGMENTS
00185 !       ----------------------------------
00186 !
00187       DDMIN = 10000.D0
00188       DO IEL=1, NELEM
00189         DO I = 1,3
00190           IF(.NOT.YESNO(ELTSEG(IEL,I)))THEN
00191             NSG = ELTSEG(IEL,I)
00192 !
00193 !           RECUPERATE NODES OF THE EDGE WITH THE GOOD ORIENTATION
00194 !           WITH RESPECT TO THE NORMAL
00195             NUBO1 = NUBO(1,NSG)
00196             NUBO2 = NUBO(2,NSG)
00197             PROD_SCAL= ((X(NUBO2)-X(NUBO1))*VNOIN(1,NSG)+
00198      &                  (Y(NUBO2)-Y(NUBO1))*VNOIN(2,NSG))
00199             IF(PROD_SCAL.LT.0.D0)THEN
00200               NUBO1 = NUBO(2,NSG)
00201               NUBO2 = NUBO(1,NSG)
00202             ENDIF
00203             INDIC =0
00204 !
00205 !   --->    INTERMEDIATE COMPUTATIONS
00206 !
00207             HI = W(1,NUBO1)
00208             IF(HI.GT.EPS) THEN
00209               UI = W(2,NUBO1) / HI
00210               VI = W(3,NUBO1) / HI
00211             ELSE
00212               INDIC = INDIC + 1
00213               UI = 0.D0
00214               VI = 0.D0
00215               HI = 0.D0
00216             ENDIF
00217 !
00218             HJ = W(1,NUBO2)
00219             IF(HJ.GT.EPS) THEN
00220               UJ = W(2,NUBO2) / HJ
00221               VJ = W(3,NUBO2) / HJ
00222             ELSE
00223               INDIC = INDIC + 1
00224               UJ = 0.D0
00225               VJ = 0.D0
00226               HJ = 0.D0
00227             ENDIF
00228 !
00229             DIJ = SQRT((X(NUBO1)-Y(NUBO2))**2+(Y(NUBO1)-Y(NUBO2))**2)
00230             IF(DIJ.LE.DDMIN) THEN
00231               DDMIN=DIJ
00232             ENDIF
00233             RNORM = VNOIN(3,NSG)
00234 !
00235             IF(INDIC.LT.2) THEN
00236               XN = VNOIN (1,NSG)
00237               YN = VNOIN (2,NSG)
00238 !
00239               CALL FLUXE(HJ,UJ,VJ,HI,UI,VI,XN,YN,RNORM,G,FLULOC)
00240               CALL FLUSRC(NUBO1,NUBO2,NSG,VNOIN,W,FLUSCE,X,Y,AIRS,NPOIN,
00241      &                 NSEG,ZF,EPS,G)
00242             ELSE
00243               FLULOC(1)= 0.D0
00244               FLULOC(2)= 0.D0
00245               FLULOC(3)= 0.D0
00246               FLUSCE(1,NUBO1)=0.D0
00247               FLUSCE(2,NUBO1)=0.D0
00248               FLUSCE(3,NUBO1)=0.D0
00249               FLUSCE(1,NUBO2)=0.D0
00250               FLUSCE(2,NUBO2)=0.D0
00251               FLUSCE(3,NUBO2)=0.D0
00252             ENDIF
00253 !
00254 !FOR PARALLELISM
00255 !
00256             IF(NCSIZE.GT.1)THEN
00257               IF(IFABOR(IEL,I).EQ.-2)THEN !THIS IS INTERFACE EDGE
00258                 ! DEMI=DEMI*SIGN(1.0D0,PROD_SCAL)
00259                 FLULOC(1)= DEMI*FLULOC(1)
00260                 FLULOC(2)= DEMI*FLULOC(2)
00261                 FLULOC(3)= DEMI*FLULOC(3)
00262 !
00263                 FLUSCE(1,NUBO1)= DEMI*FLUSCE(1,NUBO1)
00264                 FLUSCE(2,NUBO1)= DEMI*FLUSCE(2,NUBO1)
00265                 FLUSCE(3,NUBO1)= DEMI*FLUSCE(3,NUBO1)
00266                 FLUSCE(1,NUBO2)= DEMI*FLUSCE(1,NUBO2)
00267                 FLUSCE(2,NUBO2)= DEMI*FLUSCE(2,NUBO2)
00268                 FLUSCE(3,NUBO2)= DEMI*FLUSCE(3,NUBO2)
00269               ENDIF
00270             ENDIF
00271 !
00272             FLUX(NUBO1,1)=FLUX(NUBO1,1)+FLULOC(1)*RNORM+FLUSCE(1,NUBO1)
00273             FLUX(NUBO1,2)=FLUX(NUBO1,2)+FLULOC(2)*RNORM+FLUSCE(2,NUBO1)
00274             FLUX(NUBO1,3)=FLUX(NUBO1,3)+FLULOC(3)*RNORM+FLUSCE(3,NUBO1)
00275 !
00276             FLUX(NUBO2,1)=FLUX(NUBO2,1)-FLULOC(1)*RNORM+FLUSCE(1,NUBO2)
00277             FLUX(NUBO2,2)=FLUX(NUBO2,2)-FLULOC(2)*RNORM+FLUSCE(2,NUBO2)
00278             FLUX(NUBO2,3)=FLUX(NUBO2,3)-FLULOC(3)*RNORM+FLUSCE(3,NUBO2)
00279 !
00280             YESNO(NSG)=.TRUE.
00281           ENDIF
00282 
00283         ENDDO
00284       ENDDO
00285 
00286 !     FOR PARALLESM
00287       IF(NCSIZE.GT.1)THEN
00288         CALL PARCOM2(FLUX(:,1),FLUX(:,2),FLUX(:,3),NPOIN,1,2,3,MESH)
00289       ENDIF
00290 !
00291 ! 3. COMPUTE FLUXES AT THE BOUNDARIES
00292 !     START AT ZERO
00293       FLUSORT = 0.D0
00294       FLUENT  = 0.D0
00295 !
00296 ! ====>   INLET FLUX
00297 !
00298       DO K = 1 , NPTFR
00299 !
00300         NB = NBOR(K)
00301         INDIC = 0
00302 !
00303         IF(LIMPRO(K,1).EQ.KDIR.OR.LIMPRO(K,1).EQ.KDDL) THEN
00304 !
00305 !
00306 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00307 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00308 !
00309 !   --->    INTERMEDIATE COMPUTATIONS
00310 !           -------------------------------
00311 !
00312 !         IF H IMPOSED
00313           IF(LIMPRO(K,1).EQ.KDIR) THEN
00314             HJ = WINF(1,K)
00315             IF(HJ.GT.EPS) THEN
00316               UJ = WINF(2,K) / HJ
00317               VJ = WINF(3,K) / HJ
00318             ELSE
00319               HJ = 0.D0
00320               UJ = 0.D0
00321               VJ = 0.D0
00322               INDIC = INDIC+1
00323             ENDIF
00324             HI = W(1,NB)
00325             IF(HI.GT.EPS) THEN
00326               UI = W(2,NB) / HI
00327               VI = W(3,NB) / HI
00328             ELSE
00329               HI = 0.D0
00330               UI = 0.D0
00331               VI = 0.D0
00332               INDIC = INDIC+1
00333             ENDIF
00334             XN = XNEBOR(K)
00335             YN = YNEBOR(K)
00336             RNORM = SQRT( XNEBOR(K+NPTFR)**2 + YNEBOR(K+NPTFR)**2 )
00337 !
00338             IF(INDIC.LT.2) THEN
00339               CALL FLUXE(HJ,UJ,VJ,HI,UI,VI,XN,YN,RNORM,G,FLULOC)
00340               OUTFLOW = FLULOC(1)*RNORM
00341               FLUSORT = FLUSORT + OUTFLOW
00342               FLBOR%R(K)=OUTFLOW
00343             ELSE
00344               FLULOC(1)= 0.D0
00345               FLULOC(2)= 0.D0
00346               FLULOC(3)= 0.D0
00347               FLBOR%R(K)=0.0D0
00348             ENDIF
00349 !
00350 ! RA: LIMPRO .NE. KDIR
00351           ELSE
00352             HI = WINF(1,K)
00353             IF(HI.GT.EPS) THEN
00354              UI = WINF(2,K) / HI
00355              VI = WINF(3,K) / HI
00356             ELSE
00357              HI = 0.D0
00358              UI = 0.D0
00359              VI = 0.D0
00360              INDIC = INDIC+1
00361             ENDIF
00362             HJ = W(1,NB)
00363             IF(HJ.GT.EPS) THEN
00364              UJ = W(2,NB) / HJ
00365              VJ = W(3,NB) / HJ
00366             ELSE
00367              HJ = 0.D0
00368              UJ = 0.D0
00369              VJ = 0.D0
00370              INDIC = INDIC+1
00371             ENDIF
00372             XN = XNEBOR(K)
00373             YN = YNEBOR(K)
00374             RNORM = SQRT( XNEBOR(K+NPTFR)**2 + YNEBOR(K+NPTFR)**2 )
00375 !
00376             IF(INDIC.LT.2) THEN
00377              CALL FLUXE(HI,UI,VI,HJ,UJ,VJ,XN,YN,RNORM,G,FLULOC)
00378 ! RA
00379              INFLOW = FLULOC(1)*RNORM
00380              FLUENT = FLUENT + INFLOW
00381              FLBOR%R(K)=INFLOW
00382 !
00383             ELSE
00384              FLULOC(1)= 0.D0
00385              FLULOC(2)= 0.D0
00386              FLULOC(3)= 0.D0
00387              FLBOR%R(K)=0.0D0
00388             ENDIF
00389           ENDIF
00390 !
00391           FLUX(NB,1)=FLUX(NB,1)+FLULOC(1)*RNORM
00392           FLUX(NB,2)=FLUX(NB,2)+FLULOC(2)*RNORM
00393           FLUX(NB,3)=FLUX(NB,3)+FLULOC(3)*RNORM
00394 !
00395 !------
00396 ! 5. FLUX AT THE BOUNDARY
00397 !------
00398 !
00399         ELSEIF(LIMPRO(K,1).EQ.KNEU) THEN
00400 !
00401 !
00402           IF(W(1,NB).GT.0.) THEN
00403             PRI =G*( W(1,NB)*W(1,NB))/2.D0
00404             USN = (W(2,NB)*XNEBOR(K)+W(3,NB)*YNEBOR(K))/W(1,NB)
00405             CT = SQRT(G*W(1,NB))
00406 !
00407             IF((USN+2.D0*CT).LE.0.D0) THEN
00408               PRI =0.D0
00409             ELSEIF(((USN+2.D0*CT).GE.0.D0) .AND. (USN.LE.0.D0)) THEN
00410               PRI = PRI * (1.D0 +  USN / 2.D0 / CT)
00411      &                  * (1.D0 +  USN / 2.D0 / CT)
00412      &                  * (1.D0 +  USN / 2.D0 / CT)
00413      &                  * (1.D0 +  USN / 2.D0 / CT)
00414 !
00415             ELSEIF(USN.GE.0.D0) THEN
00416               PRI = PRI*(1.D0 + 2.D0 * USN / CT)
00417             ENDIF
00418 !
00419             FLUX(NB,2) = FLUX(NB,2) + XNEBOR(K+NPTFR) * PRI
00420             FLUX(NB,3) = FLUX(NB,3) + YNEBOR(K+NPTFR) * PRI
00421           ENDIF
00422 !
00423         ENDIF
00424       ENDDO ! K
00425 !
00426       DEALLOCATE(YESNO)
00427 !-----------------------------------------------------------------------
00428 !
00429       RETURN
00430       END

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