layer.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\sisyphe\layer.f
00002 !
00094                      SUBROUTINE LAYER
00095 !                    ****************
00096 !
00097      &(ZFCL_W,NLAYER,ZR,ZF,ESTRAT,ELAY,MASBAS,ACLADM,NSICLA,NPOIN,
00098      & ELAY0,VOLTOT,ES,AVAIL,CONST_ALAYER,DTS,ESTRATNEW,NLAYNEW)
00099 !
00100 !***********************************************************************
00101 ! SISYPHE   V6P2                                   21/07/2011
00102 !***********************************************************************
00103 !
00104 !
00105 !
00106 !
00107 !
00108 !
00109 !
00110 !
00111 !
00112 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00113 !| ACLADM         |-->| MEAN DIAMETER OF SEDIMENT
00114 !| AVAIL          |<->| SEDIMENT FRACTION FOR EACH LAYER, CLASS, POINT
00115 !| CONST_ALAYER   |-->| CONSTANT ACTIVE LAYER THICKNESS OR NOT
00116 !| DTS            |-->| TIME STEP FOR SUSPENSION
00117 !| ELAY           |<->| ACTIVE LAYER THICKNESS FOR EACH POINT
00118 !| ELAY0          |<->| ACTIVE LAYER THICKNESS
00119 !| ES             |<->| LAYER THICKNESSES AS DOUBLE PRECISION
00120 !| ESTRAT         |<->| ACTIVE STRATUM THICKNESS FOR EACH POINT
00121 !| ESTRATNEW      |<->| ACTIVE STRATUM THICKNESS AT TIME T+DT
00122 !| MASBAS         |-->| INTEGRAL OF TEST FUNCTIONS
00123 !| NLAYER         |<--| NUMBER OF LAYER FOR EACH POINT
00124 !| NLAYNEW        |<->| NUMBER OF LAYER AT TIME T+DT
00125 !| NPOIN          |-->| NUMBER OF POINTS
00126 !| NSICLA         |-->| NUMBER OF SIZE CLASSES FOR BED MATERIALS
00127 !| VOLTOT         |<->| TOTAL VOLUME OF SEDIMENT IN THE BED
00128 !| ZF             |-->| ELEVATION OF BOTTOM
00129 !| ZFCL_W         |-->| BED EVOLUTION FOR EACH SEDIMENT CLASS
00130 !| ZR             |-->| NON ERODABLE BED
00131 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00132 !
00133       USE BIEF
00134       USE DECLARATIONS_SISYPHE , ONLY : NOMBLAY
00135       IMPLICIT NONE
00136       INTEGER LNG,LU
00137       COMMON/INFO/LNG,LU
00138 !
00139 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00140 !
00141       TYPE (BIEF_OBJ),  INTENT(IN)    :: ZFCL_W,ZR,ZF
00142       TYPE (BIEF_OBJ),  INTENT(IN)    :: MASBAS,ACLADM
00143       INTEGER,          INTENT(IN)    :: NSICLA,NPOIN
00144       DOUBLE PRECISION, INTENT(IN)    :: DTS
00145       LOGICAL,          INTENT(IN)    :: CONST_ALAYER
00146       TYPE (BIEF_OBJ),  INTENT(INOUT) :: NLAYER,ESTRAT,ELAY
00147       DOUBLE PRECISION, INTENT(INOUT) :: ELAY0
00148       DOUBLE PRECISION, INTENT(INOUT) :: ES(NPOIN,NOMBLAY)
00149       DOUBLE PRECISION, INTENT(INOUT) :: AVAIL(NPOIN,NOMBLAY,NSICLA)
00150       DOUBLE PRECISION, INTENT(INOUT) :: VOLTOT(NSICLA),ESTRATNEW(NPOIN)
00151       INTEGER         , INTENT(INOUT) :: NLAYNEW(NPOIN)
00152 !
00153 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00154 !
00155       DOUBLE PRECISION P_DSUM
00156       EXTERNAL         P_DSUM
00157       INTEGER  P_ISUM
00158       EXTERNAL P_ISUM
00159 !
00160 !-----------------------------------------------------------------------
00161 !
00162       INTEGER I,J,K,ARRET,ARRET2
00163       DOUBLE PRECISION EVOL,HEIGH,TEST1,TEST2,AUX
00164 !
00165 !-----------------------------------------------------------------------
00166 !
00167 !     TO CHECK FRACTIONS IN THE RANGE [-ZERO,1+ZERO]
00168 !
00169       DOUBLE PRECISION ZERO
00170 !     DATA             ZERO/1.D-10/
00171 !     IN CASE OF RESTART, THE FIRST TIME STEP IS A BIT HARD BECAUSE OF
00172 !     SINGLE PRECISION, WITHOUT RESTART 1.D-10 IS OK
00173       DATA             ZERO/1.D-7/
00174 !
00175 !-----------------------------------------------------------------------
00176 !
00177       ARRET=0
00178 !
00179 !-----------------------------------------------------------------------
00180 !
00181       DO J=1,NPOIN
00182 !
00183 !       ACTIVATE BEFORE INVESTIGATING PROBLEMS IN LAYER...
00184 !       LOOK FOR 'ACTIVATE' TO SEE OTHER LINES LIKE THIS BELOW
00185 !
00186 !       IF(ELAY%R(J).LT.0.D0) THEN
00187 !         WRITE(LU,*) 'NEGATIVE ELAY IN LAYER J=',J,' ELAY=',ELAY%R(J)
00188 !         CALL PLANTE(1)
00189 !         STOP
00190 !       ENDIF
00191 !
00192         IF(.NOT.CONST_ALAYER) ELAY0 = 3.D0 * ACLADM%R(J)
00193 !
00194         NLAYNEW(J) = NLAYER%I(J)
00195 !
00196 ! QUESTION JMH, EVOLUTION HAS BEEN COMPUTED BEFORE IN ARRAY E, WHY NOT
00197 !               EVOL=E(J) ?????
00198 !               ELAY(J) = ES(J,1) WHY IS IT AN EXTRA ARRAY ??
00199 !
00200 !
00201         HEIGH = ZF%R(J)-ZR%R(J)
00202 !
00203 !       ACTIVATE BEFORE INVESTIGATING PROBLEMS IN LAYER...
00204 !
00205 !       IF(HEIGH.LT.0.D0) THEN
00206 !         WRITE(LU,*) 'BAD DATA IN LAYER J=',J,' HEIGH=',HEIGH
00207 !         CALL PLANTE(1)
00208 !         STOP
00209 !       ENDIF
00210 !
00211 !       HERE ELAY.NE.HEIGH BECAUSE ELAY IS THE ACTIVE LAYER THICKNESS
00212         EVOL  = 0.D0
00213         DO I=1,NSICLA
00214           EVOL = EVOL + ZFCL_W%ADR(I)%P%R(J)
00215         ENDDO
00216 !
00217         IF(NLAYER%I(J).GT.1) THEN
00218 !
00219           IF(EVOL.GE.0.D0) THEN
00220 !
00221 !           DEPOSITION
00222 !
00223 !           NEW HEIGHT OF LAYER 2 (IT RECEIVES EVOL TO KEEP LAYER 1 CONSTANT)
00224             ESTRATNEW(J) = ESTRAT%R(J) + EVOL
00225 !
00226             DO I=1,NSICLA
00227 !             JMH 28/04/2010. THE OLD IMPLEMENTATION CONSISTED OF FIRST
00228 !             GIVING EVOL TO LAYER 2, WITH OLD AVAIL, THEN OF RECEIVING
00229 !             THE DEPOSITION, BUT IF A CLASS DISAPPEARS IN LAYER 1,
00230 !             IT IS NOT POSSIBLE TO GIVE IT FIRST TO LAYER 2, SO NEW
00231 !             FRACTIONS MUST BE COMPUTED BEFORE GIVING EVOL TO LAYER 2
00232 !             THEN I SEE NO DIFFERENCE BETWEEN EVOLELAY0
00233 !             SO BOTH ARE TREATED BELOW, UNLIKE RELEASE 5.9.
00234 !
00235 !             1) LAYER 1 RECEIVES ZFCL_W OF CLASS I, WE COMPUTE THE
00236 !                PROVISIONAL NEW AVAIL(J,1,I) IN AUX
00237               AUX=(AVAIL(J,1,I)*ELAY0+ZFCL_W%ADR(I)%P%R(J))/
00238      &                        (ELAY0+EVOL)
00239 !
00240 !             2) LAYER 2 RECEIVES AUX*EVOL OF CLASS I (AUX MAY BE 0 HERE)
00241               AVAIL(J,2,I)=(AUX*EVOL+AVAIL(J,2,I)*ESTRAT%R(J))/
00242      &                               ESTRATNEW(J)
00243 !
00244 !             3) SEEN FROM LAYER 1, AUX*EVOL OF CLASS I HAS BEEN GIVEN
00245 !                AND THE NEW LAYER THICKNESS IS ELAY0, HENCE THE NEW AVAIL
00246               AVAIL(J,1,I)=( AVAIL(J,1,I)*ELAY0-AUX*EVOL
00247      &                          +ZFCL_W%ADR(I)%P%R(J) )/ELAY0
00248 !
00249 ! NOTE CV: CAN BE REPLACED BY
00250 !              AVAIL(J,1,I)= AUX
00251 !
00252 !             OLD (AND WRONG) FORMULATION (IN IT -AVAIL*EVOL SHOULD BE -AUX*EVOL)
00253 !             AVAIL(J,1,I)=( AVAIL(J,1,I)*(ELAY0-EVOL)
00254 !    &                       +ZFCL_W%ADR(I)%P%R(J)     )/ELAY0
00255 !
00256               IF(AVAIL(J,1,I).GT.1.D0+ZERO.OR.
00257      &           AVAIL(J,1,I).LT.-ZERO) THEN
00258                 WRITE(LU,*) 'ERROR IN LAYER CASE 1'
00259                 CALL PLANTE(1)
00260                 STOP
00261               ENDIF
00262             ENDDO
00263 !           NEW HEIGHT OF LAYER 1
00264             ELAY%R(J) = ELAY0
00265 !
00266           ELSEIF(EVOL.GT.-ELAY0) THEN
00267 ! CV: I DON'T AGREE WITH THE COMMENT BELOW
00268 !     WE'RE IN THE CASE : -ELAY0<EVOL<0 HENCE ELAY0>-EVOL>0
00269 !     THICKNESS OF THE FIRST LAYER IS THEREFORE SUFFICIENT
00270 !
00271 !           EROSION GREATER THAN LAYER 1, WE HAVE TO DESTROY A STRATUM
00272 !
00273             IF(-EVOL.GT.ESTRAT%R(J)) THEN
00274 !
00275 !  CV: USUALLY, LAYER 2 IS VERY BROAD AND TWO LAYERS ARE IN GENERAL SUFFICIENT
00276 !      HERE LAYER 2 IS DESTROYED
00277 !
00278 !             USUAL CASE (NOTE JMH : WHY NOT .GE.2 ?
00279 !
00280               IF(NLAYER%I(J).GT.2) THEN
00281 !
00282                 DO I=1,NSICLA
00283                   AVAIL(J,1,I) = (  AVAIL(J,1,I)*ELAY0
00284      &                            + ZFCL_W%ADR(I)%P%R(J)
00285      &                            + AVAIL(J,2,I)*ESTRAT%R(J)
00286      &                            - AVAIL(J,3,I)*(EVOL+ESTRAT%R(J))
00287      &                           )/ ELAY0
00288                   IF(AVAIL(J,1,I).GT.1.D0+ZERO.OR.
00289      &               AVAIL(J,1,I).LT.-ZERO) THEN
00290                     WRITE(LU,*) 'ERROR IN LAYER CASE 2'
00291                     CALL PLANTE(1)
00292                     STOP
00293                   ENDIF
00294                   AVAIL(J,2,I) = AVAIL(J,3,I)
00295 ! CV : change 9 in NOMBLAY-1 ??
00296 ! Pourquoi pas NOMBLAY?                          ?=
00297 !                  DO K=3,MIN(9,NLAYER%I(J))
00298                  DO K=3,MIN(NOMBLAY-1,NLAYER%I(J))
00299                     AVAIL(J,K,I) = AVAIL(J,K+1,I)
00300                   ENDDO
00301                 ENDDO
00302                 ELAY%R(J) = ELAY0
00303                 NLAYNEW(J) = NLAYER%I(J) - 1
00304                 ESTRATNEW(J) = ESTRAT%R(J) + EVOL + ES(J,3)
00305 !CV              DO K=3,MIN(9,NLAYER%I(J))
00306 ! Pourquoi pas NOMBLAY?
00307                 DO K=3,MIN(NOMBLAY-1,NLAYER%I(J))
00308                   ES(J,K) = ES(J,K+1)
00309                 ENDDO
00310 !
00311 !             HERE NLAYER.GT.1 AND NOT .GT.2, SO 2 !
00312 !
00313               ELSE
00314                 DO I=1,NSICLA
00315                   IF(HEIGH.GT.0.D0) THEN
00316                     AVAIL(J,1,I) = (  AVAIL(J,1,I)*ELAY0
00317      &                              + ZFCL_W%ADR(I)%P%R(J)
00318      &                              + ESTRAT%R(J)*AVAIL(J,2,I)
00319      &                              )/(ELAY0+EVOL+ESTRAT%R(J))
00320 !                   MODIF JMH 12/04/2011, WITH 2 LAYERS, LAYER 1
00321 !                   HAS ELAY%R(J)=ELAY0, SO WHY ELAY IN DENOMINATOR
00322 !                   AND ELAY0 IN NUMERATOR ? (VICIOUS!!)
00323 !    &                              )/(ELAY%R(J)+EVOL+ESTRAT%R(J))
00324                     IF(AVAIL(J,1,I).GT.1.D0+ZERO.OR.
00325      &                 AVAIL(J,1,I).LT.-ZERO) THEN
00326                       WRITE(LU,*) 'J=',J,' NLAYER%I(J)=',NLAYER%I(J)
00327                       WRITE(LU,*) 'AVAIL(J,1,I)=',AVAIL(J,1,I)
00328                       WRITE(LU,*) 'AVAIL(J,2,I)=',AVAIL(J,2,I)
00329                       WRITE(LU,*) 'ZFCL=',ZFCL_W%ADR(I)%P%R(J)
00330                       WRITE(LU,*) 'HEIGH=',HEIGH,' ELAY0=',ELAY0
00331                       WRITE(LU,*) 'ESTRAT%R(J)=',ESTRAT%R(J)
00332                       WRITE(LU,*) 'ELAY%R(J)=',ELAY%R(J)
00333                       WRITE(LU,*) 'EVOL=',EVOL
00334                       WRITE(LU,*) 'ERROR IN LAYER CASE 3'
00335                       CALL PLANTE(1)
00336                       STOP
00337                     ENDIF
00338                   ELSE
00339                     AVAIL(J,1,I) = 0.D0
00340                   ENDIF
00341                   AVAIL(J,2,I) = 0.D0
00342                 ENDDO
00343                 NLAYNEW(J) = NLAYER%I(J) - 1
00344                 ELAY%R(J) = HEIGH
00345                 ESTRATNEW(J) = 0.D0
00346               ENDIF
00347 !
00348 !           ONLY LAYER 1 ERODED
00349 !
00350             ELSE
00351               DO I=1,NSICLA
00352                 AVAIL(J,1,I) = (  AVAIL(J,1,I) * ELAY0
00353      &                          + ZFCL_W%ADR(I)%P%R(J)
00354      &                          - EVOL*AVAIL(J,2,I)    )/ELAY0
00355                 IF(AVAIL(J,1,I).GT.1.D0+ZERO.OR.
00356      &             AVAIL(J,1,I).LT.-ZERO) THEN
00357                   WRITE(LU,*) 'ERROR IN LAYER CASE 4'
00358                   CALL PLANTE(1)
00359                   STOP
00360                 ENDIF
00361               ENDDO
00362               ELAY%R(J) = ELAY0
00363               ESTRATNEW(J) = ESTRAT%R(J) + EVOL
00364             ENDIF
00365 !
00366           ELSE
00367 !
00368 !           STOPS IF EROSION IS GREATER THAN
00369 !           THICKNESS OF THE ACTIVE LAYER!
00370 !
00371             IF(LNG.EQ.1) THEN
00372               WRITE(LU,*) 'EROSION TROP FORTE AU NOEUD J=',J
00373               WRITE(LU,*) 'DIMINUER DT OU AUGMENTER ELAY0'
00374             ENDIF
00375             IF(LNG.EQ.2) THEN
00376               WRITE(LU,*) 'TOO MUCH EROSION AT POINT J=',J
00377               WRITE(LU,*) 'DECREASE DT OR INCREASE ELAY0'
00378             ENDIF
00379             WRITE(LU,*) 'EVOL=', EVOL, 'ELAY0=',ELAY0
00380             CALL PLANTE(1)
00381             STOP
00382 !
00383 !         END OF TESTS ON EVOL
00384 !
00385           ENDIF
00386 !
00387 ! THERE WAS ONLY ONE LAYER
00388 ! ------------------------
00389 !
00390         ELSE
00391 !
00392 !         IT IS NOW BIG ENOUGH TO MAKE TWO LAYERS
00393 !
00394           IF(HEIGH.GT.ELAY0) THEN
00395             NLAYNEW(J) = 2
00396             ESTRATNEW(J) = HEIGH - ELAY0
00397             ELAY%R(J) = ELAY0
00398             DO I=1,NSICLA
00399               AVAIL(J,2,I) = AVAIL(J,1,I)
00400               AVAIL(J,1,I) = (AVAIL(J,1,I) * (ELAY0-EVOL)
00401      &                     + ZFCL_W%ADR(I)%P%R(J) )/ELAY0
00402               IF(AVAIL(J,1,I).GT.1.D0+ZERO.OR.
00403      &           AVAIL(J,1,I).LT.-ZERO) THEN
00404                  WRITE(LU,*) 'ERROR IN LAYER CASE 5'
00405                  CALL PLANTE(1)
00406                  STOP
00407               ENDIF
00408             ENDDO
00409 !
00410 ! IF THERE REMAINS ONLY ONE LAYER
00411 ! -------------------------------
00412 !
00413           ELSE
00414 !           NOTE JMH: THE TRICKIEST PART...
00415 !           THE PROBLEM OF 0/0 CREATED BY THE CHOICE OF AVAIL
00416 !           AS MAIN VARIABLE...
00417             IF(ELAY%R(J)+EVOL.GT.1.D-15) THEN
00418               DO I=1,NSICLA
00419 !               AUX=AVAIL(J,1,I)
00420                 AVAIL(J,1,I) = (AVAIL(J,1,I)*ELAY%R(J)+
00421      &                          ZFCL_W%ADR(I)%P%R(J))
00422      &                          / (ELAY%R(J)+EVOL)
00423 !               IF(AVAIL(J,1,I).GT.1.D0+ZERO.OR.
00424 !    &            AVAIL(J,1,I).LT.-ZERO) THEN
00425 !                 WRITE(LU,*) 'ERROR IN LAYER CASE 6'
00426 !                 WRITE(LU,*) 'INITIAL AVAIL=',AUX
00427 !                 WRITE(LU,*) 'J=',J,' CLASS ',I
00428 !                 WRITE(LU,*) 'EVOL=',EVOL,' ELAY=',ELAY%R(J)
00429 !                 WRITE(LU,*) 'EVOL+ELAY=',EVOL+ELAY%R(J)
00430 !                 WRITE(LU,*) 'ZFCL=',ZFCL_W%ADR(I)%P%R(J)
00431 !                 WRITE(LU,*) 'DENOMINATOR=',ELAY%R(J)+EVOL
00432 !                 WRITE(LU,*) 'NUMERATOR=',AUX*ELAY%R(J)+
00433 !    &                                     ZFCL_W%ADR(I)%P%R(J)
00434 !               ENDIF
00435                 AVAIL(J,2,I) = 0.D0
00436               ENDDO
00437               IF(ELAY%R(J)+EVOL.LT.1.D-5) THEN
00438 !               PLAYING WITH ZEROES, RISK OF SUM NOT EQUAL TO 1
00439 !               ONLY BECAUSE OF TRUNCATION ERRORS, WE NORMALISE
00440                 TEST1=0.D0
00441                 DO I=1,NSICLA
00442                   AVAIL(J,1,I)=MAX(0.D0,MIN(1.D0,AVAIL(J,1,I)))
00443                   TEST1=TEST1+AVAIL(J,1,I)
00444                 ENDDO
00445                 IF((TEST1-1.D0)**2.GT.ZERO) THEN
00446                   DO I=1,NSICLA
00447                     AVAIL(J,1,I)=AVAIL(J,1,I)/MAX(TEST1,1.D-21)
00448                   ENDDO
00449                 ENDIF
00450               ENDIF
00451             ELSE
00452 !             JMH 11/04/2011
00453 !             ACLADM BEING RECOMPUTED WITH AVAIL
00454               AVAIL(J,1,1)=1.D0
00455               AVAIL(J,2,1)=1.D0
00456               DO I=2,NSICLA
00457                 AVAIL(J,1,I) = 0.D0
00458                 AVAIL(J,2,I) = 0.D0
00459               ENDDO
00460             ENDIF
00461             ELAY%R(J) = HEIGH
00462             ESTRATNEW(J) = 0.D0
00463             NLAYNEW(J) = 1
00464           ENDIF
00465         ENDIF
00466 !
00467       NLAYER%I(J) = NLAYNEW(J)
00468       ESTRAT%R(J) = ESTRATNEW(J)
00469       ES(J,1) = ELAY%R(J)
00470 !     CORRECTION JMH 12/04/2011: IN CASE OF RESTART, ES(J,2)
00471 !     WILL BE STORED IN A FILE AND LOOKED AT TO COUNT THE
00472 !     NUMBER OF LAYERS
00473 !     IF(NLAYER%I(J).GT.1) ES(J,2) = ESTRAT%R(J)
00474       ES(J,2) = ESTRAT%R(J)
00475 !
00476       TEST1 = 0.D0
00477       TEST2 = 0.D0
00478 !
00479       DO I=1,NSICLA
00480         DO K=1,NLAYER%I(J)
00481 !         CHECKS THAT AVAIL IS IN THE RANGE (-ZERO,1+ZERO)
00482           IF(AVAIL(J,K,I).GT.1.D0+ZERO.OR.AVAIL(J,K,I).LT.-ZERO) THEN
00483             WRITE(LU,*) 'ERROR ON FRACTIONS'
00484             WRITE(LU,*) 'LAYER ',K,' CLASS ',I,' POINT ',J
00485             IF(AVAIL(J,K,I).LT.0.D0) THEN
00486               WRITE(LU,*) 'AVAIL=' ,AVAIL(J,K,I)
00487             ELSE
00488               WRITE(LU,*) 'AVAIL-1=' ,AVAIL(J,K,I)-1.D0
00489             ENDIF
00490             WRITE(LU,*) 'ZFCL=',ZFCL_W%ADR(I)%P%R(J)
00491             WRITE(LU,*) 'EVOL=',EVOL,' ELAY=',ELAY%R(J)
00492             ARRET=1
00493           ELSE
00494 !           ONCE CHECKED THAT WE HAVE ONLY TRUNCATION ERRORS, CLIPS
00495             AVAIL(J,1,I)=MAX(AVAIL(J,1,I),0.D0)
00496             AVAIL(J,1,I)=MIN(AVAIL(J,1,I),1.D0)
00497           ENDIF
00498         ENDDO
00499         TEST1 = TEST1 + AVAIL(J,1,I)
00500         TEST2 = TEST2 + AVAIL(J,2,I)
00501       ENDDO
00502 !
00503 !     CHECKS THAT SUM OF AVAIL IS 1 FOR FIRST 2 LAYERS
00504 !
00505       IF(TEST1.GT.ZERO.AND.(TEST1-1.D0)**2>ZERO) THEN
00506         WRITE(LU,*) ' PROBLEM IN LAYER: J,TEST1',J,TEST1
00507         WRITE(LU,*) ' IN LAYER 1 SUM OF FRACTIONS NOT 1'
00508         ARRET=1
00509       ENDIF
00510       IF(TEST2.GT.ZERO.AND.(TEST2-1.D0)**2>ZERO) THEN
00511         WRITE(LU,*) ' PROBLEM IN LAYER: J,TEST2',J,TEST2
00512         WRITE(LU,*) ' IN LAYER 2 SUM OF FRACTIONS IS NOT 1'
00513         ARRET=1
00514       ENDIF
00515 !
00516 !     END OF LOOP ON ALL POINTS
00517 !
00518       ENDDO
00519 !
00520 !     COMPUTES THE TOTAL VOLUME OF SEDIMENTS IN THE DOMAIN
00521 !
00522       DO I = 1, NSICLA
00523         VOLTOT(I) = 0.D0
00524       ENDDO
00525       DO I=1,NSICLA
00526         DO J=1,NPOIN
00527           DO K=1,NLAYER%I(J)
00528             VOLTOT(I) = VOLTOT(I) + ES(J,K)*AVAIL(J,K,I)*MASBAS%R(J)
00529           ENDDO
00530         ENDDO
00531       ENDDO
00532 !
00533       IF(NCSIZE.GT.1) THEN
00534         DO I=1,NSICLA
00535           VOLTOT(I) = P_DSUM(VOLTOT(I))
00536         ENDDO
00537       ENDIF
00538 !
00539 !-----------------------------------------------------------------------
00540 !
00541 !     CLEAN STOP FOR ALL PROCESSORS IF PROBLEM
00542 !
00543       ARRET2=ARRET
00544       IF(NCSIZE.GT.1) ARRET2=P_ISUM(ARRET)
00545       IF(ARRET2.GT.0) THEN
00546         IF(LNG.EQ.1) WRITE(LU,*) 'ARRET APRES ERREUR DANS LAYER'
00547         IF(LNG.EQ.2) WRITE(LU,*) 'STOP AFTER AN ERROR IN LAYER'
00548         IF(ARRET.EQ.0) THEN
00549           IF(LNG.EQ.1) WRITE(LU,*) 'DANS ',ARRET2,' PROCESSEUR(S)'
00550           IF(LNG.EQ.2) WRITE(LU,*) 'IN ',ARRET2,' PROCESSOR(S)'
00551         ENDIF
00552         CALL PLANTE(1)
00553         STOP
00554       ENDIF
00555 !
00556 !-----------------------------------------------------------------------
00557 !
00558       RETURN
00559       END

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