murd3d_pos.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac3d\murd3d_pos.f
00002 !
00080                      SUBROUTINE MURD3D_POS
00081 !                    *********************
00082 !
00083      &(FC,FN,VOLU,SVOLU,VOLUN,SVOLUN,VOLU2,SVOLU2,RMASS,
00084      & TRA01,TRA02,TRA03,STRA01,STRA02,STRA03,MESH2,MESH3,
00085      & NELEM3,NPOIN3,DT,SCHCF,MSK,MASKEL,INFOR,CALFLU,FLUX,FLUEXT,
00086      & S0F,NSCE,SOURCES,FSCE,RAIN,PLUIE,TRAIN,NPOIN2,
00087      & OPTBAN,FLODEL,FLOPAR,GLOSEG,DIMGLO,NSEG,NPLAN,
00088      & T5,FLUX_REMOVED,SAVED_VOLU2,SAVED_F,OPTION,IELM3,NITMAX)
00089 !
00090 !***********************************************************************
00091 ! TELEMAC3D   V6P2                                   21/08/2010
00092 !***********************************************************************
00093 !
00094 !
00095 !
00096 !
00097 !
00098 !
00099 !
00100 !
00101 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00102 !| CALFLU         |-->| INDICATE IF FLUX IS CALCULATED FOR BALANCE
00103 !| DIMGLO         |-->| FIRST DIMENSION OF ARRAY GLOSEG
00104 !| DT             |-->| TIME STEP
00105 !| FC             |<->| VARIABLE AFTER CONVECTION
00106 !| FLODEL         |-->| FLUXES BY SEGMENT
00107 !| FLOPAR         |-->| FLUXES BY SEGMENT, ASSEMBLED IN PARALLEL
00108 !| FLUEXT         |-->| OUTPUT FLUX BY NODE
00109 !| FLUX_REMOVED   |<->| TOTAL FLUX REMOVED OF EACH POINT
00110 !| FLUX           |<->| GLOBAL FLUXES TO BE CHANGED
00111 !| FN             |-->| VARIABLE AT TIME N
00112 !| FSCE           |-->| DIRICHLET BOUNDARY CONDITIONS OF F
00113 !| GLOSEG         |-->| FIRST AND SECOND POINT OF SEGMENTS
00114 !| IELM3          |-->| TYPE OF ELEMENT (41: PRISM, ETC.)
00115 !| INFOR          |-->| INFORMATIONS FOR SOLVERS
00116 !| MASKEL         |-->| MASKING OF ELEMENTS
00117 !|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
00118 !| MESH2          |<->| 2D MESH
00119 !| MESH3          |<->| 3D MESH
00120 !| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
00121 !| NELEM3         |-->| NUMBER OF ELEMENTS IN 3D
00122 !| NITMAX         |-->| MAXIMUM NUMBER OF ITERATIONS
00123 !| NPLAN          |-->| NUMBER OF PLANES IN THE 3D MESH OF PRISMS
00124 !| NPOIN2         |-->| NUMBER OF POINTS IN 2D
00125 !| NPOIN3         |-->| NUMBER OF 3D POINTS
00126 !| NSCE           |-->| NUMBER OF GIVEN POINTS FOR SOURCES
00127 !| NSEG           |-->| NUMBER OF SEGMENTS
00128 !| OPTBAN         |-->| OPTION FOR TIDAL FLATS, IF 1, FREE SURFACE
00129 !|                |   | MODIFIED AND PIECE-WISE LINEAR
00130 !| PLUIE          |-->| RAIN IN M/S MULTIPLIED BY VOLU2D
00131 !| RAIN           |-->| IF YES, THERE IS RAIN OR EVAPORATION
00132 !| RMASS          |<->| REMAINING MASSES
00133 !| S0F            |-->| EXPLICIT SOURCE TERM
00134 !| SAVED_F        |<->| TRACER SAVED
00135 !| SAVED_VOLU2    |<->| VOLUME VOLU2 SAVED
00136 !| SCHCF          |-->| ADVECTION SCHEME FOR F
00137 !| SOURCES        |-->| SOURCES
00138 !| STRA01         |<->| STRUCTURE OF TRA01
00139 !| STRA02         |<->| STRUCTURE OF TRA02
00140 !| STRA03         |<->| STRUCTURE OF TRA03
00141 !| SVOLU          |-->| STRUCTURE OF VOLU
00142 !| SVOLU2         |<->| STRUCTURE OF VOLU2
00143 !| SVOLUN         |-->| STRUCTURE OF VOLUN
00144 !| T5             |<->| WORK ARRAY
00145 !| TRA01          |<->| WORK ARRAY OF DIMENSION NPOIN3 EQUIVALENT TO
00146 !|                |   | VOLU2 FOR CURRENT FINAL TIME
00147 !| TRA02          |<->| WORK ARRAY OF DIMENSION NPOIN3
00148 !| TRA03          |<->| WORK ARRAY OF DIMENSION NPOIN3
00149 !| TRAIN          |-->| VALUE OF TRACER IN RAIN
00150 !| VOLU           |-->| CONTROL VOLUME AT TIME N+1
00151 !| VOLU2          |<->| LIKE VOLU, BUT ASSEMBLED IN PARALLEL
00152 !| VOLUN          |-->| CONTROL VOLUME AT TIME N
00153 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00154 !
00155       USE BIEF
00156       USE DECLARATIONS_TELEMAC
00157 !
00158       IMPLICIT NONE
00159       INTEGER LNG,LU
00160       COMMON/INFO/LNG,LU
00161 !
00162 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00163 !
00164       INTEGER, INTENT(IN)             :: SCHCF,NELEM3,NPOIN3,NPOIN2
00165       INTEGER, INTENT(IN)             :: NSCE,OPTBAN,NSEG,NPLAN,DIMGLO
00166       INTEGER, INTENT(IN)             :: GLOSEG(DIMGLO,2),OPTION,IELM3
00167       INTEGER, INTENT(IN)             :: NITMAX
00168 !
00169       DOUBLE PRECISION, INTENT(INOUT) :: FC(NPOIN3)
00170       DOUBLE PRECISION, INTENT(IN)    :: FN(NPOIN3),PLUIE(NPOIN2)
00171 !
00172       DOUBLE PRECISION, INTENT(IN)    :: MASKEL(NELEM3)
00173       DOUBLE PRECISION, INTENT(IN), TARGET    :: FLUEXT(NPOIN3)
00174 !
00175       DOUBLE PRECISION, INTENT(IN)    :: VOLUN(NPOIN3), VOLU(NPOIN3)
00176       DOUBLE PRECISION, INTENT(INOUT) :: VOLU2(NPOIN3)
00177       DOUBLE PRECISION, INTENT(INOUT) :: TRA01(NPOIN3),FLUX
00178       DOUBLE PRECISION, INTENT(INOUT) :: TRA02(NPOIN3)
00179       DOUBLE PRECISION, INTENT(INOUT), TARGET :: TRA03(NPOIN3)
00180       DOUBLE PRECISION, INTENT(IN)    :: DT,TRAIN,FSCE(NSCE)
00181 !
00182       TYPE(BIEF_OBJ), INTENT(INOUT)   :: SVOLU2
00183       TYPE(BIEF_OBJ), INTENT(INOUT)   :: FLUX_REMOVED
00184       TYPE(BIEF_OBJ), INTENT(INOUT)   :: SAVED_VOLU2,SAVED_F,T5
00185       TYPE(BIEF_OBJ), INTENT(IN)      :: SOURCES,S0F,SVOLU,SVOLUN
00186       TYPE(BIEF_OBJ), INTENT(INOUT)   :: STRA01,STRA02,STRA03
00187       TYPE(BIEF_MESH), INTENT(INOUT)  :: MESH2,MESH3
00188 !
00189 !     DIMENSION OF FLODEL AND FLOPAR=NSEG2D*NPLAN+NPOIN2*NETAGE
00190       DOUBLE PRECISION, INTENT(IN)    :: FLODEL(*),FLOPAR(*)
00191       DOUBLE PRECISION, INTENT(INOUT) :: RMASS(*)
00192 !                                        SIZE IN MEMORY = 30*NELEM
00193 !                                        THIS IS ENOUGH
00194 !
00195       LOGICAL, INTENT(IN)             :: MSK,INFOR,CALFLU,RAIN
00196 !
00197 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00198 !
00199       INTEGER IPOIN,NITER,IS,IIS,I,NSEGH,NSEGV,OPT,IR
00200       INTEGER I1,I2,IPLAN,ISEG3D,I2D,I3D,IPTFR,SIZEINDIC
00201       INTEGER REMAIN_SEG,NEWREMAIN,REMAIN_TOT
00202 !
00203 !-----------------------------------------------------------------------
00204 !
00205       DOUBLE PRECISION P_DSUM
00206       EXTERNAL         P_DSUM
00207       INTEGER  P_ISUM
00208       EXTERNAL P_ISUM
00209 !
00210       DOUBLE PRECISION RINIT,C,NEWVOL,RFLUX,RFLUX_OLD
00211       DOUBLE PRECISION VOLSEG1,VOLSEG2
00212 !
00213       DOUBLE PRECISION EPS
00214       DATA EPS /1.D-6/
00215       DOUBLE PRECISION ALLOW
00216       DATA ALLOW /1.D-5/
00217       DOUBLE PRECISION REDUC
00218       DATA REDUC /1.D-9/
00219       DOUBLE PRECISION EPS_VOLUME
00220       DATA EPS_VOLUME /1.D-8/
00221 !
00222       LOGICAL TESTING
00223       DATA    TESTING/.FALSE./
00224 !
00225 !-----------------------------------------------------------------------
00226 !
00227 !     INDIC WILL BE A LIST OF SEGMENTS WITH NON ZERO FLUXES
00228 !
00229       LOGICAL DEJA
00230       DATA DEJA/.FALSE./
00231       INTEGER, ALLOCATABLE :: INDIC(:)
00232       SAVE
00233 !
00234 !
00235       CALL CPSTVC(SVOLU2,STRA01)
00236       CALL CPSTVC(SVOLU2,STRA02)
00237       CALL CPSTVC(SVOLU2,STRA03)
00238 !
00239 !     TRA03 WILL BE THE SHARED ASSEMBLED FLUEXT IN PARALLEL
00240 !
00241       IF(NCSIZE.GT.1) THEN
00242         CALL OV('X=Y     ',TRA03,FLUEXT,FLUEXT,0.D0,NPOIN3)
00243 !       TRA03 WILL BE THE ASSEMBLED AND SHARED FLUEXT
00244         CALL PARCOM(STRA03,2,MESH3)
00245 !       SHARES AFTER SUMMING (AS WILL BEEN DONE WITH VOLUMES)
00246         DO IPTFR=1,NPTIR
00247           I2D=MESH2%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00248           DO IPLAN=1,NPLAN
00249             I3D=I2D+(IPLAN-1)*NPOIN2
00250             TRA03(I3D)=TRA03(I3D)*MESH3%FAC%R(I3D)
00251           ENDDO
00252         ENDDO
00253       ENDIF
00254 !
00255       NSEGH=NSEG*NPLAN
00256       NSEGV=NPOIN2*(NPLAN-1)
00257       IF(SCHCF.EQ.ADV_LPO_TF) THEN
00258 !       HORIZONTAL AND VERTICAL SEGMENTS
00259         REMAIN_SEG=NSEGH+NSEGV
00260         OPT=1
00261       ELSEIF(SCHCF.EQ.ADV_NSC_TF) THEN
00262 !       ALL SEGMENTS
00263         IF(IELM3.EQ.41) THEN
00264           REMAIN_SEG=NSEGH+NSEGV+2*NSEG*(NPLAN-1)
00265         ELSEIF(IELM3.EQ.51) THEN
00266           REMAIN_SEG=NSEGH+NSEGV+NSEG*(NPLAN-1)
00267         ELSE
00268           WRITE(LU,*) 'UNKNOWN ELEMENT IN MURD3D_POS:',IELM3
00269           CALL PLANTE(1)
00270           STOP
00271         ENDIF
00272         OPT=2
00273       ELSE
00274         WRITE(LU,*) 'UNKNOWN SCHEME IN MURD3D_POS:',SCHCF
00275         CALL PLANTE(1)
00276         STOP
00277       ENDIF
00278 !
00279       IF(.NOT.DEJA) THEN
00280         ALLOCATE(INDIC(REMAIN_SEG))
00281         SIZEINDIC=REMAIN_SEG
00282         DEJA=.TRUE.
00283       ELSE
00284         IF(REMAIN_SEG.GT.SIZEINDIC) THEN
00285 !         LARGER SIZE OF INDIC REQUIRED (CASE OF SEVERAL CALLS WITH
00286 !         DIFFERENT SCHEMES THAT IMPLY DIFFERENT NUMBERS OF SEGMENTS
00287 !         LIKE SCHEMES LPO_TF AND NSC_TF IN PRISMS
00288           DEALLOCATE(INDIC)
00289           ALLOCATE(INDIC(REMAIN_SEG))
00290           SIZEINDIC=REMAIN_SEG
00291         ENDIF
00292       ENDIF
00293 !
00294       IF(OPTION.EQ.2) CALL CPSTVC(SVOLU2,FLUX_REMOVED)
00295 !
00296 !***********************************************************************
00297 !
00298 !     COPIES FLUXES FROM FLODEL TO ARRAY RMASS (REMAINING MASSES
00299 !     TO BE TRANSFERRED AND THEIR ADDRESS IN INDIC)
00300 !
00301       IF(OPTION.EQ.1) THEN
00302         DO I=1,REMAIN_SEG
00303           INDIC(I)=I
00304           RMASS(I)=DT*FLOPAR(I)
00305         ENDDO
00306       ELSEIF(OPTION.EQ.2) THEN
00307         DO I=1,REMAIN_SEG
00308           INDIC(I)=I
00309           RMASS(I)=-DT*FLOPAR(I)
00310         ENDDO
00311       ENDIF
00312 !
00313 !     SHARES ASSEMBLED FLUXES ON INTERFACE SEGMENTS BY:
00314 !     DIVIDING BY 2 ON INTERFACE HORIZONTAL AND CROSSED SEGMENTS
00315 !     MULTIPLYING BY FAC ON VERTICAL FLUXES
00316 !     THIS WILL GIVE THE SAME UPWINDING INFORMATION
00317 !
00318       IF(NCSIZE.GT.1) THEN
00319         CALL SHARE_3D_FLUXES(RMASS,1.D0,NPLAN,MESH2,MESH3,OPT)
00320       ENDIF
00321 !
00322 !     REMAINING FLUXES (SPLIT INTO POSITIVE AND NEGATIVE
00323 !                       SO THAT THEY SUM CORRECTLY IN PARALLEL
00324 !                       ABSOLUTE VALUES WOULD NOT SUM CORRECTLY)
00325 !
00326       RFLUX_OLD=0.D0
00327       DO I=1,REMAIN_SEG
00328         RFLUX_OLD=RFLUX_OLD+ABS(RMASS(I))
00329       ENDDO
00330       IF(NCSIZE.GT.1) RFLUX_OLD=P_DSUM(RFLUX_OLD)
00331       RINIT=RFLUX_OLD
00332       IF(TESTING) WRITE(LU,*) 'SOMME INITIALE DES ABS(FLUX)=',RINIT/DT
00333 !
00334 !     INITIAL VALUE OF TRACER = FN
00335 !
00336       CALL OV ('X=Y     ',FC,FN,FN,C,NPOIN3)
00337 !
00338 !     VOLU2 WILL BE THE VOLUME CHANGING PROGRESSIVELY FROM VOLUN TO VOLU
00339 !
00340       CALL OS ('X=Y     ',X=SVOLU2,Y=SVOLUN)
00341 !
00342 !     TAKES INTO ACCOUNT ENTERING EXTERNAL FLUXES
00343 !     THIS IS DONE WITHOUT CHANGING THE TRACER
00344 !    (BUT THE MASS OF TRACER CHANGES)
00345 !
00346       IF(NCSIZE.GT.1) THEN
00347         DO I=1,NPOIN3
00348 !         FLUEXT SHARED IN PARALLEL (HENCE SAME SIGN)
00349           VOLU2(I)=VOLU2(I)-MIN(TRA03(I),0.D0)*DT
00350         ENDDO
00351       ELSE
00352         DO I=1,NPOIN3
00353 !         FLUEXT SHARED IN PARALLEL (HENCE SAME SIGN)
00354           VOLU2(I)=VOLU2(I)-MIN(FLUEXT(I),0.D0)*DT
00355         ENDDO
00356       ENDIF
00357 !
00358       IF(CALFLU) THEN
00359         IF(NCSIZE.GT.1) THEN
00360           DO IPOIN = 1,NPOIN3
00361             FLUX = FLUX + DT*FC(IPOIN)*MIN(TRA03(IPOIN),0.D0)
00362           ENDDO
00363         ELSE
00364           DO IPOIN = 1,NPOIN3
00365             FLUX = FLUX + DT*FC(IPOIN)*MIN(FLUEXT(IPOIN),0.D0)
00366           ENDDO
00367         ENDIF
00368 !       ENTERING SOURCES
00369         IF(NSCE.GT.0) THEN
00370           DO IS=1,NSCE
00371             IIS=IS
00372 !           HERE IN PARALLEL SOURCES WITHOUT PARCOM
00373 !           ARE STORED AT ADRESSES IS+NSCE (SEE SOURCES_SINKS.F)
00374             IF(NCSIZE.GT.1) IIS=IIS+NSCE
00375             DO IPOIN=1,NPOIN3
00376               IF(SOURCES%ADR(IS)%P%R(IPOIN).GT.0.D0) THEN
00377                 FLUX=FLUX-DT*FSCE(IS)*SOURCES%ADR(IIS)%P%R(IPOIN)
00378               ENDIF
00379             ENDDO
00380           ENDDO
00381         ENDIF
00382       ENDIF
00383 !
00384 !
00385 !-----------------------------------------------------------------------
00386 !
00387       IF(NCSIZE.GT.1) CALL PARCOM(SVOLU2,2,MESH3)
00388 !
00389 !     ENTERING SOURCES (WITH FC = FSCE)
00390 !     IT WILL ALWAYS GIVE POSITIVE VOLUMES
00391       IF(NSCE.GT.0) THEN
00392         DO IS=1,NSCE
00393           DO IPOIN=1,NPOIN3
00394 !           HERE VERSION OF SOURCES ASSEMBLED IN PARALLEL
00395             IF(SOURCES%ADR(IS)%P%R(IPOIN).GT.0.D0) THEN
00396               VOLU2(IPOIN)=VOLU2(IPOIN)+DT*SOURCES%ADR(IS)%P%R(IPOIN)
00397               FC(IPOIN)=FC(IPOIN)+DT*(FSCE(IS)-FC(IPOIN))
00398      &                       *SOURCES%ADR(IS)%P%R(IPOIN)/VOLU2(IPOIN)
00399             ENDIF
00400           ENDDO
00401         ENDDO
00402       ENDIF
00403 !
00404 !     RAIN-EVAPORATION (HERE ONLY RAIN, NOT EVAPORATION, HENCE
00405 !                       VALUE OF TRACER IN RAIN TAKEN INTO ACCOUNT)
00406 !
00407       IF(RAIN) THEN
00408         DO IPOIN=1,NPOIN2
00409           IF(PLUIE(IPOIN).GT.0.D0) THEN
00410             IS=NPOIN3-NPOIN2+IPOIN
00411 !           ASSEMBLED FORM OF PLUIE NEEDED HERE
00412             VOLU2(IS)=VOLU2(IS)+DT*PLUIE(IPOIN)
00413 !           DILUTION EFFECT FOR ALL TRACERS
00414             FC(IS)=FC(IS)+DT*(TRAIN-FC(IS))*PLUIE(IPOIN)/VOLU2(IS)
00415           ENDIF
00416         ENDDO
00417       ENDIF
00418 !
00419       NITER = 0
00420 !
00421 777   CONTINUE
00422 !
00423       NITER = NITER + 1
00424 !
00425 !
00426 !
00427 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00428 !     FOR DISTRIBUTING THE VOLUMES BETWEEN SEGMENTS
00429 !
00430       IF(OPTION.EQ.2) THEN
00431 !
00432 !       FLUX_REMOVED (T6)    : TOTAL FLUX REMOVED OF EACH POINT
00433 !       SAVED_VOLU2 (T8)     : VOLUME VOLU2 SAVED
00434 !       SAVED_F (T9)         : TRACER SAVED
00435 !
00436         IF(NITER.EQ.1) THEN
00437           DO I=1,NPOIN3
00438             FLUX_REMOVED%R(I)=0.D0
00439             SAVED_VOLU2%R(I)=VOLU2(I)
00440             SAVED_F%R(I)=FC(I)
00441             T5%R(I)=FC(I)*VOLU2(I)
00442           ENDDO
00443           IF(NCSIZE.GT.1) THEN
00444 !           SHARES AFTER SUMMING (AS HAS BEEN DONE WITH FLUXES)
00445             DO IPTFR=1,NPTIR
00446               I2D=MESH2%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00447               DO IPLAN=1,NPLAN
00448                 I3D=I2D+(IPLAN-1)*NPOIN2
00449                 VOLU2(I3D)=VOLU2(I3D)*MESH3%FAC%R(I3D)
00450                 T5%R(I3D) = T5%R(I3D)*MESH3%FAC%R(I3D)
00451               ENDDO
00452             ENDDO
00453           ENDIF
00454         ELSE
00455 !         NOT ALL THE POINTS NEED TO BE INITIALISED NOW
00456           DO IR=1,REMAIN_SEG
00457             I=INDIC(IR)
00458             I1=GLOSEG(I,1)
00459             I2=GLOSEG(I,2)
00460             FLUX_REMOVED%R(I1)=0.D0
00461             FLUX_REMOVED%R(I2)=0.D0
00462 !           SAVING THE DEPTH AND TRACER
00463             SAVED_VOLU2%R(I1)=VOLU2(I1)
00464             SAVED_VOLU2%R(I2)=VOLU2(I2)
00465             SAVED_F%R(I1)=FC(I1)
00466             SAVED_F%R(I2)=FC(I2)
00467             T5%R(I1)=FC(I1)*VOLU2(I1)
00468             T5%R(I2)=FC(I2)*VOLU2(I2)
00469           ENDDO
00470 !         CANCELLING INTERFACE POINTS (SOME MAY BE ISOLATED IN A SUBDOMAIN
00471 !         AT THE TIP OF AN ACTIVE SEGMENT WHICH IS ELSEWHERE)
00472           IF(NCSIZE.GT.1) THEN
00473             DO IPTFR=1,NPTIR
00474               I2D=MESH3%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00475               DO IPLAN=1,NPLAN
00476                 I3D=(IPLAN-1)*NPOIN2+I2D
00477                 FLUX_REMOVED%R(I3D)=0.D0
00478 !               SAVING THE VOLUME AND TRACER
00479                 SAVED_VOLU2%R(I3D)=VOLU2(I3D)
00480                 SAVED_F%R(I3D)=FC(I3D)
00481                 VOLU2(I3D)=VOLU2(I3D)*MESH3%FAC%R(I3D)
00482                 T5%R(I3D) = T5%R(I3D)*MESH3%FAC%R(I3D)
00483               ENDDO
00484             ENDDO
00485           ENDIF
00486         ENDIF
00487         DO I=1,REMAIN_SEG
00488           ISEG3D=INDIC(I)
00489           I1=GLOSEG(ISEG3D,1)
00490           I2=GLOSEG(ISEG3D,2)
00491 !         POSITIVE FLUXES FROM 1 TO 2 !!!
00492           IF(RMASS(ISEG3D).GT.EPS_VOLUME) THEN
00493             FLUX_REMOVED%R(I1)=FLUX_REMOVED%R(I1)+RMASS(ISEG3D)
00494             VOLU2(I1)=0.D0
00495             T5%R(I1)=0.D0
00496           ELSEIF(RMASS(ISEG3D).LT.-EPS_VOLUME) THEN
00497             FLUX_REMOVED%R(I2)=FLUX_REMOVED%R(I2)-RMASS(ISEG3D)
00498             VOLU2(I2)=0.D0
00499             T5%R(I2)=0.D0
00500           ENDIF
00501         ENDDO
00502 !
00503         IF(NCSIZE.GT.1) CALL PARCOM(FLUX_REMOVED,2,MESH3)
00504 !
00505 !       FOR ISOLATED POINTS CONNECTED TO AN ACTIVE SEGMENT
00506 !       THAT IS IN ANOTHER SUBDOMAIN
00507         IF(NCSIZE.GT.1) THEN
00508           DO IPTFR=1,NPTIR
00509             I2D=MESH3%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00510             DO IPLAN=1,NPLAN
00511               I3D=(IPLAN-1)*NPOIN2+I2D
00512               IF(FLUX_REMOVED%R(I3D).GT.EPS_VOLUME) THEN
00513 !               ALL VOLUME SHARED
00514                 VOLU2(I3D)=0.D0
00515                 T5%R(I3D)=0.D0
00516               ENDIF
00517             ENDDO
00518           ENDDO
00519         ENDIF
00520 !
00521       ELSEIF(OPTION.EQ.1) THEN
00522 !
00523         IF(NCSIZE.GT.1) THEN
00524 !         SHARES AFTER SUMMING (AS HAS BEEN DONE WITH FLUXES)
00525           DO IPTFR=1,NPTIR
00526             I2D=MESH2%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00527             DO IPLAN=1,NPLAN
00528               I3D=I2D+(IPLAN-1)*NPOIN2
00529               VOLU2(I3D)=VOLU2(I3D)*MESH3%FAC%R(I3D)
00530             ENDDO
00531           ENDDO
00532         ENDIF
00533 !
00534       ENDIF
00535 !
00536 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00537 !
00538 !
00539 !     FROM HERE RMASS IS THE REMAINING MASSES TO BE PASSED BETWEEN POINTS
00540 !
00541 !     COMPUTES THE NEW VOLUME WITH FV SCHEME + DB
00542 !
00543 !     HORIZONTAL AND VERTICAL FLUXES TREATED TOGETHER
00544 !
00545       RFLUX=0.D0
00546       NEWREMAIN=0
00547 !
00548       IF(OPTION.EQ.1) THEN
00549 !
00550       DO I=1,REMAIN_SEG
00551         ISEG3D=INDIC(I)
00552         IF(RMASS(ISEG3D).GT.EPS_VOLUME) THEN
00553 !         FLUX FROM 2 TO 1 !!! (SEE REMARKS AND HOW RMASS INITIALISED)
00554           I1=GLOSEG(ISEG3D,1)
00555           I2=GLOSEG(ISEG3D,2)
00556           IF(RMASS(ISEG3D).GT.VOLU2(I2)) THEN
00557             RMASS(ISEG3D)=RMASS(ISEG3D)-VOLU2(I2)
00558             NEWVOL=VOLU2(I1)+VOLU2(I2)
00559             IF(NEWVOL.GT.0.D0) THEN
00560               FC(I1)=FC(I1)+VOLU2(I2)*(FC(I2)-FC(I1))/NEWVOL
00561               VOLU2(I1)=NEWVOL
00562               VOLU2(I2)=0.D0
00563             ENDIF
00564             RFLUX=RFLUX+RMASS(ISEG3D)
00565             NEWREMAIN=NEWREMAIN+1
00566             INDIC(NEWREMAIN)=ISEG3D
00567           ELSE
00568             NEWVOL=VOLU2(I1)+RMASS(ISEG3D)
00569             IF(NEWVOL.GT.0.D0) THEN
00570               FC(I1)=FC(I1)+RMASS(ISEG3D)*(FC(I2)-FC(I1))/NEWVOL
00571             ENDIF
00572             VOLU2(I1)=NEWVOL
00573             VOLU2(I2)=VOLU2(I2)-RMASS(ISEG3D)
00574           ENDIF
00575         ELSEIF(RMASS(ISEG3D).LT.-EPS_VOLUME) THEN
00576 !         FLUX FROM 1 TO 2 !!! (SEE REMARKS AND HOW RMASS INITIALISED)
00577           I1=GLOSEG(ISEG3D,1)
00578           I2=GLOSEG(ISEG3D,2)
00579           IF(-RMASS(ISEG3D).GT.VOLU2(I1)) THEN
00580             RMASS(ISEG3D)=RMASS(ISEG3D)+VOLU2(I1)
00581             NEWVOL=VOLU2(I2)+VOLU2(I1)
00582             IF(NEWVOL.GT.0.D0) THEN
00583               FC(I2)=FC(I2)+VOLU2(I1)*(FC(I1)-FC(I2))/NEWVOL
00584               VOLU2(I2)=NEWVOL
00585               VOLU2(I1)=0.D0
00586             ENDIF
00587             RFLUX=RFLUX-RMASS(ISEG3D)
00588             NEWREMAIN=NEWREMAIN+1
00589             INDIC(NEWREMAIN)=ISEG3D
00590           ELSE
00591             NEWVOL=VOLU2(I2)-RMASS(ISEG3D)
00592             IF(NEWVOL.GT.0.D0) THEN
00593               FC(I2)=FC(I2)-RMASS(ISEG3D)*(FC(I1)-FC(I2))/NEWVOL
00594             ENDIF
00595             VOLU2(I2)=NEWVOL
00596             VOLU2(I1)=VOLU2(I1)+RMASS(ISEG3D)
00597           ENDIF
00598         ENDIF
00599       ENDDO
00600 !
00601       ELSEIF(OPTION.EQ.2) THEN
00602 !
00603       DO IR=1,REMAIN_SEG
00604         I=INDIC(IR)
00605         IF(RMASS(I).GT.EPS_VOLUME) THEN
00606           I1=GLOSEG(I,1)
00607 !         FLUX FROM 1 TO 2 !!! (SEE REMARKS AND HOW RMASS INITIALISED)
00608           IF(SAVED_VOLU2%R(I1).GT.0.D0) THEN
00609             I2=GLOSEG(I,2)
00610 !           SHARING ON DEMAND: RMASS(I)/FLUX_REMOVED%R(I1) IS A PERCENTAGE
00611             VOLSEG1=SAVED_VOLU2%R(I1)*RMASS(I)/FLUX_REMOVED%R(I1)
00612 !           END OF SHARING ON DEMAND
00613             IF(RMASS(I).GE.VOLSEG1) THEN
00614 !             ALL VOLSEG1 WILL BE TRANSFERED TO POINT2
00615 !             VOLSEG1 > 0, HENCE VOLU2(I2) ALSO
00616               RMASS(I) =RMASS(I) -VOLSEG1
00617               VOLU2(I2)=VOLU2(I2)+VOLSEG1
00618 !             GROUPING H*F
00619               T5%R(I2)=T5%R(I2)+VOLSEG1*SAVED_F%R(I1)
00620 !             THIS MAY BE DONE SEVERAL TIMES FOR THE SAME POINT
00621 !             BUT THE LAST ONE WILL BE THE GOOD ONE
00622               FC(I2)=T5%R(I2)/VOLU2(I2)
00623               RFLUX=RFLUX+RMASS(I)
00624               NEWREMAIN=NEWREMAIN+1
00625               INDIC(NEWREMAIN)=I
00626             ELSE
00627               VOLSEG1=VOLSEG1-RMASS(I)
00628 !             GATHERING VOLUMES (HERE VOLU2(I2) WILL REMAIN POSITIVE)
00629               VOLU2(I2)=VOLU2(I2)+RMASS(I)
00630               VOLU2(I1)=VOLU2(I1)+VOLSEG1
00631               T5%R(I1)=T5%R(I1)+VOLSEG1*SAVED_F%R(I1)
00632               T5%R(I2)=T5%R(I2)+RMASS(I)*SAVED_F%R(I1)
00633               FC(I1)=T5%R(I1)/VOLU2(I1)
00634               FC(I2)=T5%R(I2)/VOLU2(I2)
00635             ENDIF
00636           ELSE
00637             RFLUX=RFLUX+RMASS(I)
00638             NEWREMAIN=NEWREMAIN+1
00639             INDIC(NEWREMAIN)=I
00640           ENDIF
00641         ELSEIF(RMASS(I).LT.-EPS_VOLUME) THEN
00642           I2=GLOSEG(I,2)
00643 !         FLUX FROM 2 TO 1 !!! (SEE REMARKS AND HOW RMASS INITIALISED)
00644           IF(SAVED_VOLU2%R(I2).GT.0.D0) THEN
00645             I1=GLOSEG(I,1)
00646 !           SHARING ON DEMAND
00647             VOLSEG2=-SAVED_VOLU2%R(I2)*RMASS(I)/FLUX_REMOVED%R(I2)
00648 !           END OF SHARING ON DEMAND
00649             IF(-RMASS(I).GE.VOLSEG2) THEN
00650 !             ALL VOLSEG2 WILL BE TRANSFERED TO POINT 1
00651 !             VOLSEG2 > 0, HENCE VOLU2(I1) ALSO
00652               VOLU2(I1)=VOLU2(I1)+VOLSEG2
00653               RMASS(I) =RMASS(I) +VOLSEG2
00654               T5%R(I1)=T5%R(I1)+VOLSEG2*SAVED_F%R(I2)
00655               FC(I1)=T5%R(I1)/VOLU2(I1)
00656               RFLUX=RFLUX-RMASS(I)
00657               NEWREMAIN=NEWREMAIN+1
00658               INDIC(NEWREMAIN)=I
00659             ELSE
00660               VOLSEG2=VOLSEG2+RMASS(I)
00661 !             GATHERING VOLUMES (HERE VOLU2(I1) WILL REMAIN POSITIVE)
00662               VOLU2(I1)=VOLU2(I1)-RMASS(I)
00663               VOLU2(I2)=VOLU2(I2)+VOLSEG2
00664               T5%R(I1)=T5%R(I1)-RMASS(I)*SAVED_F%R(I2)
00665               T5%R(I2)=T5%R(I2)+VOLSEG2*SAVED_F%R(I2)
00666               FC(I1)=T5%R(I1)/VOLU2(I1)
00667               FC(I2)=T5%R(I2)/VOLU2(I2)
00668             ENDIF
00669           ELSE
00670             RFLUX=RFLUX-RMASS(I)
00671             NEWREMAIN=NEWREMAIN+1
00672             INDIC(NEWREMAIN)=I
00673           ENDIF
00674         ENDIF
00675       ENDDO
00676 !
00677 !     ELSE
00678 !       UNKNOWN OPTION
00679       ENDIF
00680 !
00681       REMAIN_SEG=NEWREMAIN
00682 !
00683 !     MERGES VOLUMES AND FC AT INTERFACE POINTS
00684 !
00685       IF(NCSIZE.GT.1) THEN
00686         DO IPTFR=1,NPTIR
00687           I2D=MESH2%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00688           DO IPLAN=1,NPLAN
00689             I3D=I2D+(IPLAN-1)*NPOIN2
00690 !           ARRAY WITH VOLUME*FC AT INTERFACE POINTS
00691             TRA01(I3D)=VOLU2(I3D)*FC(I3D)
00692           ENDDO
00693         ENDDO
00694 !       SUMS VOLUME*FC AT INTERFACE POINTS
00695         CALL PARCOM(STRA01,2,MESH3)
00696 !       SUMS THE NEW POSITIVE PARTIAL VOLUMES OF INTERFACE POINTS
00697         CALL PARCOM(SVOLU2,2,MESH3)
00698         DO IPTFR=1,NPTIR
00699           I2D=MESH2%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00700           DO IPLAN=1,NPLAN
00701             I3D=I2D+(IPLAN-1)*NPOIN2
00702 !           ARRAY WITH VOLUME*F AT INTERFACE POINTS
00703             IF(VOLU2(I3D).GT.0.D0) THEN
00704               FC(I3D)=TRA01(I3D)/VOLU2(I3D)
00705             ENDIF
00706           ENDDO
00707         ENDDO
00708       ENDIF
00709 !
00710       REMAIN_TOT=REMAIN_SEG
00711       IF(NCSIZE.GT.1) THEN
00712         RFLUX=P_DSUM(RFLUX)
00713 !       WILL NOT SUM CORRECTLY IN PARALLEL, BUT ONLY TEST IF .EQ.0
00714         REMAIN_TOT=P_ISUM(REMAIN_TOT)
00715       ENDIF
00716 !
00717 !     4 POSSIBLE REASONS FOR STOPPING:
00718 !
00719 !     1) THERE IS NO REMAINING FLUX
00720 !     2) REMAINING FLUXES DO NOT CHANGE (POSITIVE AND NEGATIVE)
00721 !        WITH SOME ABSOLUTE ALLOWANCE OR REDUCTION COEFFICIENT
00722 !     3) ALL SEGMENTS HAVE BEEN TREATED
00723 !     4) MAXIMUM NUMBER OF ITERATIONS IS REACHED
00724 !
00725       IF( RFLUX.NE.0.D0                                   .AND.
00726      &    ABS(RFLUX-RFLUX_OLD).GT.MIN(RINIT*REDUC,ALLOW)  .AND.
00727      &               REMAIN_TOT.NE.0                      .AND.
00728      &               NITER.LT.NITMAX                   ) THEN
00729         RFLUX_OLD=RFLUX
00730         GO TO 777
00731       ENDIF
00732 !
00733 !     TAKES INTO ACCOUNT EXITING EXTERNAL FLUXES
00734 !     THIS IS DONE WITHOUT CHANGING THE TRACER
00735 !    (BUT THE MASS OF TRACER CHANGES)
00736 !
00737       IF(CALFLU) THEN
00738 !       EXITING FLUXES
00739         IF(NCSIZE.GT.1) THEN
00740           DO IPOIN = 1,NPOIN3
00741             FLUX = FLUX + DT*FC(IPOIN)*MAX(TRA03(IPOIN),0.D0)
00742           ENDDO
00743         ELSE
00744           DO IPOIN = 1,NPOIN3
00745             FLUX = FLUX + DT*FC(IPOIN)*MAX(FLUEXT(IPOIN),0.D0)
00746           ENDDO
00747         ENDIF
00748 !       EXITING SOURCES
00749         IF(NSCE.GT.0) THEN
00750           DO IS=1,NSCE
00751             IIS=IS
00752 !           HERE IN PARALLEL SOURCES WITHOUT PARCOM
00753 !           ARE STORED AT ADRESSES IS+NSCE (SEE SOURCES_SINKS.F)
00754             IF(NCSIZE.GT.1) IIS=IIS+NSCE
00755             DO IPOIN=1,NPOIN3
00756               IF(SOURCES%ADR(IS)%P%R(IPOIN).LT.0.D0) THEN
00757                 FLUX=FLUX
00758      &              -DT*FC(IPOIN)*SOURCES%ADR(IIS)%P%R(IPOIN)
00759               ENDIF
00760             ENDDO
00761           ENDDO
00762         ENDIF
00763       ENDIF
00764 !
00765 !     RAIN-EVAPORATION (HERE ONLY EVAPORATION, NOT RAIN, HENCE
00766 !                       VALUE OF TRACER IN RAIN NOT TAKEN INTO ACCOUNT
00767 !                       AND ASSUMED TO BE 0)
00768 !
00769       IF(RAIN) THEN
00770         DO IPOIN=1,NPOIN2
00771           IF(PLUIE(IPOIN).LT.0.D0) THEN
00772             IS=NPOIN3-NPOIN2+IPOIN
00773             VOLU2(IS)=VOLU2(IS)+DT*PLUIE(IPOIN)
00774 !           DIVISION BY 0 NOT CHECKED (BUT COULD BE A PROBLEM AS PLUIE<0)
00775 !           CONCENTRATION EFFECT FOR ALL TRACERS (WRONG FOR TEMPERATURE ??)
00776             FC(IS)=FC(IS)-DT*FC(IS)*PLUIE(IPOIN)/VOLU2(IS)
00777           ENDIF
00778         ENDDO
00779       ENDIF
00780 !
00781       IF(TESTING) THEN
00782 !       EXITING SOURCES (WITH FC UNCHANGED)
00783         IF(NSCE.GT.0) THEN
00784           DO IS=1,NSCE
00785             DO IPOIN=1,NPOIN3
00786 !             HERE VERSION OF SOURCES ASSEMBLED IN PARALLEL
00787               IF(SOURCES%ADR(IS)%P%R(IPOIN).LT.0.D0) THEN
00788                 VOLU2(IPOIN)=VOLU2(IPOIN)+
00789      &                      DT*SOURCES%ADR(IS)%P%R(IPOIN)
00790               ENDIF
00791             ENDDO
00792           ENDDO
00793         ENDIF
00794 !       EXITING FLUXES
00795         IF(NCSIZE.GT.1) THEN
00796           DO I=1,NPOIN3
00797             TRA02(I)=MAX(TRA03(I),0.D0)
00798 !           SHARED FLUEXT IN TRA03 ERASED NOW (NOT USED AFTER)
00799 !           NOW VOLU IS COPIED INTO TRA03 FOR PARCOM
00800             TRA03(I)=VOLU(I)
00801           ENDDO
00802         ELSE
00803           DO I=1,NPOIN3
00804             TRA02(I)=MAX(FLUEXT(I),0.D0)
00805             TRA03(I)=VOLU(I)
00806           ENDDO
00807         ENDIF
00808         IF(NCSIZE.GT.1) THEN
00809           CALL PARCOM(STRA02,2,MESH3)
00810           CALL PARCOM(STRA03,2,MESH3)
00811         ENDIF
00812 !!!!!
00813 !!!!!   IMPORTANT WARNING: THIS SHOULD BE DONE EVEN WITHOUT TESTING
00814 !!!!!   IF WE WANT THE CORRECT VOLU2, BUT VOLU2 HERE IS NOT USED
00815 !!!!!   AFTER.
00816 !!!!!
00817         DO I=1,NPOIN3
00818           VOLU2(I)=VOLU2(I)-TRA02(I)*DT
00819           IF(VOLU2(I).LT.0.D0) THEN
00820             WRITE(LU,*) 'POINT ',I,' VOLU2=',VOLU2(I)
00821             CALL PLANTE(1)
00822             STOP
00823           ENDIF
00824         ENDDO
00825 !!!!!
00826 !!!!!   END OF IMPORTANT NOTE
00827 !!!!!
00828 !       CHECKS EQUALITY OF ASSEMBLED VOLUMES
00829         C=0.D0
00830         IF(NCSIZE.GT.1) THEN
00831           DO I=1,NPOIN3
00832             C=C+ABS(VOLU2(I)-TRA03(I))*MESH3%FAC%R(I)
00833           ENDDO
00834           C=P_DSUM(C)
00835         ELSE
00836           DO I=1,NPOIN3
00837 !                            VOLU
00838             C=C+ABS(VOLU2(I)-TRA03(I))
00839           ENDDO
00840         ENDIF
00841         WRITE(LU,*) 'ERROR ON VOLUMES = ',C
00842       ENDIF
00843 !
00844 !-----------------------------------------------------------------------
00845 !
00846 !     EXPLICIT SOURCE TERMS
00847 !
00848       IF(S0F%TYPR.NE.'0') THEN
00849         DO IPOIN=1,NPOIN3
00850           IF(VOLU2(IPOIN).GT.EPS) THEN
00851             FC(IPOIN)=FC(IPOIN)+DT*S0F%R(IPOIN)/VOLU2(IPOIN)
00852           ENDIF
00853         ENDDO
00854       ENDIF
00855 !
00856 !-----------------------------------------------------------------------
00857 !
00858       IF(INFOR) THEN
00859         IF(LNG.EQ.1) WRITE(LU,101) SCHCF,NITER
00860         IF(LNG.EQ.2) WRITE(LU,102) SCHCF,NITER
00861       ENDIF
00862 !
00863 101   FORMAT(1X,'MURD3D_POS OPTION: ',1I4,'  ',1I4,' ITERATIONS')
00864 102   FORMAT(1X,'MURD3D_POS OPTION: ',1I4,'  ',1I4,' ITERATIONS')
00865 !
00866 !-----------------------------------------------------------------------
00867 !
00868       RETURN
00869       END

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