murd3d.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac3d\murd3d.f
00002 !
00110                      SUBROUTINE MURD3D
00111 !                    *****************
00112 !
00113      &(FC,FN,VOLU,VOLUN,VOLU2,SVOLU2,DB,XB,DIM1XB,
00114      & TRA01,TRA02,TRA03,STRA01,STRA02,STRA03,IKLE3,MESH3,
00115      & NELEM3,NPOIN3,DT,SCHCF,LV,MSK,MASKEL,INFOR,CALFLU,FLUX,FLUEXT,
00116      & S0F,NSCE,SOURCES,FSCE,RAIN,PLUIE,TRAIN,NPOIN2,MINFC,MAXFC,MASKPT,
00117      & OPTBAN,FLODEL,FLOPAR,GLOSEG,DIMGLO,NSEG,NPLAN,IELM3)
00118 !
00119 !***********************************************************************
00120 ! TELEMAC3D   V6P2                                   21/08/2010
00121 !***********************************************************************
00122 !
00123 !
00124 !
00125 !
00126 !
00127 !
00128 !
00129 !
00130 !
00131 !
00132 !
00133 !
00134 !
00135 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00136 !| CALFLU         |-->| INDICATE IF FLUX IS CALCULATED FOR BALANCE
00137 !| DB             |<->| NOT SYMMETRIC MURD MATRIX OPTION N
00138 !| DIMGLO         |-->| FIRST DIMENSION OF GLOSEG
00139 !| DIM1XB         |-->| FIRST DIMENSION OF XB
00140 !| DT             |-->| TIME STEP
00141 !| FC             |<->| VARIABLE AFTER CONVECTION
00142 !| FLODEL         |-->| FLUX BY MESH EDGES
00143 !| FLOPAR         |-->| FLUXES BY SEGMENT, ASSEMBLED IN PARALLEL
00144 !| FLUEXT         |-->| OUTPUT FLUX BY NODE
00145 !| FLUX           |<->| FLUXES TO BE CHANGED
00146 !| FN             |-->| VARIABLE AT TIME N
00147 !| FSCE           |-->| SOURCE
00148 !| GLOSEG         |-->| FIRST AND SECOND POINT OF SEGMENTS
00149 !| IELM3          |-->| TYPE OF ELEMENT (41:PRISM, ETC.)
00150 !| IKLE3          |-->| GLOBAL 3D CONNECTIVITY
00151 !| INFOR          |-->| INFORMATIONS FOR SOLVERS
00152 !| LV             |-->| VECTOR LENGTH OF THE MACHINE
00153 !| MASKEL         |-->| MASKING OF ELEMENTS
00154 !|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
00155 !| MASKPT         |-->| MASKING PER POINT.
00156 !|                |   | =1. : NORMAL   =0. : MASKED
00157 !| MAXFC          |<->| MAXIMUM VALUE FOR FC
00158 !| MESH3          |<->| 3D MESH
00159 !| MINFC          |<->| MINIMUM VALUE FOR FC
00160 !| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
00161 !| NELEM3         |-->| NUMBER OF ELEMENTS IN 3D
00162 !| NPLAN          |-->| NUMBER OF PLANES IN THE 3D MESH OF PRISMS
00163 !| NPOIN2         |-->| NUMBER OF POINTS IN 2D
00164 !| NPOIN3         |-->| NUMBER OF 3D POINTS
00165 !| NSCE           |-->| NUMBER OF GIVEN POINTS FOR SOURCES
00166 !| NSEG           |-->| NUMBER OF SEGMENTS
00167 !| OPTBAN         |-->| OPTION FOR TIDAL FLATS, IF 1, FREE SURFACE
00168 !|                |   | MODIFIED AND PIECE-WISE LINEAR
00169 !| PLUIE          |-->| RAIN IN M/S MULTIPLIED BY VOLU2D
00170 !| RAIN           |-->| IF YES, THERE IS RAIN OR EVAPORATION
00171 !| S0F            |-->| EXPLICIT SOURCE TERM
00172 !| SCHCF          |-->| ADVECTION SCHEME FOR F
00173 !| SOURCES        |-->| SOURCES
00174 !| STRA01         |<->| STRUCTURE OF TRA01
00175 !| STRA02         |<->| STRUCTURE OF TRA02
00176 !| STRA03         |<->| STRUCTURE OF TRA03
00177 !| SVOLU2         |-->| STRUCTURE OF VOLU2
00178 !| TRA01          |<->| WORK ARRAY OF DIMENSION NPOIN3 EQUIVALENT TO
00179 !|                |   | VOLU2 FOR CURRENT FINAL TIME
00180 !| TRA02          |<->| WORK ARRAY
00181 !| TRA03          |<->| WORK ARRAY
00182 !| TRAIN          |-->| VALUE OF TRACER IN RAIN
00183 !| VOLU           |-->| CONTROL VOLUME AT TIME N+1
00184 !| VOLU2          |<->| LIKE VOLU, BUT ASSEMBLED IN PARALLEL
00185 !| VOLUN          |-->| CONTROL VOLUME AT TIME N
00186 !| XB             |<->| NOT SYMMETRIC MURD MATRIX OPTION N
00187 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00188 !
00189       USE BIEF
00190       USE DECLARATIONS_TELEMAC
00191 !
00192       IMPLICIT NONE
00193       INTEGER LNG,LU
00194       COMMON/INFO/LNG,LU
00195 !
00196 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00197 !
00198       INTEGER, INTENT(IN)             :: SCHCF,NELEM3,NPOIN3,LV,NPOIN2
00199       INTEGER, INTENT(IN)             :: IELM3,DIM1XB
00200 !                                                     6 OR 4
00201       INTEGER, INTENT(IN)             :: IKLE3(NELEM3,*),NSCE,OPTBAN
00202       INTEGER, INTENT(IN)             :: NSEG,NPLAN,DIMGLO
00203       INTEGER, INTENT(IN)             :: GLOSEG(DIMGLO,2)
00204 !
00205       DOUBLE PRECISION, INTENT(INOUT) :: FC(NPOIN3)
00206       DOUBLE PRECISION, INTENT(IN)    :: FN(NPOIN3),PLUIE(NPOIN2)
00207       DOUBLE PRECISION, INTENT(IN)    :: MASKEL(NELEM3), FLUEXT(NPOIN3)
00208 !
00209       DOUBLE PRECISION, INTENT(IN)    :: VOLUN(NPOIN3), VOLU(NPOIN3)
00210       DOUBLE PRECISION, INTENT(INOUT) :: VOLU2(NPOIN3)
00211       DOUBLE PRECISION, INTENT(INOUT) :: TRA01(NPOIN3)
00212       DOUBLE PRECISION, INTENT(INOUT) :: TRA02(NPOIN3), TRA03(NPOIN3)
00213       DOUBLE PRECISION, INTENT(INOUT) :: FLUX
00214       DOUBLE PRECISION, INTENT(IN)    :: FSCE(NSCE),MASKPT(NPOIN3)
00215       DOUBLE PRECISION, INTENT(IN)    :: DT,TRAIN
00216 !
00217       TYPE(BIEF_OBJ), INTENT(INOUT)   :: SVOLU2,MINFC,MAXFC
00218       TYPE(BIEF_OBJ), INTENT(IN)      :: SOURCES,S0F
00219       TYPE(BIEF_OBJ), INTENT(INOUT)   :: STRA01,STRA02,STRA03
00220       TYPE(BIEF_MESH), INTENT(INOUT)  :: MESH3
00221       DOUBLE PRECISION, INTENT(INOUT) :: DB(NPOIN3),XB(DIM1XB,NELEM3)
00222 !
00223 !     DIMENSION OF FLODEL AND FLOPAR=NSEG2D*NPLAN+NPOIN2*NETAGE
00224       DOUBLE PRECISION, INTENT(IN)    :: FLODEL(*),FLOPAR(*)
00225 !
00226       LOGICAL, INTENT(IN)             :: MSK,INFOR,CALFLU,RAIN
00227 !
00228 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00229 !
00230       DOUBLE PRECISION DTJ,PHIP,PHIM,ALFA,C,ALFA2,DTJALFA,MINEL,MAXEL
00231       DOUBLE PRECISION M12,M13,M14,M15,M16,M23,M24,M25,M26,M34
00232       DOUBLE PRECISION M35,M36,M45,M46,M56,T1,T2,T3,T4,T5,T6
00233       DOUBLE PRECISION F1MF2,F1MF3,F1MF4,F2MF3,F2MF4,F3MF4
00234 !
00235       INTEGER IELEM,IPOIN,NITER,IS,IGUILT,IIS
00236       INTEGER I1,I2,I3,I4,I5,I6,OPT,ISEG3D,NSEGH,NSEGV
00237 !
00238 !-----------------------------------------------------------------------
00239 !
00240       DOUBLE PRECISION P_DMIN,P_DSUM
00241       EXTERNAL P_DMIN,P_DSUM
00242 !
00243       DOUBLE PRECISION EPS
00244       DATA EPS /1.D-6/
00245 !     DATA EPS /10.D0/
00246 !
00247       INTEGER TOTITER
00248       DATA TOTITER /0/
00249 !
00250 !***********************************************************************
00251 !
00252       CALL CPSTVC(SVOLU2,STRA01)
00253       CALL CPSTVC(SVOLU2,STRA02)
00254       CALL CPSTVC(SVOLU2,STRA03)
00255       CALL CPSTVC(SVOLU2,MINFC)
00256       CALL CPSTVC(SVOLU2,MAXFC)
00257 !
00258       NITER = 0
00259       DTJ = DT
00260 !
00261       CALL OV ('X=Y     ',FC,FN,FN,C,NPOIN3)
00262 !
00263       CALL OV ('X=Y     ',TRA01, VOLUN, VOLUN, C, NPOIN3)
00264       IF (NCSIZE.GT.1) CALL PARCOM(STRA01,2,MESH3)
00265 !
00266       CALL OV ('X=Y     ',VOLU2, VOLU, VOLU, C, NPOIN3)
00267       IF (NCSIZE.GT.1) CALL PARCOM(SVOLU2,2,MESH3)
00268 !
00269 10    CONTINUE
00270 !
00271 !-----------------------------------------------------------------------
00272 !
00273 !  BUILDS THE PSI SCHEME FROM THE N SCHEME IF SCHCF=5
00274 !  SEE "HYDRODYNAMICS OF FREE SURFACE FLOWS" PAGE 193
00275 !
00276       IF(SCHCF.EQ.ADV_PSI) THEN
00277 !
00278         DO IPOIN=1,NPOIN3
00279           TRA02(IPOIN)=0.D0
00280         ENDDO
00281 !
00282         IF(IELM3.EQ.41) THEN
00283 !
00284         DO IELEM = 1,NELEM3
00285 !
00286           I1 = IKLE3(IELEM,1)
00287           I2 = IKLE3(IELEM,2)
00288           I3 = IKLE3(IELEM,3)
00289           I4 = IKLE3(IELEM,4)
00290           I5 = IKLE3(IELEM,5)
00291           I6 = IKLE3(IELEM,6)
00292 !
00293           T1 = FC(I1)
00294           T2 = FC(I2)
00295           T3 = FC(I3)
00296           T4 = FC(I4)
00297           T5 = FC(I5)
00298           T6 = FC(I6)
00299 !
00300           M12 = (XB(01,IELEM)-XB(16,IELEM)) * (T1-T2)
00301           M13 = (XB(02,IELEM)-XB(17,IELEM)) * (T1-T3)
00302           M14 = (XB(03,IELEM)-XB(18,IELEM)) * (T1-T4)
00303           M15 = (XB(04,IELEM)-XB(19,IELEM)) * (T1-T5)
00304           M16 = (XB(05,IELEM)-XB(20,IELEM)) * (T1-T6)
00305           M23 = (XB(06,IELEM)-XB(21,IELEM)) * (T2-T3)
00306           M24 = (XB(07,IELEM)-XB(22,IELEM)) * (T2-T4)
00307           M25 = (XB(08,IELEM)-XB(23,IELEM)) * (T2-T5)
00308           M26 = (XB(09,IELEM)-XB(24,IELEM)) * (T2-T6)
00309           M34 = (XB(10,IELEM)-XB(25,IELEM)) * (T3-T4)
00310           M35 = (XB(11,IELEM)-XB(26,IELEM)) * (T3-T5)
00311           M36 = (XB(12,IELEM)-XB(27,IELEM)) * (T3-T6)
00312           M45 = (XB(13,IELEM)-XB(28,IELEM)) * (T4-T5)
00313           M46 = (XB(14,IELEM)-XB(29,IELEM)) * (T4-T6)
00314           M56 = (XB(15,IELEM)-XB(30,IELEM)) * (T5-T6)
00315 !
00316           PHIP = MAX( M12,0.D0) + MAX( M13,0.D0) + MAX( M14,0.D0)
00317      &         + MAX( M15,0.D0) + MAX( M16,0.D0) + MAX( M23,0.D0)
00318      &         + MAX( M24,0.D0) + MAX( M25,0.D0) + MAX( M26,0.D0)
00319      &         + MAX( M34,0.D0) + MAX( M35,0.D0) + MAX( M36,0.D0)
00320      &         + MAX( M45,0.D0) + MAX( M46,0.D0) + MAX( M56,0.D0)
00321           PHIM = MAX(-M12,0.D0) + MAX(-M13,0.D0) + MAX(-M14,0.D0)
00322      &         + MAX(-M15,0.D0) + MAX(-M16,0.D0) + MAX(-M23,0.D0)
00323      &         + MAX(-M24,0.D0) + MAX(-M25,0.D0) + MAX(-M26,0.D0)
00324      &         + MAX(-M34,0.D0) + MAX(-M35,0.D0) + MAX(-M36,0.D0)
00325      &         + MAX(-M45,0.D0) + MAX(-M46,0.D0) + MAX(-M56,0.D0)
00326 !
00327           IF(PHIP.GE.PHIM) THEN
00328             ALFA = (PHIP - PHIM) / MAX(PHIP,1.D-10)
00329             IF(T2.GT.T1) THEN
00330               TRA02(I2)=TRA02(I2)+XB(16,IELEM)*ALFA*(FC(I1)-FC(I2))
00331             ELSE
00332               TRA02(I1)=TRA02(I1)+XB(01,IELEM)*ALFA*(FC(I2)-FC(I1))
00333             ENDIF
00334             IF(T3.GT.T1) THEN
00335               TRA02(I3)=TRA02(I3)+XB(17,IELEM)*ALFA*(FC(I1)-FC(I3))
00336             ELSE
00337               TRA02(I1)=TRA02(I1)+XB(02,IELEM)*ALFA*(FC(I3)-FC(I1))
00338             ENDIF
00339             IF(T4.GT.T1) THEN
00340               TRA02(I4)=TRA02(I4)+XB(18,IELEM)*ALFA*(FC(I1)-FC(I4))
00341             ELSE
00342               TRA02(I1)=TRA02(I1)+XB(03,IELEM)*ALFA*(FC(I4)-FC(I1))
00343             ENDIF
00344             IF(T5.GT.T1) THEN
00345               TRA02(I5)=TRA02(I5)+XB(19,IELEM)*ALFA*(FC(I1)-FC(I5))
00346             ELSE
00347               TRA02(I1)=TRA02(I1)+XB(04,IELEM)*ALFA*(FC(I5)-FC(I1))
00348             ENDIF
00349             IF(T6.GT.T1) THEN
00350               TRA02(I6)=TRA02(I6)+XB(20,IELEM)*ALFA*(FC(I1)-FC(I6))
00351             ELSE
00352               TRA02(I1)=TRA02(I1)+XB(05,IELEM)*ALFA*(FC(I6)-FC(I1))
00353             ENDIF
00354             IF(T3.GT.T2) THEN
00355               TRA02(I3)=TRA02(I3)+XB(21,IELEM)*ALFA*(FC(I2)-FC(I3))
00356             ELSE
00357               TRA02(I2)=TRA02(I2)+XB(06,IELEM)*ALFA*(FC(I3)-FC(I2))
00358             ENDIF
00359             IF(T4.GT.T2) THEN
00360               TRA02(I4)=TRA02(I4)+XB(22,IELEM)*ALFA*(FC(I2)-FC(I4))
00361             ELSE
00362               TRA02(I2)=TRA02(I2)+XB(07,IELEM)*ALFA*(FC(I4)-FC(I2))
00363             ENDIF
00364             IF(T5.GT.T2) THEN
00365               TRA02(I5)=TRA02(I5)+XB(23,IELEM)*ALFA*(FC(I2)-FC(I5))
00366             ELSE
00367               TRA02(I2)=TRA02(I2)+XB(08,IELEM)*ALFA*(FC(I5)-FC(I2))
00368             ENDIF
00369             IF(T6.GT.T2) THEN
00370               TRA02(I6)=TRA02(I6)+XB(24,IELEM)*ALFA*(FC(I2)-FC(I6))
00371             ELSE
00372               TRA02(I2)=TRA02(I2)+XB(09,IELEM)*ALFA*(FC(I6)-FC(I2))
00373             ENDIF
00374             IF(T4.GT.T3) THEN
00375               TRA02(I4)=TRA02(I4)+XB(25,IELEM)*ALFA*(FC(I3)-FC(I4))
00376             ELSE
00377               TRA02(I3)=TRA02(I3)+XB(10,IELEM)*ALFA*(FC(I4)-FC(I3))
00378             ENDIF
00379             IF(T5.GT.T3) THEN
00380               TRA02(I5)=TRA02(I5)+XB(26,IELEM)*ALFA*(FC(I3)-FC(I5))
00381             ELSE
00382               TRA02(I3)=TRA02(I3)+XB(11,IELEM)*ALFA*(FC(I5)-FC(I3))
00383             ENDIF
00384             IF(T6.GT.T3) THEN
00385               TRA02(I6)=TRA02(I6)+XB(27,IELEM)*ALFA*(FC(I3)-FC(I6))
00386             ELSE
00387               TRA02(I3)=TRA02(I3)+XB(12,IELEM)*ALFA*(FC(I6)-FC(I3))
00388             ENDIF
00389             IF(T5.GT.T4) THEN
00390               TRA02(I5)=TRA02(I5)+XB(28,IELEM)*ALFA*(FC(I4)-FC(I5))
00391             ELSE
00392               TRA02(I4)=TRA02(I4)+XB(13,IELEM)*ALFA*(FC(I5)-FC(I4))
00393             ENDIF
00394             IF(T6.GT.T4) THEN
00395               TRA02(I6)=TRA02(I6)+XB(29,IELEM)*ALFA*(FC(I4)-FC(I6))
00396             ELSE
00397               TRA02(I4)=TRA02(I4)+XB(14,IELEM)*ALFA*(FC(I6)-FC(I4))
00398             ENDIF
00399             IF(T6.GT.T5) THEN
00400               TRA02(I6)=TRA02(I6)+XB(30,IELEM)*ALFA*(FC(I5)-FC(I6))
00401             ELSE
00402               TRA02(I5)=TRA02(I5)+XB(15,IELEM)*ALFA*(FC(I6)-FC(I5))
00403             ENDIF
00404           ELSE
00405             ALFA = (PHIM - PHIP) / MAX(PHIM,1.D-10)
00406             IF(T2.GT.T1) THEN
00407               TRA02(I1)=TRA02(I1)+XB(01,IELEM)*ALFA*(FC(I2)-FC(I1))
00408             ELSE
00409               TRA02(I2)=TRA02(I2)+XB(16,IELEM)*ALFA*(FC(I1)-FC(I2))
00410             ENDIF
00411             IF(T3.GT.T1) THEN
00412               TRA02(I1)=TRA02(I1)+XB(02,IELEM)*ALFA*(FC(I3)-FC(I1))
00413             ELSE
00414               TRA02(I3)=TRA02(I3)+XB(17,IELEM)*ALFA*(FC(I1)-FC(I3))
00415             ENDIF
00416             IF(T4.GT.T1) THEN
00417               TRA02(I1)=TRA02(I1)+XB(03,IELEM)*ALFA*(FC(I4)-FC(I1))
00418             ELSE
00419               TRA02(I4)=TRA02(I4)+XB(18,IELEM)*ALFA*(FC(I1)-FC(I4))
00420             ENDIF
00421             IF(T5.GT.T1) THEN
00422               TRA02(I1)=TRA02(I1)+XB(04,IELEM)*ALFA*(FC(I5)-FC(I1))
00423             ELSE
00424               TRA02(I5)=TRA02(I5)+XB(19,IELEM)*ALFA*(FC(I1)-FC(I5))
00425             ENDIF
00426             IF(T6.GT.T1) THEN
00427               TRA02(I1)=TRA02(I1)+XB(05,IELEM)*ALFA*(FC(I6)-FC(I1))
00428             ELSE
00429               TRA02(I6)=TRA02(I6)+XB(20,IELEM)*ALFA*(FC(I1)-FC(I6))
00430             ENDIF
00431             IF(T3.GT.T2) THEN
00432               TRA02(I2)=TRA02(I2)+XB(06,IELEM)*ALFA*(FC(I3)-FC(I2))
00433             ELSE
00434               TRA02(I3)=TRA02(I3)+XB(21,IELEM)*ALFA*(FC(I2)-FC(I3))
00435             ENDIF
00436             IF(T4.GT.T2) THEN
00437               TRA02(I2)=TRA02(I2)+XB(07,IELEM)*ALFA*(FC(I4)-FC(I2))
00438             ELSE
00439               TRA02(I4)=TRA02(I4)+XB(22,IELEM)*ALFA*(FC(I2)-FC(I4))
00440             ENDIF
00441             IF(T5.GT.T2) THEN
00442               TRA02(I2)=TRA02(I2)+XB(08,IELEM)*ALFA*(FC(I5)-FC(I2))
00443             ELSE
00444               TRA02(I5)=TRA02(I5)+XB(23,IELEM)*ALFA*(FC(I2)-FC(I5))
00445             ENDIF
00446             IF(T6.GT.T2) THEN
00447               TRA02(I2)=TRA02(I2)+XB(09,IELEM)*ALFA*(FC(I6)-FC(I2))
00448             ELSE
00449               TRA02(I6)=TRA02(I6)+XB(24,IELEM)*ALFA*(FC(I2)-FC(I6))
00450             ENDIF
00451             IF(T4.GT.T3) THEN
00452               TRA02(I3)=TRA02(I3)+XB(10,IELEM)*ALFA*(FC(I4)-FC(I3))
00453             ELSE
00454               TRA02(I4)=TRA02(I4)+XB(25,IELEM)*ALFA*(FC(I3)-FC(I4))
00455             ENDIF
00456             IF(T5.GT.T3) THEN
00457               TRA02(I3)=TRA02(I3)+XB(11,IELEM)*ALFA*(FC(I5)-FC(I3))
00458             ELSE
00459               TRA02(I5)=TRA02(I5)+XB(26,IELEM)*ALFA*(FC(I3)-FC(I5))
00460             ENDIF
00461             IF(T6.GT.T3) THEN
00462               TRA02(I3)=TRA02(I3)+XB(12,IELEM)*ALFA*(FC(I6)-FC(I3))
00463             ELSE
00464               TRA02(I6)=TRA02(I6)+XB(27,IELEM)*ALFA*(FC(I3)-FC(I6))
00465             ENDIF
00466             IF(T5.GT.T4) THEN
00467               TRA02(I4)=TRA02(I4)+XB(13,IELEM)*ALFA*(FC(I5)-FC(I4))
00468             ELSE
00469               TRA02(I5)=TRA02(I5)+XB(28,IELEM)*ALFA*(FC(I4)-FC(I5))
00470             ENDIF
00471             IF(T6.GT.T4) THEN
00472               TRA02(I4)=TRA02(I4)+XB(14,IELEM)*ALFA*(FC(I6)-FC(I4))
00473             ELSE
00474               TRA02(I6)=TRA02(I6)+XB(29,IELEM)*ALFA*(FC(I4)-FC(I6))
00475             ENDIF
00476             IF(T6.GT.T5) THEN
00477               TRA02(I5)=TRA02(I5)+XB(15,IELEM)*ALFA*(FC(I6)-FC(I5))
00478             ELSE
00479               TRA02(I6)=TRA02(I6)+XB(30,IELEM)*ALFA*(FC(I5)-FC(I6))
00480             ENDIF
00481           ENDIF
00482 !
00483         ENDDO ! IELEM
00484 !
00485         ELSEIF(IELM3.EQ.51) THEN
00486 !
00487         DO IELEM = 1,NELEM3
00488 !
00489           I1 = IKLE3(IELEM,1)
00490           I2 = IKLE3(IELEM,2)
00491           I3 = IKLE3(IELEM,3)
00492           I4 = IKLE3(IELEM,4)
00493 !
00494           F1MF2 = FC(I1)-FC(I2)
00495           F1MF3 = FC(I1)-FC(I3)
00496           F1MF4 = FC(I1)-FC(I4)
00497           F2MF3 = FC(I2)-FC(I3)
00498           F2MF4 = FC(I2)-FC(I4)
00499           F3MF4 = FC(I3)-FC(I4)
00500 !
00501           M12 = (XB(01,IELEM)-XB(07,IELEM)) * F1MF2
00502           M13 = (XB(02,IELEM)-XB(08,IELEM)) * F1MF3
00503           M14 = (XB(03,IELEM)-XB(09,IELEM)) * F1MF4
00504           M23 = (XB(04,IELEM)-XB(10,IELEM)) * F2MF3
00505           M24 = (XB(05,IELEM)-XB(11,IELEM)) * F2MF4
00506           M34 = (XB(06,IELEM)-XB(12,IELEM)) * F3MF4
00507 !
00508           PHIP = MAX( M12,0.D0) + MAX( M13,0.D0) + MAX( M14,0.D0)
00509      &         + MAX( M23,0.D0) + MAX( M24,0.D0) + MAX( M34,0.D0)
00510           PHIM = MAX(-M12,0.D0) + MAX(-M13,0.D0) + MAX(-M14,0.D0)
00511      &         + MAX(-M23,0.D0) + MAX(-M24,0.D0) + MAX(-M34,0.D0)
00512 !
00513           IF(PHIP.GE.PHIM) THEN
00514             ALFA = (PHIP - PHIM) / MAX(PHIP,1.D-10)
00515             IF(F1MF2.LT.0.D0) THEN
00516               TRA02(I2)=TRA02(I2)+XB(07,IELEM)*ALFA*F1MF2
00517             ELSE
00518               TRA02(I1)=TRA02(I1)-XB(01,IELEM)*ALFA*F1MF2
00519             ENDIF
00520             IF(F1MF3.LT.0.D0) THEN
00521               TRA02(I3)=TRA02(I3)+XB(08,IELEM)*ALFA*F1MF3
00522             ELSE
00523               TRA02(I1)=TRA02(I1)-XB(02,IELEM)*ALFA*F1MF3
00524             ENDIF
00525             IF(F1MF4.LT.0.D0) THEN
00526               TRA02(I4)=TRA02(I4)+XB(09,IELEM)*ALFA*F1MF4
00527             ELSE
00528               TRA02(I1)=TRA02(I1)-XB(03,IELEM)*ALFA*F1MF4
00529             ENDIF
00530             IF(F2MF3.LT.0.D0) THEN
00531               TRA02(I3)=TRA02(I3)+XB(10,IELEM)*ALFA*F2MF3
00532             ELSE
00533               TRA02(I2)=TRA02(I2)-XB(04,IELEM)*ALFA*F2MF3
00534             ENDIF
00535             IF(F2MF4.LT.0.D0) THEN
00536               TRA02(I4)=TRA02(I4)+XB(11,IELEM)*ALFA*F2MF4
00537             ELSE
00538               TRA02(I2)=TRA02(I2)-XB(05,IELEM)*ALFA*F2MF4
00539             ENDIF
00540             IF(F3MF4.LT.0.D0) THEN
00541               TRA02(I4)=TRA02(I4)+XB(12,IELEM)*ALFA*F3MF4
00542             ELSE
00543               TRA02(I3)=TRA02(I3)-XB(06,IELEM)*ALFA*F3MF4
00544             ENDIF
00545           ELSE
00546             ALFA = (PHIM - PHIP) / MAX(PHIM,1.D-10)
00547             IF(F1MF2.LT.0.D0) THEN
00548               TRA02(I1)=TRA02(I1)-XB(01,IELEM)*ALFA*F1MF2
00549             ELSE
00550               TRA02(I2)=TRA02(I2)+XB(07,IELEM)*ALFA*F1MF2
00551             ENDIF
00552             IF(F1MF3.LT.0.D0) THEN
00553               TRA02(I1)=TRA02(I1)-XB(02,IELEM)*ALFA*F1MF3
00554             ELSE
00555               TRA02(I3)=TRA02(I3)+XB(08,IELEM)*ALFA*F1MF3
00556             ENDIF
00557             IF(F1MF4.LT.0.D0) THEN
00558               TRA02(I1)=TRA02(I1)-XB(03,IELEM)*ALFA*F1MF4
00559             ELSE
00560               TRA02(I4)=TRA02(I4)+XB(09,IELEM)*ALFA*F1MF4
00561             ENDIF
00562             IF(F2MF3.LT.0.D0) THEN
00563               TRA02(I2)=TRA02(I2)-XB(04,IELEM)*ALFA*F2MF3
00564             ELSE
00565               TRA02(I3)=TRA02(I3)+XB(10,IELEM)*ALFA*F2MF3
00566             ENDIF
00567             IF(F2MF4.LT.0.D0) THEN
00568               TRA02(I2)=TRA02(I2)-XB(05,IELEM)*ALFA*F2MF4
00569             ELSE
00570               TRA02(I4)=TRA02(I4)+XB(11,IELEM)*ALFA*F2MF4
00571             ENDIF
00572             IF(F3MF4.LT.0.D0) THEN
00573               TRA02(I3)=TRA02(I3)-XB(06,IELEM)*ALFA*F3MF4
00574             ELSE
00575               TRA02(I4)=TRA02(I4)+XB(12,IELEM)*ALFA*F3MF4
00576             ENDIF
00577           ENDIF
00578 !
00579         ENDDO ! IELEM
00580 !
00581         ELSE
00582           WRITE(LU,*) 'ELEMENT ',IELM3,' NOT COMPUTED IN MURD3D'
00583           CALL PLANTE(1)
00584           STOP
00585         ENDIF
00586 !
00587 !-----------------------------------------------------------------------
00588 !
00589       ELSEIF(SCHCF.EQ.ADV_NSC) THEN
00590 !
00591 !  COMPUTES THE FLUX TERMS PER UNIT OF TIME:
00592 !
00593 !  I.E. : SUM ON J (LAMBDA(I,J)*(FC(J)-FC(I))
00594 !
00595 !  DECOMPOSED INTO : -(SUM ON J (LAMBDA(I,J)*(FC(I)) WHICH IS DB * FC(I)
00596 !                    +(SUM ON J (LAMBDA(I,J)*(FC(J)) WHICH IS XB * FC
00597 !
00598 !     CALL MV0606
00599 !    &('X=AY    ', TRA02, DB, 'Q', XB, 'Q', FC, C, IKLE3(1,1),
00600 !    & IKLE3(1,2),IKLE3(1,3),IKLE3(1,4),IKLE3(1,5),IKLE3(1,6),
00601 !    & NPOIN3, NELEM3, NELEM3,
00602 !    & W1(1,1),W1(1,2),W1(1,3),W1(1,4),W1(1,5),W1(1,6))
00603 !     CALL ASSVEC
00604 !    &(TRA02,IKLE3,NPOIN3,NELEM3,NELEM3,41,W1,.FALSE.,LV,MSK,MASKEL)
00605 !
00606 !     REPLACES CALL MV0606 AND CALL ASSVEC
00607 !
00608       CALL OV ('X=YZ    ',TRA02,DB,FC,C,NPOIN3)
00609 !
00610       IF(IELM3.EQ.41) THEN
00611 !
00612       DO IELEM = 1 , NELEM3
00613 !
00614         I1 = IKLE3(IELEM,1)
00615         I2 = IKLE3(IELEM,2)
00616         I3 = IKLE3(IELEM,3)
00617         I4 = IKLE3(IELEM,4)
00618         I5 = IKLE3(IELEM,5)
00619         I6 = IKLE3(IELEM,6)
00620 !
00621         TRA02(I1) = TRA02(I1) + XB(01,IELEM) * FC(I2)
00622      &                        + XB(02,IELEM) * FC(I3)
00623      &                        + XB(03,IELEM) * FC(I4)
00624      &                        + XB(04,IELEM) * FC(I5)
00625      &                        + XB(05,IELEM) * FC(I6)
00626 !
00627         TRA02(I2) = TRA02(I2) + XB(16,IELEM) * FC(I1)
00628      &                        + XB(06,IELEM) * FC(I3)
00629      &                        + XB(07,IELEM) * FC(I4)
00630      &                        + XB(08,IELEM) * FC(I5)
00631      &                        + XB(09,IELEM) * FC(I6)
00632 !
00633         TRA02(I3) = TRA02(I3) + XB(17,IELEM) * FC(I1)
00634      &                        + XB(21,IELEM) * FC(I2)
00635      &                        + XB(10,IELEM) * FC(I4)
00636      &                        + XB(11,IELEM) * FC(I5)
00637      &                        + XB(12,IELEM) * FC(I6)
00638 !
00639         TRA02(I4) = TRA02(I4) + XB(18,IELEM) * FC(I1)
00640      &                        + XB(22,IELEM) * FC(I2)
00641      &                        + XB(25,IELEM) * FC(I3)
00642      &                        + XB(13,IELEM) * FC(I5)
00643      &                        + XB(14,IELEM) * FC(I6)
00644 !
00645         TRA02(I5) = TRA02(I5) + XB(19,IELEM) * FC(I1)
00646      &                        + XB(23,IELEM) * FC(I2)
00647      &                        + XB(26,IELEM) * FC(I3)
00648      &                        + XB(28,IELEM) * FC(I4)
00649      &                        + XB(15,IELEM) * FC(I6)
00650 !
00651         TRA02(I6) = TRA02(I6) + XB(20,IELEM) * FC(I1)
00652      &                        + XB(24,IELEM) * FC(I2)
00653      &                        + XB(27,IELEM) * FC(I3)
00654      &                        + XB(29,IELEM) * FC(I4)
00655      &                        + XB(30,IELEM) * FC(I5)
00656 !
00657       ENDDO
00658 !
00659       ELSEIF(IELM3.EQ.51) THEN
00660 !
00661       DO IELEM = 1 , NELEM3
00662 !
00663         I1 = IKLE3(IELEM,1)
00664         I2 = IKLE3(IELEM,2)
00665         I3 = IKLE3(IELEM,3)
00666         I4 = IKLE3(IELEM,4)
00667 !
00668         TRA02(I1) = TRA02(I1) + XB(01,IELEM) * FC(I2)
00669      &                        + XB(02,IELEM) * FC(I3)
00670      &                        + XB(03,IELEM) * FC(I4)
00671 !
00672         TRA02(I2) = TRA02(I2) + XB(04,IELEM) * FC(I3)
00673      &                        + XB(05,IELEM) * FC(I4)
00674      &                        + XB(07,IELEM) * FC(I1)
00675 !
00676         TRA02(I3) = TRA02(I3) + XB(06,IELEM) * FC(I4)
00677      &                        + XB(08,IELEM) * FC(I1)
00678      &                        + XB(10,IELEM) * FC(I2)
00679 !
00680         TRA02(I4) = TRA02(I4) + XB(09,IELEM) * FC(I1)
00681      &                        + XB(11,IELEM) * FC(I2)
00682      &                        + XB(12,IELEM) * FC(I3)
00683 !
00684       ENDDO
00685 !
00686       ELSE
00687         WRITE(LU,*) 'ELEMENT ',IELM3,' NOT COMPUTED IN MURD3D'
00688         CALL PLANTE(1)
00689         STOP
00690       ENDIF
00691 !
00692 !     END OF THE REPLACEMENT OF CALL MV0606
00693 !
00694       ELSEIF(SCHCF.EQ.ADV_LPO) THEN
00695 !
00696         IF(IELM3.NE.41) THEN
00697           WRITE(LU,*) 'MURD3D: ELEMENT ',IELM3,' NOT IMPLEMENTED'
00698           WRITE(LU,*) '        WITH SCHEME ',SCHCF
00699           CALL PLANTE(1)
00700           STOP
00701         ENDIF
00702 !
00703         NSEGH=NSEG*NPLAN
00704         NSEGV=(NPLAN-1)*NPOIN2
00705 !
00706 !       COMPUTES DB AND TRA02 IN UPWIND EXPLICIT FINITE VOLUMES
00707 !       POSSIBLE OPTIMISATION: DB IS NOT COMPUTED IF OPT=2
00708 !
00709         DO IPOIN=1,NPOIN3
00710           DB(IPOIN)=0.D0
00711           TRA02(IPOIN)=0.D0
00712         ENDDO
00713 !
00714 !       HORIZONTAL AND VERTICAL FLUXES ONLY
00715 !       BEWARE : FLUXES FROM POINT 2 TO 1 IN SEGMENT
00716 !                WITH POINTS 1 AND 2 (SEE PRECON AND FLUX3D)
00717 !
00718         DO ISEG3D = 1,NSEGH+NSEGV
00719           I1=GLOSEG(ISEG3D,1)
00720           I2=GLOSEG(ISEG3D,2)
00721           IF(FLOPAR(ISEG3D).GT.0.D0) THEN
00722             DB(I1)   = DB(I1)   -FLODEL(ISEG3D)
00723             TRA02(I1)= TRA02(I1)-FLODEL(ISEG3D)*(FC(I1)-FC(I2))
00724           ELSEIF(FLOPAR(ISEG3D).LT.0.D0) THEN
00725             DB(I2)   = DB(I2)   +FLODEL(ISEG3D)
00726             TRA02(I2)= TRA02(I2)+FLODEL(ISEG3D)*(FC(I2)-FC(I1))
00727           ENDIF
00728         ENDDO
00729 !
00730       ELSE
00731         WRITE(LU,*) 'MURD3D: UNKNOWN ADVECTION SCHEME: ',SCHCF
00732         CALL PLANTE(1)
00733         STOP
00734       ENDIF
00735 !
00736       IF(NCSIZE.GT.1) CALL PARCOM(STRA02,2,MESH3)
00737 !
00738 !-----------------------------------------------------------------------
00739 !
00740 !  COMPUTES THE LIMITING SUB-TIMESTEP :
00741 !     MULTIPLIES THE REMAINING SUB TIMESTEP BY DTJ
00742 !     AND ADDS TO VOLU NEGATIVE TERM (MASS AT N+1)
00743 !
00744 !     OPT=1 : OLD CRITERION
00745 !     OPT=2 : NEW CRITERION LESS RESTRICTIVE
00746       OPT=2
00747 !
00748       IF(OPT.EQ.2) THEN
00749 !
00750 !     COMPUTES THE LOCAL EXTREMA
00751 !
00752       DO IPOIN=1,NPOIN3
00753         MINFC%R(IPOIN)=FC(IPOIN)
00754         MAXFC%R(IPOIN)=FC(IPOIN)
00755       ENDDO
00756 !
00757       IF(IELM3.EQ.41) THEN
00758 !
00759         DO IELEM=1,NELEM3
00760           I1=IKLE3(IELEM,1)
00761           I2=IKLE3(IELEM,2)
00762           I3=IKLE3(IELEM,3)
00763           I4=IKLE3(IELEM,4)
00764           I5=IKLE3(IELEM,5)
00765           I6=IKLE3(IELEM,6)
00766           MINEL=MIN(FC(I1),FC(I2),FC(I3),FC(I4),FC(I5),FC(I6))
00767           MAXEL=MAX(FC(I1),FC(I2),FC(I3),FC(I4),FC(I5),FC(I6))
00768           MINFC%R(I1)=MIN(MINFC%R(I1),MINEL)
00769           MINFC%R(I2)=MIN(MINFC%R(I2),MINEL)
00770           MINFC%R(I3)=MIN(MINFC%R(I3),MINEL)
00771           MINFC%R(I4)=MIN(MINFC%R(I4),MINEL)
00772           MINFC%R(I5)=MIN(MINFC%R(I5),MINEL)
00773           MINFC%R(I6)=MIN(MINFC%R(I6),MINEL)
00774           MAXFC%R(I1)=MAX(MAXFC%R(I1),MAXEL)
00775           MAXFC%R(I2)=MAX(MAXFC%R(I2),MAXEL)
00776           MAXFC%R(I3)=MAX(MAXFC%R(I3),MAXEL)
00777           MAXFC%R(I4)=MAX(MAXFC%R(I4),MAXEL)
00778           MAXFC%R(I5)=MAX(MAXFC%R(I5),MAXEL)
00779           MAXFC%R(I6)=MAX(MAXFC%R(I6),MAXEL)
00780         ENDDO
00781 !
00782       ELSEIF(IELM3.EQ.51) THEN
00783 !
00784         DO IELEM=1,NELEM3
00785           I1=IKLE3(IELEM,1)
00786           I2=IKLE3(IELEM,2)
00787           I3=IKLE3(IELEM,3)
00788           I4=IKLE3(IELEM,4)
00789           MINEL=MIN(FC(I1),FC(I2),FC(I3),FC(I4))
00790           MAXEL=MAX(FC(I1),FC(I2),FC(I3),FC(I4))
00791           MINFC%R(I1)=MIN(MINFC%R(I1),MINEL)
00792           MINFC%R(I2)=MIN(MINFC%R(I2),MINEL)
00793           MINFC%R(I3)=MIN(MINFC%R(I3),MINEL)
00794           MINFC%R(I4)=MIN(MINFC%R(I4),MINEL)
00795           MAXFC%R(I1)=MAX(MAXFC%R(I1),MAXEL)
00796           MAXFC%R(I2)=MAX(MAXFC%R(I2),MAXEL)
00797           MAXFC%R(I3)=MAX(MAXFC%R(I3),MAXEL)
00798           MAXFC%R(I4)=MAX(MAXFC%R(I4),MAXEL)
00799         ENDDO
00800 !
00801       ELSE
00802         WRITE(LU,*) 'ELEMENT ',IELM3,' NOT COMPUTED IN MURD3D'
00803         CALL PLANTE(1)
00804         STOP
00805       ENDIF
00806 !
00807 !     IN PARALLEL MODE: GLOBAL EXTREMA
00808 !
00809       IF(NCSIZE.GT.1) THEN
00810         CALL PARCOM(MAXFC,3,MESH3)
00811         CALL PARCOM(MINFC,4,MESH3)
00812       ENDIF
00813 !
00814 !     NEW COMPUTATION OF TRA03
00815 !
00816       DO IPOIN=1,NPOIN3
00817         IF(TRA02(IPOIN).GT.0.D0) THEN
00818           TRA03(IPOIN)=VOLU2(IPOIN)-DTJ*TRA02(IPOIN)/
00819      &                 MAX(MAXFC%R(IPOIN)-FC(IPOIN),1.D-12)
00820         ELSE
00821           TRA03(IPOIN)=VOLU2(IPOIN)+DTJ*TRA02(IPOIN)/
00822      &                 MAX(FC(IPOIN)-MINFC%R(IPOIN),1.D-12)
00823         ENDIF
00824       ENDDO
00825 !     POSITIVE SOURCES CHANGE THE MONOTONICITY CRITERION
00826       IF(NSCE.GT.0) THEN
00827         DO IS=1,NSCE
00828           DO IPOIN=1,NPOIN3
00829             IF(SOURCES%ADR(IS)%P%R(IPOIN).GT.0.D0) THEN
00830               TRA03(IPOIN)=TRA03(IPOIN)
00831      &                    -DTJ*SOURCES%ADR(IS)%P%R(IPOIN)
00832 !                                          WITH PARCOM
00833             ENDIF
00834           ENDDO
00835         ENDDO
00836       ENDIF
00837 !
00838       ELSEIF(OPT.EQ.1) THEN
00839 !
00840       IF(SCHCF.EQ.ADV_PSI) THEN
00841         WRITE(LU,*) 'MURD3D: OPT=1 IS NOT ALLOWED WITH PSI-SCHEME'
00842         WRITE(LU,*) '        BECAUSE DIAGONAL DB IS NOT COMPUTED'
00843         CALL PLANTE(1)
00844         STOP
00845       ENDIF
00846 !
00847 !     TRA03 : COEFFICIENT OF FC(I), THAT MUST REMAIN POSITIVE FOR MONOTONICITY
00848       CALL OV('X=Y+CZ  ',TRA03, VOLU, DB, DTJ, NPOIN3)
00849 !
00850 !     POSITIVE SOURCES CHANGE THE MONOTONICITY CRITERION
00851 !
00852       IF(NSCE.GT.0) THEN
00853         DO IS=1,NSCE
00854           DO IPOIN=1,NPOIN3
00855             IF(SOURCES%ADR(IS)%P%R(IPOIN).GT.0.D0) THEN
00856               TRA03(IPOIN)=TRA03(IPOIN)
00857      &                    -DTJ*SOURCES%ADR(IS+NSCE)%P%R(IPOIN)
00858 !                                          WITHOUT PARCOM
00859             ENDIF
00860           ENDDO
00861         ENDDO
00862       ENDIF
00863 !
00864       IF(NCSIZE.GT.1) CALL PARCOM(STRA03,2,MESH3)
00865 !
00866       ELSE
00867         WRITE(LU,*) 'MURD3D: OPT=',OPT,' NOT IMPLEMENTED'
00868         CALL PLANTE(1)
00869         STOP
00870       ENDIF
00871 !
00872 !     IF MONOTONICITY IS NOT ENSURED, REDUCTION OF TIME-STEP BY FACTOR ALFA
00873 !     THE MINIMUM ON ALL POINTS WILL BE TAKEN
00874 !
00875       ALFA = 1.D0
00876       IGUILT=1
00877       IF(OPTBAN.EQ.2) THEN
00878         DO IPOIN = 1,NPOIN3
00879           IF(TRA03(IPOIN).LT.0.D0.AND.MASKPT(IPOIN).GT.0.5D0) THEN
00880             IF(ABS(TRA01(IPOIN)-TRA03(IPOIN)).GT.EPS) THEN
00881 !             CONSIDERING THAT THE NEW TIME-STEP WILL BE ALFA*DTJ
00882 !             VOLU WILL BE AT THAT TIME ALFA*VOLU(N+1)+(1-ALFA)*VOLU(IN TRA01)
00883 !             MAXIMUM POSSIBLE ALFA IS SUCH THAT VOLU+DB*ALFA*DTJ=0
00884 !             HENCE THE FOLLOWING FORMULA :
00885               ALFA2=TRA01(IPOIN)/(TRA01(IPOIN)-TRA03(IPOIN))
00886               IF(ALFA.GT.ALFA2) THEN
00887                 ALFA = ALFA2
00888                 IGUILT=IPOIN
00889               ENDIF
00890             ENDIF
00891           ENDIF
00892         ENDDO
00893       ELSE
00894         DO IPOIN = 1,NPOIN3
00895 !                                          TIDAL FLATS : VOLUN=0
00896           IF(TRA03(IPOIN).LT.0.D0.AND.TRA01(IPOIN).GT.EPS) THEN
00897             IF(ABS(TRA01(IPOIN)-TRA03(IPOIN)).GT.EPS) THEN
00898 !             CONSIDERING THAT THE NEW TIME-STEP WILL BE ALFA*DTJ
00899 !             VOLU WILL BE AT THAT TIME ALFA*VOLU(N+1)+(1-ALFA)*VOLU(IN TRA01)
00900 !             MAXIMUM POSSIBLE ALFA IS SUCH THAT VOLU+DB*ALFA*DTJ=0
00901 !             HENCE THE FOLLOWING FORMULA :
00902               ALFA2=TRA01(IPOIN)/(TRA01(IPOIN)-TRA03(IPOIN))
00903               IF(ALFA.GT.ALFA2) THEN
00904                 ALFA = ALFA2
00905                 IGUILT=IPOIN
00906               ENDIF
00907             ENDIF
00908           ENDIF
00909         ENDDO
00910       ENDIF
00911 !
00912       IF(NCSIZE.GT.1) ALFA = P_DMIN(ALFA)
00913       DTJALFA=DTJ*ALFA
00914 !
00915 !     COMPUTES VOLU AFTER AN EXTRA ALFA*DTJ
00916 !
00917       CALL OV('X=CX    ',TRA01,TRA01,TRA01,1.D0-ALFA,NPOIN3)
00918       CALL OV('X=X+CY  ',TRA01,VOLU2,VOLU2,     ALFA,NPOIN3)
00919 !
00920       IF(CALFLU) THEN
00921         DO IPOIN = 1,NPOIN3
00922           FLUX = FLUX + DTJALFA*FC(IPOIN)*FLUEXT(IPOIN)
00923         ENDDO
00924         IF(NSCE.GT.0) THEN
00925           DO IS=1,NSCE
00926             IIS=IS
00927 !           HERE IN PARALLEL SOURCES WITHOUT PARCOM
00928 !           ARE STORED AT ADRESSES IS+NSCE (SEE SOURCES_SINKS.F)
00929             IF(NCSIZE.GT.1) IIS=IIS+NSCE
00930             DO IPOIN=1,NPOIN3
00931               IF(SOURCES%ADR(IS)%P%R(IPOIN).GT.0.D0) THEN
00932                 FLUX=FLUX
00933      &              -DTJALFA*FSCE(IS)*SOURCES%ADR(IIS)%P%R(IPOIN)
00934               ELSE
00935                 FLUX=FLUX
00936      &              -DTJALFA*FC(IPOIN)*SOURCES%ADR(IIS)%P%R(IPOIN)
00937               ENDIF
00938             ENDDO
00939           ENDDO
00940         ENDIF
00941       ENDIF
00942 !
00943 !     ADVECTION DURING ALFA*DTJ
00944 !
00945 !     SOURCES (BUT WHEN INTAKE, FSCE=FC)
00946 !
00947       IF(NSCE.GT.0) THEN
00948         DO IS=1,NSCE
00949           IF(OPTBAN.EQ.2) THEN
00950             DO IPOIN=1,NPOIN3
00951               IF(MASKPT(IPOIN).GT.0.5D0.AND.TRA01(IPOIN).GT.EPS) THEN
00952                 FC(IPOIN)=FC(IPOIN)+DTJALFA*(FSCE(IS)-FC(IPOIN))
00953      &          *MAX(SOURCES%ADR(IS)%P%R(IPOIN),0.D0)/TRA01(IPOIN)
00954               ENDIF
00955             ENDDO
00956           ELSE
00957             DO IPOIN=1,NPOIN3
00958               IF(TRA01(IPOIN).GT.EPS) THEN
00959                 FC(IPOIN)=FC(IPOIN)+DTJALFA*(FSCE(IS)-FC(IPOIN))
00960      &          *MAX(SOURCES%ADR(IS)%P%R(IPOIN),0.D0)/TRA01(IPOIN)
00961               ENDIF
00962             ENDDO
00963           ENDIF
00964         ENDDO
00965       ENDIF
00966 !
00967 !     RAIN (NOTE: SHOULD BE TAKEN INTO ACCOUNT IN STABILITY CRITERION)
00968 !     VALUE OF TRACER IN RAIN TAKEN INTO ACCOUNT ONLY IF RAIN POSITIVE
00969 !     NOT IN CASE OF EVAPORATION, HENCE THE MAX(PLUIE,0)
00970 !
00971       IF(RAIN) THEN
00972         IF(OPTBAN.EQ.2) THEN
00973           DO IPOIN=1,NPOIN2
00974             IS=NPOIN3-NPOIN2+IPOIN
00975             IF(MASKPT(IPOIN).GT.0.5D0.AND.TRA01(IS).GT.EPS) THEN
00976               C=TRAIN*MAX(PLUIE(IPOIN),0.D0)-FC(IS)*PLUIE(IPOIN)
00977               FC(IS)=FC(IS)+DTJALFA*C/TRA01(IS)
00978             ENDIF
00979           ENDDO
00980         ELSE
00981           DO IPOIN=1,NPOIN2
00982             IS=NPOIN3-NPOIN2+IPOIN
00983             IF(TRA01(IS).GT.EPS) THEN
00984               C=TRAIN*MAX(PLUIE(IPOIN),0.D0)-FC(IS)*PLUIE(IPOIN)
00985               FC(IS)=FC(IS)+DTJALFA*C/TRA01(IS)
00986             ENDIF
00987           ENDDO
00988         ENDIF
00989       ENDIF
00990 !
00991 !     FLUXES
00992 !
00993       IF(S0F%TYPR.NE.'0') THEN
00994         DO IPOIN=1,NPOIN3
00995           IF(OPTBAN.EQ.2) THEN
00996           IF(MASKPT(IPOIN).GT.0.5D0.AND.TRA01(IPOIN).GT.EPS) THEN
00997             FC(IPOIN)=FC(IPOIN)+DTJALFA*
00998      &                     (S0F%R(IPOIN)+TRA02(IPOIN))/TRA01(IPOIN)
00999           ENDIF
01000           ELSE
01001           IF(TRA01(IPOIN).GT.EPS) THEN
01002             FC(IPOIN)=FC(IPOIN)+DTJALFA*
01003      &                     (S0F%R(IPOIN)+TRA02(IPOIN))/TRA01(IPOIN)
01004           ENDIF
01005           ENDIF
01006         ENDDO
01007       ELSE
01008         IF(OPTBAN.EQ.2) THEN
01009         DO IPOIN=1,NPOIN3
01010           IF(MASKPT(IPOIN).GT.0.5D0.AND.TRA01(IPOIN).GT.EPS) THEN
01011             FC(IPOIN)=FC(IPOIN)+DTJALFA*TRA02(IPOIN)/TRA01(IPOIN)
01012           ENDIF
01013         ENDDO
01014         ELSE
01015         DO IPOIN=1,NPOIN3
01016           IF(TRA01(IPOIN).GT.EPS) THEN
01017             FC(IPOIN)=FC(IPOIN)+DTJALFA*TRA02(IPOIN)/TRA01(IPOIN)
01018           ENDIF
01019         ENDDO
01020         ENDIF
01021       ENDIF
01022 !
01023 !     DTJ WAS THE REMAINING TIME, ALFA*DTJ HAS BEEN DONE, THE REST IS:
01024       DTJ = DTJ * (1.D0-ALFA)
01025       NITER = NITER + 1
01026       TOTITER = TOTITER + 1
01027       IF(NITER.GE.100) THEN
01028         WRITE (LU,*) 'MURD3D: ITERATION NO. REACHED ',NITER,', STOP.'
01029         WRITE (LU,*) 'ALFA = ',ALFA
01030         WRITE (LU,*) 'GUILTY POINT = ',IGUILT
01031         CALL PLANTE(1)
01032         STOP
01033       ENDIF
01034 !
01035       IF(ALFA.LT.1.D0) GOTO 10
01036 !
01037 !-----------------------------------------------------------------------
01038 !
01039       IF(INFOR) THEN
01040         IF(LNG.EQ.1) WRITE(LU,101) SCHCF,NITER
01041         IF(LNG.EQ.2) WRITE(LU,102) SCHCF,NITER
01042       ENDIF
01043 !
01044 101   FORMAT(' MURD3D OPTION : ',1I2,'    ',1I4,' ITERATIONS')
01045 102   FORMAT(' MURD3D OPTION: ',1I2,'    ',1I4,' ITERATIONS')
01046 !
01047 !-----------------------------------------------------------------------
01048 !
01049       RETURN
01050       END

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