bedload_solvs_vf.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\sisyphe\bedload_solvs_vf.f
00002 !
00084                      SUBROUTINE BEDLOAD_SOLVS_VF
00085 !                    ***************************
00086 !
00087      &(MESH,QSX,QSY,LIMTEC,UNSV2D,EBOR,BREACH,NSEG,NPTFR,NPOIN,
00088      & KENT,KDIR,KDDL,DT,T10,ZFCL,FLUX,CSF_SABLE,FLBCLA,AVA,LIQBOR,QBOR,
00089      & NUBO,VNOIN)
00090 !
00091 !***********************************************************************
00092 ! SISYPHE   V7P0                                      03/06/2014
00093 !***********************************************************************
00094 !
00095 !
00096 !
00097 !
00098 !
00099 !
00100 !
00101 !
00102 !
00103 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00104 !| AVA            |-->| PERCENTAGE OF CLASS IN SEDIMENT
00105 !| BREACH         |<->| INDICATOR FOR NON ERODIBLE BED
00106 !| DT             |-->| TIME STEP
00107 !| EBOR           |<->| BOUNDARY CONDITION FOR BED EVOLUTION (DIRICHLET)
00108 !| FLBCLA         |<->| FLUXES AT BOUNDARY FOR THE CLASS
00109 !| FLUX           |<->| SEDIMENT FLUX
00110 !| KDIR           |-->| CONVENTION FOR DIRICHLET VALUE
00111 !| KDDL           |-->| CONVENTION FOR DEGREE OF FREEDOM
00112 !| KENT           |-->| CONVENTION FOR PRESCRIBED VALUE AT ENTRANCE
00113 !| LIMTEC         |-->| TECHNICAL BOUNDARY CONDITIONS FOR BED EVOLUTION
00114 !| LIQBOR         |-->| TYPE OF BOUNDARY CONDITIONS FOR BEDLOAD DISCHARGE
00115 !| MESH           |<->| MESH STRUCTURE
00116 !| NPOIN          |-->| NUMBER OF POINTS
00117 !| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
00118 !| NSEG           |-->| NUMBER OF SEGMENTS
00119 !| NUBO           |-->| GLOBAL NUMBER OF EDGE EXTREMITIES
00120 !| QBOR           |-->| PRESCRIBED BEDLOAD DISCHARGES
00121 !| QSX            |<->| BEDLOAD TRANSPORT RATE X-DIRECTION
00122 !| QSY            |<->| BEDLOAD TRANSPORT RATE Y-DIRECTION
00123 !| T10            |<->| WORK BIEF_OBJ STRUCTURE
00124 !| UNSV2D         |-->| INVERSE OF INTEGRALS OF TEST FUNCTIONS
00125 !| VNOIN          |-->| OUTWARD UNIT NORMALS
00126 !| ZFCL           |<->| BEDLOAD EVOLUTION FOR EACH SEDIMENT CLASS
00127 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00128 !
00129       USE INTERFACE_SISYPHE, EX_BEDLOAD_SOLVS_VF => BEDLOAD_SOLVS_VF
00130       USE BIEF
00131       IMPLICIT NONE
00132       INTEGER LNG,LU
00133       COMMON/INFO/LNG,LU
00134 !
00135 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00136 !
00137       TYPE(BIEF_MESH),  INTENT(INOUT) :: MESH
00138       TYPE(BIEF_OBJ),   INTENT(IN)    :: QSX,QSY,LIMTEC,UNSV2D,EBOR
00139       TYPE(BIEF_OBJ),   INTENT(IN)    :: BREACH,LIQBOR,QBOR
00140       INTEGER,          INTENT(IN)    :: NSEG,NPTFR,NPOIN,KENT,KDIR,KDDL
00141       DOUBLE PRECISION, INTENT(IN)    :: DT,CSF_SABLE
00142       DOUBLE PRECISION, INTENT(IN)    :: AVA(NPOIN)
00143       TYPE(BIEF_OBJ),   INTENT(INOUT) :: T10,FLBCLA,ZFCL,FLUX
00144 ! RA
00145       INTEGER, INTENT(IN)             :: NUBO(2,NSEG)
00146       DOUBLE PRECISION, INTENT(IN)    :: VNOIN(3,NSEG)
00147 !
00148 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00149 !
00150       INTEGER          :: ISEGIN,K,N,IEL,IEL1,IEL2
00151       DOUBLE PRECISION :: QSMOY1,QSMOY2,QSP,VNOIN1,VNOIN2,RNORM,XN,YN
00152       DOUBLE PRECISION :: ZFCLDIR,PROD_SCAL
00153 !
00154 !======================================================================!
00155 !======================================================================!
00156 !                               PROGRAM                                !
00157 !======================================================================!
00158 !======================================================================!
00159 !
00160 !     INTIALISES THE DIVERGENCE
00161 !
00162       CALL OS('X=0     ', X=FLUX)
00163 !
00164 !     DETERMINES THE OUTGOING FLUX FOR EACH CELL
00165 !
00166       DO ISEGIN = 1, NSEG
00167 !
00168         IEL1 = NUBO(1,ISEGIN)
00169         IEL2 = NUBO(2,ISEGIN)
00170 !
00171         ! II.1 - SEGMENT LENGTH (RNORM)
00172         ! ----------------------------------
00173         VNOIN1 = VNOIN(1,ISEGIN)
00174         VNOIN2 = VNOIN(2,ISEGIN)
00175         RNORM  = VNOIN(3,ISEGIN)
00176         PROD_SCAL= (MESH%X%R(IEL2)-MESH%X%R(IEL1))*VNOIN1+
00177      &             (MESH%Y%R(IEL2)-MESH%Y%R(IEL1))*VNOIN2
00178         IF(PROD_SCAL.LT.0.D0)THEN
00179           IEL1 = NUBO(2,ISEGIN)
00180           IEL2 = NUBO(1,ISEGIN)
00181         ENDIF
00182         ! II.2 - QS FOR THE SEGMENT, BROKEN UP ACCORDING TO X AND Y
00183         ! ---------------------------------------------
00184         QSMOY1 = 0.5D0*(QSX%R(IEL1) + QSX%R(IEL2))
00185         QSMOY2 = 0.5D0*(QSY%R(IEL1) + QSY%R(IEL2))
00186         ! II.3 - PROJECTS QS FOR THE SEGMENT ONTO THE SEGMENT NORMAL
00187         ! ------------------------------------------------------------
00188         QSP = VNOIN1*QSMOY1 + VNOIN2*QSMOY2
00189         ! II.4 - UPWIND SCHEME ON NODES WITH A "PROBLEM"
00190         ! ----------------------------------------------
00191         IF(BREACH%I(IEL1).EQ.1.AND.QSP.GT.0.D0) THEN
00192           QSMOY1 = QSX%R(IEL1)
00193           QSMOY2 = QSY%R(IEL1)
00194         ENDIF
00195         IF(BREACH%I(IEL2).EQ.1.AND.QSP.LT.0.D0) THEN
00196           QSMOY1 = QSX%R(IEL2)
00197           QSMOY2 = QSY%R(IEL2)
00198         ENDIF
00199         QSP = VNOIN1*QSMOY1 + VNOIN2*QSMOY2
00200         ! II.5 - INTEGRATES BY THE SEGMENT LENGTH
00201         ! ---------------------------------------------
00202         FLUX%R(IEL1) = FLUX%R(IEL1) + RNORM*QSP
00203         FLUX%R(IEL2) = FLUX%R(IEL2) - RNORM*QSP
00204       ENDDO
00205 !
00206 !     BOUNDARIES
00207 !
00208       DO K = 1 , NPTFR
00209         IEL = MESH%NBOR%I(K)
00210 !       VECTOR NORMAL TO A BOUNDARY NODE
00211 !       VERSION WHICH IS NOT NORMED
00212         XN = MESH%XNEBOR%R(K+NPTFR)
00213         YN = MESH%YNEBOR%R(K+NPTFR)
00214 !
00215 !       ADDING BOUNDARY FLUX AT OPEN BOUNDARIES
00216 !       QBOR HAS PRIORITY HERE, SO IN CASE OF LIQBOR=KENT
00217 !       LIMTEC=KDIR WILL NOT BE CONSIDERED, SEE BEDLOAD_SOLVS_FE
00218 !       FOR MORE EXPLANATIONS
00219 !
00220         IF(LIQBOR%I(K).EQ.KENT) THEN
00221           FLBCLA%R(K) = QBOR%R(K)
00222           FLUX%R(IEL) = FLUX%R(IEL) + FLBCLA%R(K)
00223         ELSEIF(LIMTEC%I(K).EQ.KDDL.OR.LIMTEC%I(K).EQ.KDIR) THEN
00224 !         IF KDIR WILL BE UPDATED LATER
00225           FLBCLA%R(K)= QSX%R(IEL)*XN + QSY%R(IEL)*YN
00226 !         ADDS THE CONTRIBUTION OF THE FLUX ON THE BOUNDARY SEGMENT
00227           FLUX%R(IEL) = FLUX%R(IEL) + FLBCLA%R(K)
00228         ELSE
00229 !         NO SEDIMENT FLUX ACCROSS SOLID BOUNDARIES
00230           FLBCLA%R(K)=0.D0
00231         ENDIF
00232       ENDDO
00233 !
00234 !     ASSEMBLING IN PARALLEL
00235 !
00236       IF(NCSIZE.GT.1) CALL PARCOM(FLUX, 2, MESH)
00237 !
00238 !     SOLVING, NEGATIVE SIGN BECAUSE OUTGOING FLUX IS POSITIVE
00239 !     NOTE JMH: FLUX MUST BE HERE DIV(QS)
00240 !
00241       CALL OS('X=CYZ   ', X=ZFCL, Y=FLUX, Z=UNSV2D, C=-DT)
00242 !
00243       DO K=1,NPTFR
00244 !                                  PRIORITY OF LIQBOR, SEE ABOVE
00245         IF(LIMTEC%I(K).EQ.KDIR.AND.LIQBOR%I(K).NE.KENT) THEN
00246           N = MESH%NBOR%I(K)
00247 !         ZFCLDIR: DIRICHLET VALUE OF EVOLUTION, ZFCL WILL BE DIVIDED BY
00248 !         CSF_SABLE AFTER, AND THEN IT WILL BE AVA(N)*EBOR...
00249           ZFCLDIR = AVA(N)*EBOR%R(K)*CSF_SABLE
00250 !         CORRECTION OF BOUNDARY FLUX TO GET ZFCLDIR
00251           FLBCLA%R(K)=FLBCLA%R(K)-(ZFCLDIR-ZFCL%R(N))/(DT*UNSV2D%R(N))
00252 !         ZFCLDIR FINALLY IMPOSED
00253           ZFCL%R(N) = ZFCLDIR
00254         ENDIF
00255 !
00256       ENDDO
00257 !
00258 !======================================================================!
00259 !======================================================================!
00260 !
00261       RETURN
00262       END

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