maxslope.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\sisyphe\maxslope.f
00002 !
00071                      SUBROUTINE MAXSLOPE
00072 !                    *******************
00073 !
00074      &(SLOPE,ZF,ZR,XEL,YEL,NELEM,NELMAX,NPOIN,IKLE,EVOL,UNSV2D,MESH,
00075      & ZFCL_MS,AVAIL,NOMBLAY,NSICLA)
00076 !
00077 !***********************************************************************
00078 ! SISYPHE   V6P3                                   21/07/2011
00079 !***********************************************************************
00080 !
00081 !
00082 !
00083 !
00084 !
00085 !
00086 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00087 !| EVOL           |<->| WORK ARRAY, THEN EVOLUTION DUE TO SLIDE
00088 !| IKLE           |-->| CONNECTIVITY TABLE
00089 !| MESH           |-->| MESH STRUCTURE
00090 !| NELEM          |-->| NUMBER OF ELEMENTS
00091 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00092 !| NPOIN          |-->| NUMBER OF POINTS IN THE MESH
00093 !| SLOPE          |-->| MAXIMUM SLOPE IN DEGREES
00094 !| UNSV2D         |-->| INVERSE OF INTEGRAL OF BASES
00095 !| XEL,YEL        |-->| MESH COORDINATES PER ELEMENT
00096 !| ZF             |<->| BOTTOM THAT WILL BE MODIFIED
00097 !| ZR             |-->| NON ERODABLE BED
00098 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00099 !
00100       USE BIEF
00101 !
00102       USE INTERFACE_SISYPHE, EX_MAXSLOPE => MAXSLOPE
00103 !
00104       IMPLICIT NONE
00105       INTEGER LNG,LU
00106       COMMON/INFO/LNG,LU
00107 !
00108 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00109 !
00110       INTEGER, INTENT(IN) :: NELEM,NELMAX,NPOIN,NOMBLAY,NSICLA
00111       INTEGER, INTENT(IN) :: IKLE(NELMAX,3)
00112 !
00113       DOUBLE PRECISION, INTENT(IN   ) :: SLOPE
00114       DOUBLE PRECISION, INTENT(INOUT) :: ZF(NPOIN)
00115       DOUBLE PRECISION, INTENT(IN)    :: ZR(NPOIN)
00116       DOUBLE PRECISION, INTENT(IN)    :: XEL(NELMAX,3),YEL(NELMAX,3)
00117       DOUBLE PRECISION, INTENT(INOUT) :: AVAIL(NPOIN,NOMBLAY,NSICLA)
00118 !
00119       TYPE(BIEF_OBJ), INTENT(INOUT)   :: EVOL,ZFCL_MS
00120       TYPE(BIEF_OBJ), INTENT(IN)      :: UNSV2D
00121       TYPE(BIEF_MESH) :: MESH
00122 !
00123 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00124 !
00125       INTEGER IELEM,I1,I2,I3,I,IG1,IG2,IR1,IR2,J
00126       DOUBLE PRECISION X2,X3,Y2,Y3,Z2,Z3,A,B,L,ZC,DEUXSURF,TANSL
00127       DOUBLE PRECISION Q(3),QG1,QG2,QR1,QR2
00128 !
00129       LOGICAL CASE2
00130 !
00131       INTRINSIC SQRT,MIN,MAX,TAN
00132 !
00133 !-----------------------------------------------------------------------
00134 !
00135       TANSL=TAN(3.141592653589D0*SLOPE/180.D0)
00136 !
00137 !     INITIALISES THE RIGHT-HAND SIDE EVOL TO ZERO
00138 !
00139       CALL CPSTVC(UNSV2D,EVOL)
00140       CALL OS('X=0     ',X=EVOL)
00141 !
00142 !     ONE CLASS VERSION
00143 !
00144       IF(NSICLA.EQ.1) THEN
00145 !
00146 !       LOOP ON ELEMENTS
00147 !
00148         DO IELEM=1,NELEM
00149 !
00150           I1=IKLE(IELEM,1)
00151           I2=IKLE(IELEM,2)
00152           I3=IKLE(IELEM,3)
00153 !
00154           X2=XEL(IELEM,2)
00155           X3=XEL(IELEM,3)
00156           Y2=YEL(IELEM,2)
00157           Y3=YEL(IELEM,3)
00158           Z2=ZF(I2)-ZF(I1)
00159           Z3=ZF(I3)-ZF(I1)
00160 !
00161 !         TWICE THE TRIANGLE AREA
00162 !
00163           DEUXSURF=X2*Y3-X3*Y2
00164 !
00165 !         AVERAGE BOTTOM IN THE ELEMENT
00166 !
00167           ZC=(ZF(I1)+ZF(I2)+ZF(I3))/3.D0
00168 !
00169 !         COMPONENTS OF BOTTOM GRADIENT
00170 !
00171           A=(Z2*Y3-Z3*Y2)/DEUXSURF
00172           B=(Z3*X2-Z2*X3)/DEUXSURF
00173 !
00174 !         CORRECTING FACTOR ON SLOPE
00175 !
00176           L=MIN(1.D0,TANSL/MAX(SQRT(A**2+B**2),1.D-8))
00177 !
00178 !         L LIMITED DUE TO NON-ERODABLE BEDS : ZF MUST NOT GO BELOW ZR
00179 !
00180           IF(ZF(I1).GT.ZC) L=MAX(L,(ZR(I1)-ZC)/MAX(ZF(I1)-ZC,1.D-8))
00181           IF(ZF(I2).GT.ZC) L=MAX(L,(ZR(I2)-ZC)/MAX(ZF(I2)-ZC,1.D-8))
00182           IF(ZF(I3).GT.ZC) L=MAX(L,(ZR(I3)-ZC)/MAX(ZF(I3)-ZC,1.D-8))
00183 !
00184 !         BUILDS THE RIGHT-HAND SIDE
00185 !
00186 !         HERE THE EVOLUTIONS ARE MULTIPLIED BY SURFAC/3
00187 !         BECAUSE THE REAL EVOLUTION TAKING INTO ACCOUNT OTHER ELEMENTS
00188 !         WILL NEED A FACTOR (SURFAC/3)/(INTEGRAL OF BASIS)
00189 !
00190           EVOL%R(I1)=EVOL%R(I1)+(1.D0-L)*(ZC-ZF(I1))*DEUXSURF/6.D0
00191           EVOL%R(I2)=EVOL%R(I2)+(1.D0-L)*(ZC-ZF(I2))*DEUXSURF/6.D0
00192           EVOL%R(I3)=EVOL%R(I3)+(1.D0-L)*(ZC-ZF(I3))*DEUXSURF/6.D0
00193 !
00194         ENDDO
00195 !
00196       ELSE
00197 !
00198 !       MULTI-CLASS VERSION
00199 !
00200 !       INITIALING TO 0. THE EVOLUTIONS DUE TO SLIDE FOR EACH CLASS
00201 !
00202         DO I=1,NSICLA
00203           CALL OS('X=0     ',X=ZFCL_MS%ADR(I)%P)
00204         ENDDO
00205 !
00206 !       LOOP ON ELEMENTS
00207 !
00208         DO IELEM=1,NELEM
00209 !
00210           I1=IKLE(IELEM,1)
00211           I2=IKLE(IELEM,2)
00212           I3=IKLE(IELEM,3)
00213 !
00214           X2=XEL(IELEM,2)
00215           X3=XEL(IELEM,3)
00216           Y2=YEL(IELEM,2)
00217           Y3=YEL(IELEM,3)
00218           Z2=ZF(I2)-ZF(I1)
00219           Z3=ZF(I3)-ZF(I1)
00220 !
00221 !         TWICE THE TRIANGLE AREA
00222 !
00223           DEUXSURF=X2*Y3-X3*Y2
00224 !
00225 !         AVERAGE BOTTOM IN THE ELEMENT
00226 !
00227           ZC=(ZF(I1)+ZF(I2)+ZF(I3))/3.D0
00228 !
00229 !         COMPONENTS OF BOTTOM GRADIENT
00230 !
00231           A=(Z2*Y3-Z3*Y2)/DEUXSURF
00232           B=(Z3*X2-Z2*X3)/DEUXSURF
00233 !
00234 !         CORRECTING FACTOR ON SLOPE
00235 !
00236           L=MIN(1.D0,TANSL/MAX(SQRT(A**2+B**2),1.D-8))
00237 !
00238 !         L LIMITED DUE TO NON-ERODABLE BEDS : ZF MUST NOT GO BELOW ZR
00239 !
00240           IF(ZF(I1).GT.ZC) L=MAX(L,(ZR(I1)-ZC)/MAX(ZF(I1)-ZC,1.D-8))
00241           IF(ZF(I2).GT.ZC) L=MAX(L,(ZR(I2)-ZC)/MAX(ZF(I2)-ZC,1.D-8))
00242           IF(ZF(I3).GT.ZC) L=MAX(L,(ZR(I3)-ZC)/MAX(ZF(I3)-ZC,1.D-8))
00243 !
00244 !         BUILDS THE RIGHT-HAND SIDE
00245 !
00246 !         HERE THE EVOLUTIONS ARE MULTIPLIED BY SURFAC/3
00247 !         BECAUSE THE REAL EVOLUTION TAKING INTO ACCOUNT OTHER ELEMENTS
00248 !         WILL NEED A FACTOR (SURFAC/3)/(INTEGRAL OF BASIS)
00249 !
00250 !         FIRST IN TERMS OF QUANTITIES BROUGHT TO POINTS
00251 !
00252           Q(1)=(1.D0-L)*(ZC-ZF(I1))*DEUXSURF/6.D0
00253           Q(2)=(1.D0-L)*(ZC-ZF(I2))*DEUXSURF/6.D0
00254           Q(3)=(1.D0-L)*(ZC-ZF(I3))*DEUXSURF/6.D0
00255 !
00256           EVOL%R(I1)=EVOL%R(I1)+Q(1)
00257           EVOL%R(I2)=EVOL%R(I2)+Q(2)
00258           EVOL%R(I3)=EVOL%R(I3)+Q(3)
00259 !
00260 !         TAKING INTO ACCOUNT THE QUANTITIES TO UPDATE ZFCL_MS
00261 !         IG1 AND IG2 : POINTS THAT GIVE
00262 !         IR1 AND IR2 : POINTS THAT RECEIVE
00263 !         CASE2: TWO POINTS GIVE TO THE THIRD ONE (THE OTHER CASE IS
00264 !                ONE POINT GIVES TO THE TWO OTHERS)
00265           CASE2=.FALSE.
00266 !
00267 !         PARAMETERISING TO REDUCE THE 6 CASES TO 2
00268 !
00269           IF(Q(1).GE.0.D0) THEN
00270             IF(Q(2).GE.0.D0) THEN
00271 !             3 GIVES TO 1 AND 2
00272               IG1=I3
00273               QG1=Q(3)
00274               IR1=I1
00275               QR1=Q(1)
00276               IR2=I2
00277               QR2=Q(2)
00278             ELSE
00279               IF(Q(3).GE.0.D0) THEN
00280 !               2 GIVES TO 1 AND 3
00281                 IG1=I2
00282                 QG1=Q(2)
00283                 IR1=I1
00284                 QR1=Q(1)
00285                 IR2=I3
00286                 QR2=Q(3)
00287               ELSE
00288 !               2 AND 3 GIVE TO 1
00289                 IG1=I2
00290                 QG1=Q(2)
00291                 IG2=I3
00292                 QG2=Q(3)
00293                 IR1=I1
00294                 QR1=Q(1)
00295                 CASE2=.TRUE.
00296               ENDIF
00297             ENDIF
00298           ELSE
00299             IF(Q(2).GT.0.D0) THEN
00300               IF(Q(3).GT.0.D0) THEN
00301 !               1 GIVES TO 2 AND 3
00302                 IG1=I1
00303                 QG1=Q(1)
00304                 IR1=I2
00305                 QR1=Q(2)
00306                 IR2=I3
00307                 QR2=Q(3)
00308               ELSE
00309 !               1 AND 3 GIVE TO 2
00310                 IG1=I1
00311                 QG1=Q(1)
00312                 IG2=I3
00313                 QG2=Q(3)
00314                 IR1=I2
00315                 QR1=Q(2)
00316                 CASE2=.TRUE.
00317               ENDIF
00318             ELSE
00319 !             1 AND 2 GIVE TO 3
00320               IG1=I1
00321               QG1=Q(1)
00322               IG2=I2
00323               QG2=Q(2)
00324               IR1=I3
00325               QR1=Q(3)
00326               CASE2=.TRUE.
00327             ENDIF
00328           ENDIF
00329 !
00330           IF(CASE2) THEN
00331 !
00332 !           THE TWO DONNORS CASE : IG1 AND IG2 GIVE TO IR1
00333 !           ZFCL_MS IS HERE VOLUMES
00334 !
00335             DO I=1,NSICLA
00336               ZFCL_MS%ADR(I)%P%R(IG1)=ZFCL_MS%ADR(I)%P%R(IG1)
00337      &                               +QG1*AVAIL(IG1,1,I)
00338               ZFCL_MS%ADR(I)%P%R(IG2)=ZFCL_MS%ADR(I)%P%R(IG2)
00339      &                               +QG2*AVAIL(IG2,1,I)
00340               ZFCL_MS%ADR(I)%P%R(IR1)=ZFCL_MS%ADR(I)%P%R(IR1)
00341      &                               -QG1*AVAIL(IG1,1,I)
00342      &                               -QG2*AVAIL(IG2,1,I)
00343             ENDDO
00344 !
00345           ELSE
00346 !
00347 !           THE ONE DONNOR CASE : IG1 GIVES TO IR1 AND IR2
00348 !           ZFCL_MS IS HERE VOLUMES
00349 !
00350             DO I=1,NSICLA
00351               ZFCL_MS%ADR(I)%P%R(IG1)=ZFCL_MS%ADR(I)%P%R(IG1)
00352      &                               +QG1*AVAIL(IG1,1,I)
00353               ZFCL_MS%ADR(I)%P%R(IR1)=ZFCL_MS%ADR(I)%P%R(IR1)
00354      &                               +QR1*AVAIL(IG1,1,I)
00355               ZFCL_MS%ADR(I)%P%R(IR2)=ZFCL_MS%ADR(I)%P%R(IR2)
00356      &                               +QR2*AVAIL(IG1,1,I)
00357             ENDDO
00358 !
00359           ENDIF
00360 !
00361         ENDDO
00362 !
00363 !       ADDING VOLUMES IN PARALLEL
00364 !
00365         IF(NCSIZE.GT.1) THEN
00366           DO I=1,NSICLA
00367             CALL PARCOM(ZFCL_MS%ADR(I)%P,2,MESH)
00368           ENDDO
00369         ENDIF
00370 !
00371 !       FINAL DIVISION BY THE INTEGRAL OF BASES: VOLUMES CHANGED INTO
00372 !       BED VARIATIONS (LIKE EVOL BELOW)
00373 !
00374         DO I=1,NSICLA
00375           DO J=1,NPOIN
00376             ZFCL_MS%ADR(I)%P%R(J)=ZFCL_MS%ADR(I)%P%R(J)*UNSV2D%R(J)
00377           ENDDO
00378         ENDDO
00379 !
00380       ENDIF
00381 !
00382 !-----------------------------------------------------------------------
00383 !
00384 !     FINAL RESOLUTION
00385 !
00386       IF(NCSIZE.GT.1) THEN
00387         CALL PARCOM(EVOL,2,MESH)
00388       ENDIF
00389 !
00390 !     FINAL DIVISION BY THE INTEGRAL OF BASES: QUANTITIES CHANGED INTO
00391 !     ELEVATIONS
00392 !
00393       DO I=1,NPOIN
00394         EVOL%R(I)=EVOL%R(I)*UNSV2D%R(I)
00395       ENDDO
00396 !
00397 !-----------------------------------------------------------------------
00398 !
00399       RETURN
00400       END

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