flusec_sisyphe.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\sisyphe\flusec_sisyphe.f
00002 !
00075                      SUBROUTINE FLUSEC_SISYPHE
00076 !                    *************************
00077 !
00078      &(U,V,H,QSXC,QSYC,CHARR,QSXS,QSYS,SUSP,
00079      & IKLE,NELMAX,NELEM,X,Y,DT,NCP,CTRLSC,INFO,TPS)
00080 !
00081 !***********************************************************************
00082 ! SISYPHE   V7P0                                   21/07/2011
00083 !***********************************************************************
00084 !
00085 !
00086 !
00087 !
00088 !
00089 !
00090 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00091 !| CHARR          |-->| LOGICAL, BEDLOAD
00092 !| CTRLSC         |-->| CONTROL SECTION
00093 !| DT             |-->| TIME STEP
00094 !| H              |-->| WATER DEPTH
00095 !| IKLE           |-->| CONNECTIVITY TABLE
00096 !| INFO           |-->| IF YES, PRINT
00097 !| NCP            |-->| TWO TIMES THE NUMBER OF CONTROL SECTIONS
00098 !| NELEM          |-->| NUMBER OF ELEMENTS
00099 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00100 !| QSXC           |<->| BEDLOAD TRANSPORT RATE X-DIRECTION
00101 !| QSXS           |<->| SUSPENSION TRANSPORT RATE X-DIRECTION
00102 !| QSYC           |<->| BEDLOAD TRANSPORT RATE Y-DIRECTION
00103 !| QSYS           |<->| SUSPENSION TRANSPORT RATE Y-DIRECTION
00104 !| SUSP           |-->| LOGICAL, SUSPENSION
00105 !| TPS            |-->| TEMPS
00106 !| U,V            |-->| VELOCITY FIELD COMPONENTS
00107 !| X,Y            |-->| NODES COORDINATES
00108 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00109 !
00110       USE BIEF
00111       USE DECLARATIONS_SISYPHE, ONLY: CHAIN,MESH
00112       USE INTERFACE_SISYPHE, EX_FLUSEC_SISYPHE => FLUSEC_SISYPHE
00113 !
00114       IMPLICIT NONE
00115       INTEGER LNG,LU
00116       COMMON/INFO/LNG,LU
00117 !
00118 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00119 !
00120       INTEGER, INTENT(IN)          :: NELMAX,NELEM,NCP
00121       INTEGER, INTENT(IN)          :: IKLE(NELMAX,*)
00122       INTEGER, INTENT(IN)          :: CTRLSC(NCP)
00123       DOUBLE PRECISION, INTENT(IN) :: X(*),Y(*),TPS,DT
00124       LOGICAL, INTENT(IN)          :: INFO,SUSP,CHARR
00125       TYPE(BIEF_OBJ), INTENT(IN)   :: U,V,H,QSXC,QSYC,QSXS,QSYS
00126 !
00127 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00128 !
00129       INTEGER NSEMAX,ERR
00130       PARAMETER(NSEMAX=100)
00131 !
00132       INTEGER IELEM,I1,I2,I3,ELBEST,IGBEST,ILBEST
00133       INTEGER ILPREC,ISEG,ISEC,NSEC,PT,DEP,ARR
00134 !
00135       DOUBLE PRECISION DIST,DIST1,DIST2,DIST3
00136       DOUBLE PRECISION H1,H2,X1,Y1,X2,Y2,UN1,UN2,NX,NY,SUR6
00137 !
00138       DOUBLE PRECISION, ALLOCATABLE :: FLX(:),VOLNEG(:),VOLPOS(:)
00139       DOUBLE PRECISION, ALLOCATABLE :: FLXS(:),FLXC(:)
00140       DOUBLE PRECISION, ALLOCATABLE :: VOLNEGS(:),VOLPOSS(:)
00141       DOUBLE PRECISION, ALLOCATABLE :: VOLNEGC(:),VOLPOSC(:)
00142       INTEGER, ALLOCATABLE :: NSEG(:),LISTE(:,:,:)
00143 !
00144       LOGICAL DEJA
00145       DATA DEJA/.FALSE./
00146       LOGICAL :: OLD_METHOD=.FALSE.
00147 !
00148       SAVE LISTE,DEJA,NSEG,VOLNEG,VOLPOS,FLX,FLXS,VOLNEGS,VOLPOSS
00149       SAVE FLXC,VOLNEGC,VOLPOSC,OLD_METHOD
00150 !
00151 !-----------------------------------------------------------------------
00152 !
00153       SUR6 = 1.D0/6.D0
00154       NSEC = NCP/2
00155 !
00156 !  LOOKS FOR WAYS THAT JOIN THE POINT COUPLES:
00157 !
00158       IF(.NOT.DEJA) THEN
00159 !
00160 !     DYNAMICALLY ALLOCATES FLX, VOLNEG, VOLPOS, ETC.
00161 !
00162       ALLOCATE(FLX(NSEC)           ,STAT=ERR)
00163       ALLOCATE(VOLNEG(NSEC)        ,STAT=ERR)
00164       ALLOCATE(VOLPOS(NSEC)        ,STAT=ERR)
00165       ALLOCATE(NSEG(NCP)           ,STAT=ERR)
00166       ALLOCATE(LISTE(NCP,NSEMAX,2) ,STAT=ERR)
00167 !     S FOR SUSPENSION, C FOR BEDLOAD
00168       ALLOCATE(FLXS(NSEC)          ,STAT=ERR)
00169       ALLOCATE(VOLNEGS(NSEC)       ,STAT=ERR)
00170       ALLOCATE(VOLPOSS(NSEC)       ,STAT=ERR)
00171       ALLOCATE(FLXC(NSEC)          ,STAT=ERR)
00172       ALLOCATE(VOLNEGC(NSEC)       ,STAT=ERR)
00173       ALLOCATE(VOLPOSC(NSEC)       ,STAT=ERR)
00174 !
00175       IF(ERR.NE.0) THEN
00176         IF(LNG.EQ.1) WRITE(LU,100) ERR
00177         IF(LNG.EQ.2) WRITE(LU,200) ERR
00178 100     FORMAT(1X,'FLUSEC : ERREUR A L''ALLOCATION DE MEMOIRE : ',/,1X,
00179      &            'CODE D''ERREUR : ',1I6)
00180 200     FORMAT(1X,'FLUSEC: ERROR DURING ALLOCATION OF MEMORY: ',/,1X,
00181      &            'ERROR CODE: ',1I6)
00182       ENDIF
00183 !
00184       IF (.NOT.ALLOCATED(CHAIN)) OLD_METHOD=.TRUE.
00185 !
00186       DO ISEC=1,NSEC
00187         FLX(ISEC)=0.0D0
00188         VOLNEG(ISEC) =0.D0
00189         VOLPOS(ISEC) =0.D0
00190         VOLNEGS(ISEC)=0.D0
00191         VOLPOSS(ISEC)=0.D0
00192         VOLNEGC(ISEC)=0.D0
00193         VOLPOSC(ISEC)=0.D0
00194       ENDDO
00195 !
00196       DO ISEC =1,NSEC
00197 !
00198 !!JAJ #### IN THE SERIAL CASE, OR "CLASSICAL" IN PARALLEL,
00199 !     FOLLOW THE ALGORITHM OF FINDING SEGMENT CHAINS
00200 !
00201 ! NOTE: IF YOU CHANGE THE ALGORITHM, CHANGE IT IN PARTEL AS WELL
00202 !
00203         IF (NCSIZE.LE.1 .OR. OLD_METHOD) THEN
00204 !
00205         DEP = CTRLSC(1+2*(ISEC-1))
00206         ARR = CTRLSC(2+2*(ISEC-1))
00207         IF(NCSIZE.GT.1) THEN
00208           DEP=GLOBAL_TO_LOCAL_POINT(DEP,MESH)
00209           ARR=GLOBAL_TO_LOCAL_POINT(ARR,MESH)
00210           IF(DEP.EQ.0.AND.ARR.EQ.0) THEN
00211             NSEG(ISEC)=0
00212             EXIT
00213           ENDIF
00214           IF((DEP.EQ.0.AND.ARR.NE.0).OR.(DEP.NE.0.AND.ARR.EQ.0)) THEN
00215             NSEG(ISEC)=-1
00216             EXIT
00217           ENDIF
00218         ENDIF
00219         PT = DEP
00220         ISEG = 0
00221         DIST=(X(DEP)-X(ARR))**2+(Y(DEP)-Y(ARR))**2
00222 10      CONTINUE
00223 !
00224         DO IELEM =1,NELEM
00225 !
00226           I1 = IKLE(IELEM,1)
00227           I2 = IKLE(IELEM,2)
00228           I3 = IKLE(IELEM,3)
00229 !         IF THE ELEMENT CONTAINS THE CURRENT POINT:
00230           IF(PT.EQ.I1.OR.PT.EQ.I2.OR.PT.EQ.I3) THEN
00231             DIST1 = (X(I1)-X(ARR))**2 + (Y(I1)-Y(ARR))**2
00232             DIST2 = (X(I2)-X(ARR))**2 + (Y(I2)-Y(ARR))**2
00233             DIST3 = (X(I3)-X(ARR))**2 + (Y(I3)-Y(ARR))**2
00234             IF(DIST1.LT.DIST) THEN
00235               DIST = DIST1
00236               ELBEST = IELEM
00237               IGBEST = I1
00238               ILBEST = 1
00239               IF(I1.EQ.PT) ILPREC = 1
00240               IF(I2.EQ.PT) ILPREC = 2
00241               IF(I3.EQ.PT) ILPREC = 3
00242             ENDIF
00243             IF(DIST2.LT.DIST) THEN
00244               DIST = DIST2
00245               ELBEST = IELEM
00246               IGBEST = I2
00247               ILBEST = 2
00248               IF(I1.EQ.PT) ILPREC = 1
00249               IF(I2.EQ.PT) ILPREC = 2
00250               IF(I3.EQ.PT) ILPREC = 3
00251             ENDIF
00252             IF(DIST3.LT.DIST) THEN
00253               DIST = DIST3
00254               ELBEST = IELEM
00255               IGBEST = I3
00256               ILBEST = 3
00257               IF(I1.EQ.PT) ILPREC = 1
00258               IF(I2.EQ.PT) ILPREC = 2
00259               IF(I3.EQ.PT) ILPREC = 3
00260             ENDIF
00261           ENDIF
00262 !
00263         ENDDO !IELEM
00264 !
00265         IF(IGBEST.EQ.PT) THEN
00266           IF(LNG.EQ.1) WRITE(LU,32)
00267           IF(LNG.EQ.2) WRITE(LU,33)
00268 32        FORMAT(1X,'FLUSEC : BLOCAGE DE L''ALGORITHME')
00269 33        FORMAT(1X,'FLUSEC : ALGORITHM FAILED')
00270           CALL PLANTE(1)
00271           STOP
00272         ELSE
00273           PT = IGBEST
00274           ISEG = ISEG + 1
00275           IF(ISEG.GT.NSEMAX) THEN
00276             IF(LNG.EQ.1) THEN
00277               WRITE(LU,*) 'FLUSEC : TROP DE SEGMENTS DANS UNE'
00278               WRITE(LU,*) '         SECTION. AUGMENTER NSEMAX'
00279             ENDIF
00280             IF(LNG.EQ.2) THEN
00281               WRITE(LU,*) 'FLUSEC: TOO MANY SEGMENTS IN A   '
00282               WRITE(LU,*) '        SECTION. INCREASE  NSEMAX'
00283             ENDIF
00284             CALL PLANTE(1)
00285             STOP
00286           ENDIF
00287           LISTE(ISEC,ISEG,1) = IKLE(ELBEST,ILPREC)
00288           LISTE(ISEC,ISEG,2) = IKLE(ELBEST,ILBEST)
00289           IF(IGBEST.NE.ARR) GO TO 10
00290         ENDIF
00291 !
00292         NSEG(ISEC) = ISEG
00293 !
00294 !JAJ #### THIS PART TO BE DONE IN THE PARALLEL CASE; FILL LISTE
00295 ! WITH READY SEGMENT CHAINS PROVIDED BY PARTEL: SEE READ_SECTIONS
00296 ! NOTE: FUTURE OPTIMISATION - USE CHAIN STRUCTURE IN THE WHOLE ROUTINE
00297 !
00298         ELSE
00299 !
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           WRITE(LU,*) 'CHAIN@SISYPHE -> LISTE@SISYPHE:'
00307           WRITE(LU,*) 'ISEC,NSEG(ISEC): ',ISEC,NSEG(ISEC)
00308           DO ISEG=1,NSEG(ISEC)
00309             WRITE(LU,*) LISTE(ISEC,ISEG,:)
00310           END DO
00311         ENDIF
00312 !
00313       ENDDO
00314 !
00315 !     IF(.NOT.DEJA) THEN
00316       ENDIF
00317 !
00318 !-----------------------------------------------------------------------
00319 !
00320       DEJA = .TRUE.
00321 !
00322 !-----------------------------------------------------------------------
00323 !
00324       DO ISEC = 1 , NSEC
00325 !
00326       FLX(ISEC)  = 0.D0
00327       FLXS(ISEC) = 0.D0
00328       FLXC(ISEC) = 0.D0
00329 !
00330       IF(NSEG(ISEC).GE.1) THEN
00331 !
00332 !     COMPUTES THE FLUX DIRECTLY, REGARDLESS OF THE WEAK FORM
00333 !     OF THE IMPERMEABILITY CONDITION
00334 !
00335       DO ISEG = 1 , NSEG(ISEC)
00336         I1 = LISTE(ISEC,ISEG,1)
00337         I2 = LISTE(ISEC,ISEG,2)
00338         X1 = X(I1)
00339         X2 = X(I2)
00340         Y1 = Y(I1)
00341         Y2 = Y(I2)
00342         H1 = H%R(I1)
00343         H2 = H%R(I2)
00344         NX = Y1-Y2
00345         NY = X2-X1
00346         UN1= U%R(I1)*NX + V%R(I1)*NY
00347         UN2= U%R(I2)*NX + V%R(I2)*NY
00348         FLX(ISEC) = FLX(ISEC) + ((H1+H2)*(UN1+UN2)+H2*UN2+H1*UN1)*SUR6
00349         IF(SUSP) THEN
00350           UN1= QSXS%R(I1)*NX + QSYS%R(I1)*NY
00351           UN2= QSXS%R(I2)*NX + QSYS%R(I2)*NY
00352           FLXS(ISEC) = FLXS(ISEC) + 0.5D0*(UN1+UN2)
00353         ENDIF
00354         IF(CHARR) THEN
00355           UN1= QSXC%R(I1)*NX + QSYC%R(I1)*NY
00356           UN2= QSXC%R(I2)*NX + QSYC%R(I2)*NY
00357           FLXC(ISEC) = FLXC(ISEC) + 0.5D0*(UN1+UN2)
00358         ENDIF
00359       ENDDO
00360 !
00361       IF(FLX(ISEC).GT.0.D0) THEN
00362         VOLPOS(ISEC) = VOLPOS(ISEC) + FLX(ISEC)*DT
00363       ELSE
00364         VOLNEG(ISEC) = VOLNEG(ISEC) + FLX(ISEC)*DT
00365       ENDIF
00366 !
00367       IF(SUSP) THEN
00368         IF(FLXS(ISEC).GT.0.D0) THEN
00369           VOLPOSS(ISEC) = VOLPOSS(ISEC) + FLXS(ISEC)*DT
00370         ELSE
00371           VOLNEGS(ISEC) = VOLNEGS(ISEC) + FLXS(ISEC)*DT
00372         ENDIF
00373       ENDIF
00374 !
00375       IF(CHARR) THEN
00376         IF(FLXC(ISEC).GT.0.D0) THEN
00377           VOLPOSC(ISEC) = VOLPOSC(ISEC) + FLXC(ISEC)*DT
00378         ELSE
00379           VOLNEGC(ISEC) = VOLNEGC(ISEC) + FLXC(ISEC)*DT
00380         ENDIF
00381       ENDIF
00382 !
00383 !     IF(NSEG(ISEC).GT.1)...
00384       ENDIF
00385 !
00386       ENDDO
00387 !
00388 !-----------------------------------------------------------------------
00389 !
00390 !     PRINTS THE RESULTS / !JAJ HERE ALLREDUCES FOR VALUES
00391 !
00392       CALL FLUXPR_SISYPHE(NSEC,CTRLSC,FLX,VOLNEG,VOLPOS,INFO,TPS,
00393      &                    NSEG,NCSIZE,
00394      &                    FLXS,VOLNEGS,VOLPOSS,SUSP,
00395      &                    FLXC,VOLNEGC,VOLPOSC,CHARR)
00396 !
00397 !-----------------------------------------------------------------------
00398 !
00399 !      WRITE(LU,*) '-> LEAVING FLUSEC_SISYPHE'
00400       RETURN
00401       END

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