bedload_nerbed_vf.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\sisyphe\bedload_nerbed_vf.f
00002 !
00079                      SUBROUTINE BEDLOAD_NERBED_VF !
00080 !                    ******************************
00081 !
00082      &(MESH,LIEBOR,KSORT,ELAY,V2DPAR,QSX,QSY,AVA,NPOIN,NSEG,NPTFR,
00083      & DT,QS,T1,T2,T3,BREACH,CSF_SABLE,NUBO,VNOIN)
00084 !
00085 !***********************************************************************
00086 ! SISYPHE   V7P0                                   03/06/2014
00087 !***********************************************************************
00088 !
00089 !
00090 !
00091 !
00092 !
00093 !
00094 !
00095 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00096 !| AVA            |-->| PERCENT AVAILABLE
00097 !| BREACH         |<->| INDICATOR FOR NON ERODIBLE BED (FINITE VOLUMES SCHEMES)
00098 !| DT             |-->| TIME STEP
00099 !| ELAY           |<->| THICKNESS OF SURFACE LAYER
00100 !| KSORT          |-->| CONVENTION FOR FREE OUTPUT
00101 !| LIEBOR         |<->| PHYSICAL BOUNDARY CONDITIONS FOR BED EVOLUTION
00102 !| MESH           |<->| MESH STRUCTURE
00103 !| NPOIN          |-->| NUMBER OF POINTS
00104 !| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
00105 !| NSEG           |-->| NUMBER OF SEGMENTS PER CONTROL SECTION
00106 !| NUBO           |-->| GLOBAL NUMBER OF EDGE EXTREMITIES
00107 !| QS             |<->| BEDLOAD TRANSPORT RATE
00108 !| QSX            |-->| SOLID DISCHARGE X
00109 !| QSY            |-->| SOLID DISCHARGE Y
00110 !| T1             |<->| WORK BIEF_OBJ STRUCTURE
00111 !| T2             |<->| WORK BIEF_OBJ STRUCTURE
00112 !| T3             |<->| WORK BIEF_OBJ STRUCTURE
00113 !| V2DPAR         |-->| INTEGRAL OF TEST FUNCTIONS, ASSEMBLED IN PARALLEL
00114 !| VNOIN          |-->| OUTWARD UNIT NORMALS
00115 !| CSF_SABLE      |-->| VOLUME CONCENTRATION OF SAND (1-POROSITY)
00116 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00117 !
00118       USE INTERFACE_SISYPHE, EX_BEDLOAD_NERBED_VF => BEDLOAD_NERBED_VF
00119       USE BIEF
00120       IMPLICIT NONE
00121       INTEGER LNG,LU
00122       COMMON/INFO/LNG,LU
00123 !
00124 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00125 !
00126       TYPE(BIEF_MESH),  INTENT(INOUT) :: MESH
00127       TYPE(BIEF_OBJ),   INTENT(IN)    :: LIEBOR
00128       TYPE(BIEF_OBJ),   INTENT(IN)    :: QSX, QSY
00129       INTEGER,          INTENT(IN)    :: NPOIN, NSEG, NPTFR,KSORT
00130       DOUBLE PRECISION, INTENT(IN)    :: DT
00131       TYPE(BIEF_OBJ),   INTENT(INOUT) :: QS, T1, T2, T3
00132       TYPE(BIEF_OBJ),   INTENT(INOUT) :: BREACH
00133       DOUBLE PRECISION, INTENT(IN)    :: ELAY(NPOIN),V2DPAR(NPOIN)
00134       DOUBLE PRECISION, INTENT(IN)    :: AVA(NPOIN), CSF_SABLE
00135 ! RA
00136       INTEGER, INTENT(IN)             :: NUBO(2,NSEG)
00137       DOUBLE PRECISION, INTENT(IN)    :: VNOIN(3,NSEG)
00138 !
00139 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00140 !
00141       INTEGER          :: I, K
00142       INTEGER          :: IEL, IEL1, IEL2, ISEGIN
00143       DOUBLE PRECISION :: QSP1, QSP2, QSPC
00144       DOUBLE PRECISION :: XN, YN, TEMP,PROD_SCAL
00145       DOUBLE PRECISION :: VNOIN1, VNOIN2, RNORM
00146 !
00147 !======================================================================!
00148 !======================================================================!
00149 !                               PROGRAM                                !
00150 !======================================================================!
00151 !======================================================================!
00152 !
00153       ! ****************** !
00154       ! I - INITIALISES    !
00155       ! ****************** !
00156 !
00157       ! BREACH INDICATES IF NON ERODABLE BED WILL BE REACHED
00158       ! DURING TIME STEP FOR THIS POINT
00159       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00160 !
00161 !     GIVE T1 AND T2 THE SAME STRUCTURE AS QS
00162 !
00163       CALL CPSTVC(QS,T1)
00164       CALL CPSTVC(QS,T2)
00165 !
00166       DO IEL = 1, NPOIN
00167         BREACH%I(IEL) = 0
00168         T1%R(IEL)=0.D0
00169         T2%R(IEL)=0.D0
00170       ENDDO
00171 !
00172       ! ************************************************* !
00173       ! II - DETERMINES THE OUTGOING FLUX FOR EACH CELL   ! (_IMP_)
00174       ! ************************************************* !
00175       ! THE PRINCIPLE IS THAT QS IS CALCULATED FOR EACH SEGMENT AS
00176       ! HALF THE SUM OF THE QS AT THE CENTERS OF THE ELEMENTS WHICH
00177       ! SEGMENT FORMS THE BOUNDARY. IT IS THEN PROJECTED ON THE NORMAL
00178       ! TO THE SEGMENT, MULTIPLIED BY THE LENGTH OF THE SEGMENT, AND
00179       ! THIS FLUX TERM IS ADDED (OR SUBTRACTED) TO THE TWO ELEMENTS.
00180 !
00181       DO ISEGIN = 1, NSEG
00182 !
00183         IEL1 = NUBO(1,ISEGIN)
00184         IEL2 = NUBO(2,ISEGIN)
00185 !
00186         ! II.1 - SEGMENT LENGTH (RNORM)
00187         ! ----------------------------------
00188         VNOIN1 = VNOIN(1,ISEGIN)
00189         VNOIN2 = VNOIN(2,ISEGIN)
00190         RNORM  = VNOIN(3,ISEGIN)
00191         PROD_SCAL= (MESH%X%R(IEL2)-MESH%X%R(IEL1))*VNOIN1+
00192      &             (MESH%Y%R(IEL2)-MESH%Y%R(IEL1))*VNOIN2
00193         IF(PROD_SCAL.LT.0.D0)THEN
00194           IEL1 = NUBO(2,ISEGIN)
00195           IEL2 = NUBO(1,ISEGIN)
00196         ENDIF
00197 !
00198         ! II.2 - PROJECTS QS FOR THE SEGMENT ONTO THE SEGMENT NORMAL
00199         ! ------------------------------------------------------------
00200         QSP1 = VNOIN1*QSX%R(IEL1) + VNOIN2*QSY%R(IEL1)
00201         QSP2 = VNOIN1*QSX%R(IEL2) + VNOIN2*QSY%R(IEL2)
00202         QSPC = (QSP1+QSP2)*0.5D0
00203 !
00204         ! II.3 - QS SUCH AS THE OUTGOING FLUX IS MAXIMUM
00205         ! ----------------------------------------------
00206         T1%R(IEL1) = T1%R(IEL1) + RNORM*MAX(QSPC,QSP1,0.D0)
00207         T1%R(IEL2) = T1%R(IEL2) - RNORM*MIN(QSPC,QSP2,0.D0)
00208 !
00209         IF(QSPC > 0.D0) THEN
00210           T2%R(IEL1) = T2%R(IEL1) + RNORM*QSP1
00211         ELSEIF(QSPC < 0.D0) THEN
00212           T2%R(IEL2) = T2%R(IEL2) - RNORM*QSP2
00213         ENDIF
00214 !
00215       ENDDO
00216 !
00217       ! ************************************** !
00218       ! III - LOOP ON THE BOUNDARY NODES       !
00219       ! ************************************** !
00220 !
00221       DO K = 1, NPTFR
00222         IEL = MESH%NBOR%I(K)
00223 !
00224         ! III.1 - FREE EVOLUTION: SEDIMENTS ARE FREE TO LEAVE
00225         ! ---------------------------------------------------------
00226         IF (LIEBOR%I(K) == KSORT) THEN
00227 !
00228           ! XNEBOR (*+NPTFR) AND YNEBOR (*+NPTFR)
00229           ! CONTAIN THE VECTOR NORMAL TO A BOUNDARY NODE
00230           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00231           XN   = MESH%XNEBOR%R(K+NPTFR)
00232           YN   = MESH%YNEBOR%R(K+NPTFR)
00233           TEMP = QSX%R(IEL)*XN + QSY%R(IEL)*YN
00234           IF (TEMP > 0.D0) THEN
00235             T1%R(IEL) = T1%R(IEL) + TEMP
00236             T2%R(IEL) = T2%R(IEL) + TEMP
00237           ENDIF
00238 !
00239         ENDIF
00240 !
00241         ! III.2 - FOR A SOLID BOUNDARY: NOTHING TO PROGRAM
00242         !         BECAUSE THE SEDIMENT FLUX IS ZERO HERE
00243         !         FOR IMPOSED EVOLUTION : SEE BEDLOAD_SOLVS_VF.F
00244         ! --------------------------------------------------------
00245       ENDDO
00246 !
00247       IF(NCSIZE > 1) THEN
00248         CALL PARCOM(T1, 2, MESH)
00249         CALL PARCOM(T2, 2, MESH)
00250       ENDIF
00251 !
00252       ! ************************************************ !
00253       ! IV - COMPUTES THE MAXIMUM FLUX AUTHORISED PER CELL!
00254       ! ************************************************ !
00255 !
00256       DO I = 1, NPOIN
00257 !
00258         T3%R(I)=ELAY(I)*V2DPAR(I)*AVA(I)* CSF_SABLE/DT
00259         IF (T3%R(I) < 0.D0) T3%R(I) = 0.D0
00260 !
00261         ! IF THE OUTGOING FLUX IS TOO LARGE, QS IS CAPPED AT THE NODE
00262         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00263         IF(T1%R(I) > T3%R(I)) THEN
00264           BREACH%I(I) = 1
00265           IF(T2%R(I) > T3%R(I)) THEN
00266             QS%R(I) = QS%R(I)*T3%R(I)/T2%R(I)
00267           ENDIF
00268         ENDIF
00269 !
00270       ENDDO
00271 !
00272 !======================================================================!
00273 !======================================================================!
00274 !
00275       RETURN
00276       END

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