flusec_telemac2d.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\flusec_telemac2d.f
00002 !
00076                      SUBROUTINE FLUSEC_TELEMAC2D
00077 !                    ***************************
00078 !
00079      &(U,V,H,IKLE,XEL,YEL,NELMAX,NELEM,X,Y,DT,NCP,CTRLSC,INFO,TPS,
00080      & MSKSEC,BM1,BM2,T1,HPROP,MESH,S,CV1,IFABOR,COMFLU,CUMFLO)
00081 !
00082 !***********************************************************************
00083 ! TELEMAC2D   V7P0                                   21/08/2010
00084 !***********************************************************************
00085 !
00086 !
00087 !
00088 !
00089 !
00090 !
00091 !
00092 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00093 !| BM1            |<->| WORK MATRIX
00094 !| BM2            |<->| WORK MATRIX
00095 !| COMFLU         |-->| KEYWORD: COMPATIBLE COMPUTATION OF FLUXES
00096 !| CTRLSC         |-->| DATA ON CONTROL SECTIONS
00097 !| CUMFLO         |-->| KEYWORD: PRINTING CUMULATED FLOWRATES
00098 !| CV1            |<->| WORK ARRAY IN A BIEF_OBJ STRUCTURE
00099 !| DT             |-->| TIME STEP
00100 !| H              |-->| WATER DEPTH
00101 !| HPROP          |-->| PROPAGATION DEPTH
00102 !| IFABOR         |-->| ELEMENTS BEHIND THE EDGES OF A TRIANGLE
00103 !|                |   | IF NEGATIVE OR ZERO, THE EDGE IS A LIQUID
00104 !|                |   | BOUNDARY
00105 !| IKLE           |-->| CONNECTIVITY TABLE
00106 !| INFO           |-->| IF YES, PRINT
00107 !| MESH           |-->| MESH STRUCTURE
00108 !| MSKSEC         |<->| BLOCK OF WORK ARRAYS, PER SECTION
00109 !| NCP            |-->| TWO TIMES THE NUMBER OF CONTROL SECTIONS
00110 !| NELEM          |-->| NUMBER OF ELEMENTS
00111 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00112 !| S              |-->| VOID STRUCTURE
00113 !| T1             |<->| WOR kARRAY IN A BIEF_OBJ STRUCTURE
00114 !| TPS            |-->| TIME IN SECONDS
00115 !| U,V            |-->| VELOCITY FIELD COMPONENTS
00116 !| X              |-->| ABSCISSAE OF POINTS IN THE MESH
00117 !| XEL            |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
00118 !| Y              |-->| ORDINATES OF POINTS IN THE MESH
00119 !| YEL            |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
00120 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00121 !
00122       USE BIEF
00123       USE DECLARATIONS_TELEMAC2D, ONLY: CHAIN
00124 !
00125       IMPLICIT NONE
00126       INTEGER LNG,LU
00127       COMMON/INFO/LNG,LU
00128 !
00129 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00130 !
00131       INTEGER, INTENT(IN) :: NELMAX,NELEM
00132       INTEGER, INTENT(INOUT) :: NCP ! CAN BE CHANGED
00133       INTEGER, INTENT(INOUT) :: CTRLSC(NCP) ! CAN BE CHANGED
00134       INTEGER, INTENT(IN) :: IKLE(NELMAX,*)
00135       INTEGER, INTENT(IN) :: IFABOR(NELMAX,*)
00136       DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*),TPS,DT
00137       DOUBLE PRECISION, INTENT(IN) :: XEL(NELMAX,*),YEL(NELMAX,*)
00138       LOGICAL, INTENT(IN) :: INFO,COMFLU,CUMFLO
00139       TYPE(BIEF_OBJ), INTENT(IN)    :: HPROP,S,U,V,H
00140       TYPE(BIEF_OBJ), INTENT(INOUT) :: BM1,BM2,T1,MSKSEC,CV1
00141       TYPE(BIEF_MESH) :: MESH
00142 !
00143 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00144 !
00145       INTEGER NSEMAX,ERR
00146       PARAMETER(NSEMAX=200)
00147 !
00148       INTEGER IELEM,IEL,I1,I2,I3,J1,J2,J3,ELBEST,IGBEST,ILBEST,IELMH
00149       INTEGER ILPREC,ISEG,ISEC,NSEC,PT,DEP,ARR,IELMU
00150 !
00151       DOUBLE PRECISION DIST,DIST1,DIST2,DIST3
00152       DOUBLE PRECISION H1,H2,X1,Y1,X2,Y2,UN1,UN2,NX,NY,SUR6
00153 !
00154       LOGICAL OK
00155 !
00156       DOUBLE PRECISION, ALLOCATABLE :: FLX(:),VOLNEG(:),VOLPOS(:)
00157       INTEGER, ALLOCATABLE :: NSEG(:),LISTE(:,:,:)
00158 !
00159       LOGICAL DEJA
00160       DATA DEJA/.FALSE./
00161       LOGICAL :: OLD_METHOD=.FALSE.
00162 !
00163       SAVE LISTE,DEJA,NSEG,VOLNEG,VOLPOS,FLX,OLD_METHOD
00164 !
00165 !-----------------------------------------------------------------------
00166 !
00167       SUR6=1.D0/6.D0
00168       NSEC = NCP/2
00169 !
00170 !  SEARCHES PATHS JOINING POINT PAIRS:
00171 !
00172       IF(.NOT.DEJA) THEN
00173 !
00174 !     DYNAMIC ALLOCATION OF FLX, VOLNEG, VOLPOS
00175 !
00176       ALLOCATE(FLX(NCP)             ,STAT=ERR)
00177       ALLOCATE(VOLNEG(NCP)          ,STAT=ERR)
00178       ALLOCATE(VOLPOS(NCP)          ,STAT=ERR)
00179       ALLOCATE(NSEG(NCP)            ,STAT=ERR)
00180       ALLOCATE(LISTE(NCP,NSEMAX,2)  ,STAT=ERR)
00181 !
00182       IF(ERR.NE.0) THEN
00183         IF(LNG.EQ.1) WRITE(LU,100) ERR
00184         IF(LNG.EQ.2) WRITE(LU,200) ERR
00185 100     FORMAT(1X,'FLUSEC : ERREUR A L''ALLOCATION DE MEMOIRE : ',/,1X,
00186      &            'CODE D''ERREUR : ',1I6)
00187 200     FORMAT(1X,'FLUSEC: ERROR DURING ALLOCATION OF MEMORY: ',/,1X,
00188      &            'ERROR CODE: ',1I6)
00189       ENDIF
00190 !
00191       IF (.NOT.ALLOCATED(CHAIN)) OLD_METHOD=.TRUE.
00192       DO ISEC =1,NSEC
00193 !!JAJ #### IN THE SERIAL CASE, OR "CLASSICAL" IN PARALLEL,
00194 !     FOLLOW THE ALGORITHM OF FINDING SEGMENT CHAINS
00195 !
00196 ! NOTE: IF YOU CHANGE THE ALGORITHM, CHANGE IT IN PARTEL AS WELL
00197 !
00198         IF (NCSIZE.LE.1 .OR. OLD_METHOD) THEN
00199 !
00200         DEP = CTRLSC(1+2*(ISEC-1))
00201         ARR = CTRLSC(2+2*(ISEC-1))
00202         VOLNEG(ISEC)=0.D0
00203         VOLPOS(ISEC)=0.D0
00204         IF(NCSIZE.GT.1) THEN
00205           DEP=GLOBAL_TO_LOCAL_POINT(DEP,MESH)
00206           ARR=GLOBAL_TO_LOCAL_POINT(ARR,MESH)
00207           IF(DEP.EQ.0.AND.ARR.EQ.0) THEN
00208             NSEG(ISEC)=0
00209             CYCLE
00210           ENDIF
00211           IF((DEP.EQ.0.AND.ARR.NE.0).OR.(DEP.NE.0.AND.ARR.EQ.0)) THEN
00212             NSEG(ISEC)=-1
00213             CYCLE
00214           ENDIF
00215         ENDIF
00216         PT = DEP
00217         ISEG = 0
00218         DIST=(X(DEP)-X(ARR))**2+(Y(DEP)-Y(ARR))**2
00219 10      CONTINUE
00220 !
00221         DO IELEM =1,NELEM
00222 !
00223           I1 = IKLE(IELEM,1)
00224           I2 = IKLE(IELEM,2)
00225           I3 = IKLE(IELEM,3)
00226 !         IF THE ELEMENT CONTAINS THE SELECTED NODE:
00227           IF(PT.EQ.I1.OR.PT.EQ.I2.OR.PT.EQ.I3) THEN
00228             DIST1 = (X(I1)-X(ARR))**2 + (Y(I1)-Y(ARR))**2
00229             DIST2 = (X(I2)-X(ARR))**2 + (Y(I2)-Y(ARR))**2
00230             DIST3 = (X(I3)-X(ARR))**2 + (Y(I3)-Y(ARR))**2
00231             IF(DIST1.LT.DIST) THEN
00232               DIST = DIST1
00233               ELBEST = IELEM
00234               IGBEST = I1
00235               ILBEST = 1
00236               IF(I1.EQ.PT) ILPREC = 1
00237               IF(I2.EQ.PT) ILPREC = 2
00238               IF(I3.EQ.PT) ILPREC = 3
00239             ENDIF
00240             IF(DIST2.LT.DIST) THEN
00241               DIST = DIST2
00242               ELBEST = IELEM
00243               IGBEST = I2
00244               ILBEST = 2
00245               IF(I1.EQ.PT) ILPREC = 1
00246               IF(I2.EQ.PT) ILPREC = 2
00247               IF(I3.EQ.PT) ILPREC = 3
00248             ENDIF
00249             IF(DIST3.LT.DIST) THEN
00250               DIST = DIST3
00251               ELBEST = IELEM
00252               IGBEST = I3
00253               ILBEST = 3
00254               IF(I1.EQ.PT) ILPREC = 1
00255               IF(I2.EQ.PT) ILPREC = 2
00256               IF(I3.EQ.PT) ILPREC = 3
00257             ENDIF
00258           ENDIF
00259 !
00260         ENDDO ! IELEM
00261 !
00262         IF(IGBEST.EQ.PT) THEN
00263           IF(LNG.EQ.1) WRITE(LU,32)
00264           IF(LNG.EQ.2) WRITE(LU,33)
00265 32        FORMAT(1X,'FLUSEC_TELEMAC2D : BLOCAGE DE L''ALGORITHME')
00266 33        FORMAT(1X,'FLUSEC_TELEMAC2D : ALGORITHM FAILED')
00267           CALL PLANTE(1)
00268           STOP
00269         ELSE
00270           PT = IGBEST
00271           ISEG = ISEG + 1
00272           IF(ISEG.GT.NSEMAX) THEN
00273             IF(LNG.EQ.1) THEN
00274               WRITE(LU,*) 'FLUSEC_TELEMAC2D : TROP DE SEGMENTS DANS UNE'
00275               WRITE(LU,*) '                   SECTION. AUGMENTER NSEMAX'
00276             ENDIF
00277             IF(LNG.EQ.2) THEN
00278               WRITE(LU,*) 'FLUSEC_TELEMAC2D: TOO MANY SEGMENTS IN A   '
00279               WRITE(LU,*) '                  SECTION. INCREASE  NSEMAX'
00280             ENDIF
00281             CALL PLANTE(1)
00282             STOP
00283           ENDIF
00284           LISTE(ISEC,ISEG,1) = IKLE(ELBEST,ILPREC)
00285           LISTE(ISEC,ISEG,2) = IKLE(ELBEST,ILBEST)
00286           IF(IGBEST.NE.ARR) GO TO 10
00287         ENDIF
00288 !
00289         NSEG(ISEC) = ISEG
00290 !
00291 !JAJ #### THIS PART TO BE DONE IN THE PARALLEL CASE; FILL LISTE
00292 ! WITH READY SEGMENT CHAINS PROVIDED BY PARTEL: SEE READ_SECTIONS
00293 ! NOTE: FUTURE OPTIMISATION - USE CHAIN STRUCTURE IN THE WHOLE ROUTINE
00294 !
00295         ELSE
00296 !
00297           FLX(ISEC)=0.0D0
00298           VOLNEG(ISEC)=0.0D0
00299           VOLPOS(ISEC)=0.0D0
00300           NSEG(ISEC) = CHAIN(ISEC)%NSEG
00301           LISTE(ISEC,:,:)=0
00302           DO ISEG=1,NSEG(ISEC)
00303             LISTE(ISEC,ISEG,1) = CHAIN(ISEC)%LISTE(ISEG,1)
00304             LISTE(ISEC,ISEG,2) = CHAIN(ISEC)%LISTE(ISEG,2)
00305           END DO
00306         ENDIF
00307 !
00308 !       IF COMPATIBLE COMPUTATION OF FLUXES=YES
00309         IF(COMFLU.AND.NSEG(ISEC).GE.1) THEN
00310 !
00311 !       LOOKING AT ALL ELEMENTS TOUCHING THE SECTION WITH 2 POINTS
00312 !       MARKING THEM WITH 1 ON ONE SIDE AND -1 ON THE OTHER SIDE
00313 !
00314         DO IEL=1,NELEM
00315           MSKSEC%ADR(ISEC)%P%R(IEL)=0.D0
00316           J1=IKLE(IEL,1)
00317           J2=IKLE(IEL,2)
00318           J3=IKLE(IEL,3)
00319           DO ISEG=1,NSEG(ISEC)
00320             I1 = LISTE(ISEC,ISEG,1)
00321             I2 = LISTE(ISEC,ISEG,2)
00322 !           LEFT SIDE
00323             IF    ( (J1.EQ.I1.AND.J2.EQ.I2) .OR.
00324      &              (J2.EQ.I1.AND.J3.EQ.I2) .OR.
00325      &              (J3.EQ.I1.AND.J1.EQ.I2)      ) THEN
00326               MSKSEC%ADR(ISEC)%P%R(IEL)=1.D0
00327 !           RIGHT SIDE
00328             ELSEIF( (J1.EQ.I2.AND.J2.EQ.I1) .OR.
00329      &              (J2.EQ.I2.AND.J3.EQ.I1) .OR.
00330      &              (J3.EQ.I2.AND.J1.EQ.I1)      ) THEN
00331               MSKSEC%ADR(ISEC)%P%R(IEL)=-1.D0
00332             ENDIF
00333           ENDDO
00334         ENDDO
00335 !
00336 !       OTHER TRIANGLES WITH ONLY 1 POINT TOUCHING THE SECTION
00337 !       LOOKING AT NEIGHBOURS TO FIND THEIR SIDE
00338 !
00339         DO
00340           OK=.TRUE.
00341           DO IEL=1,NELEM
00342             J1=IKLE(IEL,1)
00343             J2=IKLE(IEL,2)
00344             J3=IKLE(IEL,3)
00345             DO ISEG=1,NSEG(ISEC)
00346               I1 = LISTE(ISEC,ISEG,1)
00347               I2 = LISTE(ISEC,ISEG,2)
00348               IF((J1.EQ.I1.OR.J2.EQ.I1.OR.J3.EQ.I1.OR.
00349      &            J1.EQ.I2.OR.J2.EQ.I2.OR.J3.EQ.I2).AND.
00350      &            ABS(MSKSEC%ADR(ISEC)%P%R(IEL)).LT.0.5D0) THEN
00351 !               LOOKING AT NEIGHBOURS
00352                 IF(IFABOR(IEL,1).GT.0) THEN
00353                   IELEM=IFABOR(IEL,1)
00354                   IF(ABS(MSKSEC%ADR(ISEC)%P%R(IELEM)).GT.0.5D0) THEN
00355                     MSKSEC%ADR(ISEC)%P%R(IEL)=
00356      &                           MSKSEC%ADR(ISEC)%P%R(IELEM)
00357                   ENDIF
00358                 ENDIF
00359                 IF(IFABOR(IEL,2).GT.0) THEN
00360                   IELEM=IFABOR(IEL,2)
00361                   IF(ABS(MSKSEC%ADR(ISEC)%P%R(IELEM)).GT.0.5D0) THEN
00362                     MSKSEC%ADR(ISEC)%P%R(IEL)=
00363      &                           MSKSEC%ADR(ISEC)%P%R(IELEM)
00364                   ENDIF
00365                 ENDIF
00366                 IF(IFABOR(IEL,3).GT.0) THEN
00367                   IELEM=IFABOR(IEL,3)
00368                   IF(ABS(MSKSEC%ADR(ISEC)%P%R(IELEM)).GT.0.5D0) THEN
00369                     MSKSEC%ADR(ISEC)%P%R(IEL)=
00370      &                           MSKSEC%ADR(ISEC)%P%R(IELEM)
00371                   ENDIF
00372                 ENDIF
00373 !               CORRECTION JMH 11/01/2008
00374 !               IF(MSKSEC%ADR(ISEC)%P%R(IEL).LT.0.5D0) OK=.FALSE.
00375                 IF(ABS(MSKSEC%ADR(ISEC)%P%R(IEL)).LT.0.5D0) OK=.FALSE.
00376               ENDIF
00377             ENDDO
00378           ENDDO
00379           IF(OK) EXIT
00380         ENDDO
00381 !
00382 !       IF(COMFLU.AND.NSEG(ISEC).GE.1) THEN
00383         ENDIF
00384 !
00385       ENDDO ! ISEC
00386 !
00387 !     IF(.NOT.DEJA) THEN
00388       ENDIF
00389 !
00390 !-----------------------------------------------------------------------
00391 !
00392       DEJA = .TRUE.
00393 !
00394 !-----------------------------------------------------------------------
00395 !
00396       IELMH=HPROP%ELM
00397       IELMU=U%ELM
00398 !
00399       DO ISEC = 1 , NSEC
00400       FLX(ISEC) = 0.D0
00401       IF(NSEG(ISEC).GE.1) THEN
00402 !
00403       IF(COMFLU) THEN
00404 !
00405 !     COMPUTING THE FLUX AS IN THE CONTINUITY EQUATION
00406 !     (HOWEVER IMPLICITATION SHOULD PERHAPS BE ALSO USED)
00407 !
00408       CALL MATRIX(BM1,'M=N     ','MATFGR         X',IELMH,IELMU,
00409      &            1.D0,HPROP,S,S,S,S,S,MESH,.TRUE.,MSKSEC%ADR(ISEC)%P)
00410       CALL MATRIX(BM2,'M=N     ','MATFGR         Y',IELMH,IELMU,
00411      &            1.D0,HPROP,S,S,S,S,S,MESH,.TRUE.,MSKSEC%ADR(ISEC)%P)
00412 !
00413       CALL MATVEC( 'X=AY    ',CV1,BM1,U,0.D0,MESH)
00414       CALL MATVEC( 'X=X+AY  ',CV1,BM2,V,0.D0,MESH)
00415 !
00416 !     SUMMING CV1 FOR ALL POINTS OF THE SECTION, THIS IS THE FLUX !
00417 !     (OBTAINED BY CONTINUITY EQUATION AND AN INTEGRATION BY PARTS)
00418 !
00419       DO ISEG = 1 , NSEG(ISEC)
00420         I1   = LISTE(ISEC,ISEG,1)
00421         FLX(ISEC) = FLX(ISEC) + CV1%R(I1)
00422       ENDDO
00423 !     LAST SEGMENT, ADDING THE LAST POINT
00424       I2   = LISTE(ISEC,NSEG(ISEC),2)
00425       FLX(ISEC) = FLX(ISEC) + CV1%R(I2)
00426 !
00427 !     WHEN BOTH UPWIND AND DOWNSTREAM ELEMENTS ARE TAKEN INTO ACCOUNT
00428 !     WITH DIFFERENT SIGNS, WE GET TWICE THE FLUX
00429 !
00430       FLX(ISEC)=FLX(ISEC)*0.5D0
00431 !
00432       ELSE
00433 !
00434 !       COMPUTING THE FLUX DIRECTLY, REGARDLESS OF THE WEAK FORM
00435 !       OF THE IMPERMEABILITY CONDITION
00436 !
00437         DO ISEG = 1 , NSEG(ISEC)
00438           I1 = LISTE(ISEC,ISEG,1)
00439           I2 = LISTE(ISEC,ISEG,2)
00440           X1 = X(I1)
00441           X2 = X(I2)
00442           Y1 = Y(I1)
00443           Y2 = Y(I2)
00444           H1 = H%R(I1)
00445           H2 = H%R(I2)
00446           NX = Y1-Y2
00447           NY = X2-X1
00448           UN1= U%R(I1)*NX + V%R(I1)*NY
00449           UN2= U%R(I2)*NX + V%R(I2)*NY
00450           FLX(ISEC) = FLX(ISEC) + ((H1+H2)*(UN1+UN2)+H2*UN2+H1*UN1)*SUR6
00451         ENDDO
00452 !
00453       ENDIF
00454 !
00455       IF(FLX(ISEC).GT.0.D0) THEN
00456         VOLPOS(ISEC) = VOLPOS(ISEC) + FLX(ISEC)*DT
00457       ELSE
00458         VOLNEG(ISEC) = VOLNEG(ISEC) + FLX(ISEC)*DT
00459       ENDIF
00460 !
00461 !     IF(NSEG(ISEC).GT.1)...
00462       ENDIF
00463 !
00464       ENDDO ! ISEC
00465 !
00466 !-----------------------------------------------------------------------
00467 !
00468 !     PRINTING THE RESULTS / !JAJ HERE ALLREDUCES FOR VALUES
00469 !
00470       CALL FLUXPR_TELEMAC2D
00471      &  (NSEC,CTRLSC,FLX,VOLNEG,VOLPOS,INFO,TPS,
00472      &   NSEG,NCSIZE,CUMFLO)
00473 !
00474 !-----------------------------------------------------------------------
00475 !
00476       RETURN
00477       END

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