hyd_waf.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\hyd_waf.f
00002 !
00059                         SUBROUTINE HYD_WAF
00060 !                       *******************
00061 
00062      &(NS,NSEG,NELEM,NUBO,G,W,ZF,VNOCL,DT,DTHAUT,
00063      & X,Y,CE,ELTSEG,NEISEG)
00064 !
00065 !***********************************************************************
00066 ! TELEMAC 2D VERSION V6P3                                     15/01/2013
00067 !***********************************************************************
00068 !
00069 !
00070 !     FUNCTION  : COMPUTES ALL THE FLUXES FOR INTERNAL INTERFACES USING
00071 !                 WAF FLUX.
00072 !
00073 !
00074 !
00075 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00076 ! |  CE            |<-->|  FLUX  INCREMENTS AT INTERNAL FACES          |                          |
00077 ! |  DT            | -->|  TIME STEP                                   |
00078 ! |  ELTSEG        | -->|  SEGMENT NUMBERS PER ELEMENT                 |
00079 ! |  G             | -->|  GRAVITY                                     |
00080 ! |  NELEM         | -->|  NUMBER OF TOTAL ELEMENTS                    |
00081 ! |  NS            | -->|  NUMBER OF TOTAL MESH NODES                  |
00082 ! |  NSEG          | -->|  NUMBER OF TOTAL MESH EDGES                  |
00083 ! |  NUBO          ! -->!  GLOBAL NUMBER OF EDGE EXTREMITIES           |
00084 ! |  NEISEG        | -->| LEFT & RIGHT NEIGHBOUR SEGMENTS OF A SEGMENT |
00085 ! !  VNOCL         | -->|  OUTWARD UNIT NORMALS                        |
00086 ! !                |    |   (2 FIRST COMPONENTS) AND                   |
00087 ! !                |    |   SEGMENT LENGTH  (THIRD COMPONENT)          |
00088 ! |  W             | -->|  W(1,IS) = H,  W(2,IS)=U  ,W(3,IS)=V         |
00089 ! |  X,Y           | -->|  X AND Y COORDINATES                         |
00090 ! |  ZF            | -->|  BATHYMETRIES                                |
00091 ! !________________|____|______________________________________________|
00092 !  MODE: -->( UNCHANGEABLE INPUT ),<--(OUTPUT),<-->(CHANGEABLE INPUT)
00093 !-----------------------------------------------------------------------
00094 !     - CALLINF SUBROUTINE(S) : RESOLU
00095 !
00096 !***********************************************************************
00097 !
00098       USE INTERFACE_TELEMAC2D, EX_HYD_WAF => HYD_WAF
00099       USE BIEF_DEF,ONLY: NCSIZE
00100 !
00101       IMPLICIT NONE
00102       INTEGER LNG,LU
00103       COMMON/INFO/LNG,LU
00104 !
00105 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00106 !
00107       INTEGER, INTENT(IN) :: NS,NSEG,NELEM
00108       INTEGER, INTENT(IN) :: NUBO(2,NSEG),NEISEG(2,NSEG)
00109       INTEGER, INTENT(IN)             :: ELTSEG(NELEM,3)
00110       DOUBLE PRECISION, INTENT(IN)    :: ZF(NS),VNOCL(3,NSEG)
00111       DOUBLE PRECISION, INTENT(IN)    :: X(NS),Y(NS)
00112       DOUBLE PRECISION, INTENT(IN)    :: G,W(3,NS)
00113       DOUBLE PRECISION, INTENT(IN)    :: DTHAUT(*)
00114       DOUBLE PRECISION, INTENT(INOUT) :: CE(NS,3),DT
00115 !
00116 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00117 !
00118       INTEGER    :: NSG,NUBO1,NUBO2,IVAR,I,IDRY,IEL
00119       INTEGER    :: NUBOL,NUBOR,SEG1,SEG2,IER
00120 !
00121       DOUBLE PRECISION :: H1,H2,EPS,FLX(4)
00122       DOUBLE PRECISION :: ZF1,ZF2,XNN,YNN,RNN
00123       DOUBLE PRECISION :: V21,V22,V31,V32,DX
00124       DOUBLE PRECISION :: HIJ,HJI,DZIJ,DZJI
00125       DOUBLE PRECISION :: HGZI,HGZJ,HDXZ1,HDYZ1,HDXZ2,HDYZ2
00126       DOUBLE PRECISION :: HL_UP,HR_UP
00127       DOUBLE PRECISION :: VL_UP,VR_UP
00128       DOUBLE PRECISION :: PROD_SCAL
00129       DOUBLE PRECISION :: PSIL_UP,PSIR_UP
00130       LOGICAL, ALLOCATABLE  :: YESNO(:)
00131 !  XI = X/T, AT THE (X,T) DIAGRAMM, WE COMPUETE THE FLUX AT XI = 0
00132       DOUBLE PRECISION,PARAMETER :: XI=0.0D0
00133 ! PSI1 AND PSI2 ARE TRACER
00134 ! NOT ADAPTED YET
00135       DOUBLE PRECISION           :: PSI1,PSI2
00136 !
00137 !**************************************************************
00138       EPS=1.E-6
00139 ! NO TRACER UNTIL NOW ...
00140       PSI1 = 0.0D0
00141       PSI2 = 0.0D0
00142       ALLOCATE(YESNO(NSEG),STAT=IER)
00143       IF(IER.NE.0)THEN
00144         IF(LNG.EQ.1)WRITE(LU,*)'FLUX_TCH: ERREUR D''ALLOCATION'
00145         IF(LNG.EQ.2)WRITE(LU,*)'FLUX_TCH: ALLOCATION ERROR '
00146         CALL PLANTE(1)
00147         STOP
00148       ENDIF
00149 !**************************************************************
00150 ! INITIALIZATION OF CE
00151       DO I=1,3
00152         DO IVAR=1,NS
00153           CE(IVAR,I) = 0.D0
00154         ENDDO
00155       ENDDO
00156 ! INITIALIZATION OF YESNO
00157       DO I=1,NSEG
00158         YESNO(I)=.FALSE.
00159       ENDDO
00160 !
00161 !-----------------------------------------------------------------------
00162 !     LOOP OVER GLOBAL LIST OF EDGES
00163 !    *******************************
00164 !
00165       DO IEL=1,NELEM
00166         DO I = 1,3
00167           IF(.NOT.YESNO(ELTSEG(IEL,I)))THEN
00168             NSG = ELTSEG(IEL,I)
00169 !           INDICATOR FOR DRY CELLS
00170             IDRY=0
00171 !           INITIALIZATION
00172             FLX(1) = 0.0D0
00173             FLX(2) = 0.0D0
00174             FLX(3) = 0.0D0
00175             FLX(4) = 0.0D0
00176 !
00177 !           RECUPERATE NODES OF THE EDGE WITH THE GOOD ORIENTATION
00178 !           WITH RESPECT TO THE NORMAL
00179             NUBO1 = NUBO(1,NSG)
00180             NUBO2 = NUBO(2,NSG)
00181             PROD_SCAL= ((X(NUBO2)-X(NUBO1))*VNOCL(1,NSG)+
00182      &                  (Y(NUBO2)-Y(NUBO1))*VNOCL(2,NSG))
00183             IF(PROD_SCAL.LT.0.D0)THEN
00184               NUBO1 = NUBO(2,NSG)
00185               NUBO2 = NUBO(1,NSG)
00186             ENDIF
00187 !           THEIR BATHYMETRIES
00188             ZF1   =    ZF(NUBO1)
00189             ZF2   =    ZF(NUBO2)
00190 !           MEAN DISTANCE BETWEEN THEM (FOR CFL)
00191             DX    = 0.5D0*(DTHAUT(NUBO1)+DTHAUT(NUBO2))
00192 !           NORMAL COORDINATES NX, NY AND SEGMENT LENGTH
00193             XNN       = VNOCL(1,NSG)
00194             YNN       = VNOCL(2,NSG)
00195             RNN       = VNOCL(3,NSG)
00196 !
00197 !           WATER DEPTH
00198 !
00199             H1=W(1,NUBO1)
00200             H2=W(1,NUBO2)
00201 !******************************************************
00202 !           HYDROSTATIC RECONSTRUCTION !!!
00203 !
00204 !           BATHY AT THE INTERFACE
00205 !
00206             DZIJ = MAX(0.D0,ZF2-ZF1)
00207             HIJ  = MAX(0.D0,H1- DZIJ)
00208 !******************************************************
00209 !           HYDROSTATIC RECONSTRUCTION !!!
00210 !
00211             DZJI = MAX(0.D0,ZF1-ZF2)
00212             HJI  = MAX(0.D0,H2- DZJI)
00213 !******************************************************
00214 !
00215 !           VELOCITY COMPONENTS
00216 !
00217             IF(H1.GT.EPS)THEN
00218               V21 = W(2,NUBO1)/H1
00219               V31 = W(3,NUBO1)/H1
00220             ELSE
00221               V21=0.0D0
00222               V31=0.0D0
00223               IDRY=IDRY+1
00224             ENDIF
00225 !
00226             IF(H2.GT.EPS)THEN
00227               V22 = W(2,NUBO2)/H2
00228               V32 = W(3,NUBO2)/H2
00229             ELSE
00230               V22=0.0D0
00231               V32=0.0D0
00232               IDRY=IDRY+1
00233             ENDIF
00234 !
00235 !           SEGMENT NEIGHBORS (FOR LIMITER)
00236 !
00237             SEG1 = NEISEG(1,NSG)
00238             SEG2 = NEISEG(2,NSG)
00239 !             VERIFY THAT WE HAVE THE GOOD NEIGHBORS
00240             IF((SEG1.LE.0.OR.SEG2.LE.0).AND.NCSIZE.LE.1)THEN
00241               WRITE(LU,*)'PROBLEM TO FIND SEGMENT NEIGHBORS'
00242               WRITE(LU,*)'WE ARE IN HYD_WAF.F'
00243               WRITE(LU,*)'SEGMENT OF INTEREST  :',NSG
00244               WRITE(LU,*)'NEIGHBORS ARE  :',SEG1,SEG2
00245               CALL PLANTE(1)
00246               STOP
00247             ENDIF
00248 !
00249             NUBOL=0
00250             NUBOR=0
00251 !
00252             IF(NUBO(1,SEG1).EQ.NUBO1) THEN
00253               NUBOL = NUBO(2,SEG1)
00254             ELSE
00255               NUBOL = NUBO(1,SEG1)
00256             ENDIF
00257             IF(NUBO(1,SEG2).EQ.NUBO2) THEN
00258               NUBOR = NUBO(2,SEG2)
00259             ELSE
00260               NUBOR = NUBO(1,SEG2)
00261             ENDIF
00262 !             VERIFY THAT WE HAVE THE GOOD NEIGHBORS
00263             IF(NUBOL.LE.0.OR.NUBOR.LE.0) THEN
00264               WRITE(LU,*)'PROBLEM TO FIND NEIGHBOR'
00265               WRITE(LU,*)'WE ARE IN HYD_WAF.F'
00266               WRITE(LU,*)'NODES ARE  :',NUBO1,NUBO2
00267               WRITE(LU,*)'NEIGHBORS ARE  :', NUBOR,NUBOL
00268               CALL PLANTE(1)
00269               STOP
00270             ENDIF
00271 !
00272 !           WATER DEPTH, VELOCITY, TRACER OF NEIGHBORS
00273 !
00274             HL_UP   = W(1,NUBOL)
00275             HR_UP   = W(1,NUBOR)
00276             IF(HL_UP.GT.EPS)THEN
00277               VL_UP =  W(3,NUBOL)/HL_UP
00278             ELSE
00279               VL_UP = 0.0D0
00280             ENDIF
00281             IF(HR_UP.GT.EPS)THEN
00282               VR_UP = W(3,NUBOR)/HR_UP
00283             ELSE
00284               VR_UP = 0.0D0
00285             ENDIF
00286             PSIL_UP = PSI1
00287             PSIR_UP = PSI2
00288 !
00289 !           LOCAL FLUX COMPUTATION
00290 !
00291 !           AT LEAST ONE WET CELL
00292             IF(IDRY.LT.2)THEN
00293 !           CALL FLUX_WAF(XI,H1,H2,V21,V22,V31,V32,PSI1,PSI2,
00294               CALL FLUX_WAF(XI,HIJ,HJI,V21,V22,V31,V32,PSI1,PSI2,
00295      &                      HL_UP,HR_UP,VL_UP,VR_UP,PSIL_UP,PSIR_UP,
00296      &                      XNN,YNN,DT,DX,FLX)
00297 !
00298 !             GEOMETRIC SOURCE TERMS
00299 !
00300               HGZI =0.5D0*RNN*(HIJ+H1)*(HIJ-H1)
00301               HGZJ =0.5D0*RNN*(HJI+H2)*(HJI-H2)
00302 !
00303               HDXZ1  = G*XNN*HGZI
00304               HDYZ1  = G*YNN*HGZI
00305 !
00306               HDXZ2  = G*XNN*HGZJ
00307               HDYZ2  = G*YNN*HGZJ
00308 !
00309 !             FLUX INCREMENT
00310 !
00311               CE(NUBO1,1) = CE(NUBO1,1) - RNN*FLX(1)
00312               CE(NUBO1,2) = CE(NUBO1,2) - RNN*FLX(2) + HDXZ1
00313               CE(NUBO1,3) = CE(NUBO1,3) - RNN*FLX(3) + HDYZ1
00314 !
00315               CE(NUBO2,1) = CE(NUBO2,1) + RNN*FLX(1)
00316               CE(NUBO2,2) = CE(NUBO2,2) + RNN*FLX(2) - HDXZ2
00317               CE(NUBO2,3) = CE(NUBO2,3) + RNN*FLX(3) - HDYZ2
00318             ENDIF
00319 !
00320             YESNO(NSG)=.TRUE.
00321           ENDIF
00322         ENDDO
00323       ENDDO
00324 !
00325       DEALLOCATE(YESNO)
00326 !-----------------------------------------------------------------------
00327       RETURN
00328       END

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