cvtrvf_pos.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\cvtrvf_pos.f
00002 !
00097                      SUBROUTINE CVTRVF_POS
00098 !                    *********************
00099 !
00100      &(F,FN,FSCEXP,DIFT,CONV,H,HN,HPROP,UDEL,VDEL,DM1,ZCONV,SOLSYS,
00101      & VISC,VISC_S,SM,SMH,YASMH,SMI,YASMI,FBOR,MASKTR,MESH,
00102      & T1,T2,T3,T4,T5,T6,T7,T8,HNT,HT,AGGLOH,TE1,DT,ENTET,BILAN,
00103      & OPDTRA,MSK,MASKEL,S,MASSOU,OPTSOU,LIMTRA,KDIR,KDDL,NPTFR,FLBOR,
00104      & YAFLBOR,V2DPAR,UNSV2D,IOPT,FLBORTRA,MASKPT,GLOSEG1,GLOSEG2,NBOR,
00105      & OPTION,FLULIM,YAFLULIM,RAIN,PLUIE,TRAIN,GIVEN_FLUX,FLUX_GIVEN,
00106      & NITMAX)
00107 !
00108 !***********************************************************************
00109 ! BIEF   V7P0                                     13/06/2014
00110 !***********************************************************************
00111 !
00112 !
00113 !
00114 !
00115 !
00116 !
00117 !
00118 !
00119 !
00120 !
00121 !
00122 !
00123 !
00124 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00125 !| AGGLOH         |-->| MASS-LUMPING UTILISE DANS L'EQUATION DE CONTINUITE
00126 !| BILAN          |-->| LOGICAL TRIGGERING A MASS BALANCE INFORMATION
00127 !| CONV           |-->| LOGICAL, IF YES THERE IS ADVECTION OF F
00128 !| DIFT           |-->| LOGICAL, IF YES THERE IS DIFFUSION OF F
00129 !| DM1            |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
00130 !|                |   | IS DM1*GRAD(ZCONV), SEE SOLSYS.
00131 !| DT             |-->| TIME STEP
00132 !| ENTET          |-->| LOGICAL, IF YES INFORMATION IS GIVEN ON MASS
00133 !|                |   | CONSERVATION.
00134 !| F              |<--| F AT TIME T(N+1)
00135 !| FBOR           |-->| DIRICHLET CONDITIONS ON F.
00136 !| FLBOR          |-->| FLUXES AT BOUNDARIES
00137 !| FLBORTRA       |<->| TRACER FLUXES AT BOUNDARIES
00138 !| FLUX_GIVEN     |-->| IF GIVEN_FLUX=YES, THE FLUX IS GIVEN IN
00139 !|                |   | GIVEN_FLUX
00140 !| FN             |-->| F AT TIME T(N)
00141 !| FSCEXP         |-->| EXPLICIT PART OF THE SOURCE TERM
00142 !|                |   | EQUAL TO ZERO EVERYWHERE BUT ON SOURCES
00143 !|                |   | WHERE THERE IS FSCE - (1-TETAT) FN
00144 !|                |   | SEE DIFSOU
00145 !| GIVEN_FLUX     |-->| IF GIVEN_FLUX=YES, THE FLUX IS GIVEN IN
00146 !|                |   | GIVEN_FLUX AND WILL NOT BE COMPUTED HERE
00147 !| GLOSEG1        |-->| FIRST POINT OF SEGMENTS
00148 !| GLOSEG2        |-->| SECOND POINT OF SEGMENTS
00149 !| HNT,HT         |<--| WORK ARRAYS (MODIFIED DEPTHS TO TAKE MASS-LUMPING
00150 !|                |   | INTO ACCOUNT)
00151 !| HPROP          |-->| PROPAGATION DEPTH (DONE IN CVDFTR).
00152 !| IOPT           |-->| OPTIONS FOR COMPUTATION (NUMBER BETWEEN 0 AND 13)
00153 !|                |   | THE TENS (IOPT2, I.E. 0 OR 1):
00154 !|                |   | 0: UCONV OBEYS THE CONTINUITY EQUATION
00155 !|                |   | 1: UCONV DOES NOT OBEY THE CONTINUITY EQUATION
00156 !|                |   | THE UNITS (IOPT1, I.E. 0 TO 3): VARIANT FOR FLUXES
00157 !|                |   | 0: CONSTANT PER ELEMENT = 0
00158 !|                |   | 1: CHI-TUAN PHAM'S CONSTANT
00159 !|                |   | 2: N SCHEME
00160 !|                |   | 3: PSI SCHEME
00161 !| KDDL           |-->| CONVENTION FOR DEGREE OF FREEDOM
00162 !| KDIR           |-->| CONVENTION FOR DIRICHLET POINT
00163 !| LIMTRA         |-->| BOUNDARY CONDITIONS ON BOOUNDARY POINTS
00164 !| MASKEL         |-->| MASKING OF ELEMENTS
00165 !|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
00166 !| MASKPT         |-->| MASKING PER POINT.
00167 !| MASSOU         |-->| MASS OF TRACER ADDED BY SOURCE TERM
00168 !|                |   | SEE DIFSOU
00169 !| MESH           |-->| MESH STRUCTURE
00170 !| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
00171 !| NBOR           |-->| GLOBAL NUMBERS OF BOUNDARY POINTS
00172 !| NITMAX         |-->| MAXIMUM NUMBER OF ITERATIONS
00173 !| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
00174 !| OPDTRA         |-->| OPTION FOR THE DIFFUSION OF TRACERS
00175 !| OPTION         |-->| OPTION OF ALGORITHM FOR EDGE-BASED ADVECTION
00176 !|                |   | 1: FAST BUT SENSITIVE TO SEGMENT NUMBERING
00177 !|                |   | 2: INDEPENDENT OF SEGMENT NUMBERING
00178 !| OPTSOU         |-->| TYPE OF SOURCES
00179 !|                |   | 1: NORMAL
00180 !|                |   | 2: DIRAC
00181 !| PLUIE          |-->| RAIN OR EVAPORATION, IN M/S
00182 !| RAIN           |-->| IF YES: RAIN OR EVAPORATION
00183 !| S              |-->| VOID STRUCTURE
00184 !| SM             |-->| SOURCE TERMS.
00185 !| SMH            |-->| SOURCE TERM IN CONTINUITY EQUATION
00186 !| SMI            |-->| IMPLICIT SOURCE TERM
00187 !| SOLSYS         |-->| 1 OR 2. IF 2 ADVECTION FIELD IS UCONV + DM1*GRAD(ZCONV)
00188 !| T1             |<->| WORK BIEF_OBJ STRUCTURE
00189 !| T2             |<->| WORK BIEF_OBJ STRUCTURE
00190 !| T3             |<->| WORK BIEF_OBJ STRUCTURE
00191 !| T4             |<->| WORK BIEF_OBJ STRUCTURE
00192 !| T5             |<->| WORK BIEF_OBJ STRUCTURE
00193 !| T6             |<->| WORK BIEF_OBJ STRUCTURE
00194 !| T7             |<->| WORK BIEF_OBJ STRUCTURE
00195 !| T8             |<->| WORK BIEF_OBJ STRUCTURE
00196 !| TE1            |<->| WORK BIEF_OBJ STRUCTURE FOR ELEMENTS
00197 !| TRAIN          |-->| VALUE OF TRACER IN THE RAIN
00198 !| UDEL           |-->| X-COMPONENT OF ADVECTION VELOCITY
00199 !| UNSV2D         |-->| INVERSE OF V2DPAR
00200 !| V2DPAR         |-->| INTEGRAL OF 2D TEST FUNCTIONS, ASSEMBLED
00201 !|                |   | IN PARALLEL.
00202 !| VDEL           |-->| X-COMPONENT OF ADVECTION VELOCITY
00203 !| VISC           |-->| VISCOSITY COEFFICIENTS ALONG X,Y AND Z .
00204 !|                |   | IF P0 : PER ELEMENT
00205 !|                |   | IF P1 : PERR POINT
00206 !| VISC_S         |<->| WORK ARRAY FOR SAVING VISC
00207 !| YAFLBOR        |-->| IF YES FLBOR IS GIVEN
00208 !| YASMH          |-->| IF YES, SMH MUST BE TAKEN INTO ACCOUNT
00209 !| YASMI          |-->| IF YES, SMI MUST BE TAKEN INTO ACCOUNT
00210 !| ZCONV          |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
00211 !|                |   | IS DM1*GRAD(ZCONV), SEE SOLSYS.
00212 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00213 !
00214       USE BIEF, EX_CVTRVF_POS => CVTRVF_POS
00215 !
00216       IMPLICIT NONE
00217       INTEGER LNG,LU
00218       COMMON/INFO/LNG,LU
00219 !
00220 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00221 !
00222       INTEGER, INTENT(IN)             :: OPDTRA,OPTSOU,KDIR,NPTFR,SOLSYS
00223       INTEGER, INTENT(IN)             :: KDDL,IOPT,OPTION,NITMAX
00224       INTEGER, INTENT(IN)             :: GLOSEG1(*),GLOSEG2(*)
00225       INTEGER, INTENT(IN)             :: NBOR(NPTFR)
00226       INTEGER, INTENT(INOUT)          :: LIMTRA(NPTFR)
00227 !                                                               NSEG
00228       DOUBLE PRECISION, INTENT(IN)    :: DT,AGGLOH,TRAIN,FLULIM(*)
00229       DOUBLE PRECISION, INTENT(INOUT) :: MASSOU
00230       LOGICAL, INTENT(IN)             :: BILAN,CONV,YASMH,YAFLBOR,RAIN
00231       LOGICAL, INTENT(IN)             :: DIFT,MSK,ENTET,YASMI,YAFLULIM
00232       LOGICAL, INTENT(IN)             :: FLUX_GIVEN
00233       TYPE(BIEF_OBJ), INTENT(IN)      :: MASKEL,H,HN,DM1,ZCONV,MASKPT
00234       TYPE(BIEF_OBJ), INTENT(IN)      :: V2DPAR,UNSV2D,HPROP
00235       TYPE(BIEF_OBJ), INTENT(INOUT)   :: F,SM,HNT,HT
00236       TYPE(BIEF_OBJ), INTENT(IN)      :: UDEL,VDEL,FN,SMI,SMH
00237       TYPE(BIEF_OBJ), INTENT(INOUT)   :: TE1,FLBORTRA,FBOR
00238       TYPE(BIEF_OBJ), INTENT(INOUT)   :: T1,T2,T3,T4,T5,T6,T7,T8
00239       TYPE(BIEF_OBJ), INTENT(IN)      :: FSCEXP,S,MASKTR
00240       TYPE(BIEF_OBJ), INTENT(INOUT)   :: FLBOR
00241       TYPE(BIEF_OBJ), INTENT(IN)      :: VISC_S,VISC,PLUIE,GIVEN_FLUX
00242       TYPE(BIEF_MESH)                 :: MESH
00243 !
00244 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00245 !
00246       DOUBLE PRECISION P_DSUM,P_DMIN,P_DMAX
00247       EXTERNAL         P_DSUM,P_DMIN,P_DMAX
00248 !
00249       INTEGER I,IOPT1,IOPT2,NPOIN,IPTFR,I1,I2,NITER,REMAIN_SEG,NEWREMAIN
00250       INTEGER IR,N
00251 !
00252 !-----------------------------------------------------------------------
00253 !
00254       DOUBLE PRECISION C,CPREV,CINIT,HFL1,HFL2,TET,H1N,H2N,HSEG1,HSEG2
00255       CHARACTER(LEN=16) FORMUL
00256       DOUBLE PRECISION, POINTER, DIMENSION(:) :: FXMAT
00257       LOGICAL TESTING
00258       DATA TESTING/.FALSE./
00259       DOUBLE PRECISION EPS_FLUX
00260       DATA             EPS_FLUX/1.D-15/
00261 !
00262 !-----------------------------------------------------------------------
00263 !
00264 !     INDIC WILL BE A LIST OF SEGMENTS WITH NON ZERO FLUXES
00265 !
00266       LOGICAL DEJA
00267       DATA DEJA/.FALSE./
00268       INTEGER, ALLOCATABLE          :: INDIC(:)
00269       SAVE
00270       IF(.NOT.DEJA) THEN
00271         ALLOCATE(INDIC(MESH%NSEG))
00272         DEJA=.TRUE.
00273       ENDIF
00274 !
00275 !-----------------------------------------------------------------------
00276 !
00277       FXMAT=>MESH%MSEG%X%R(1:MESH%NSEG)
00278 !
00279 !-----------------------------------------------------------------------
00280 !
00281       NPOIN=H%DIM1
00282 !
00283 !     EXTRACTING OPTIONS, AND CONTROL
00284 !
00285       IOPT2=IOPT/10
00286       IOPT1=IOPT-10*IOPT2
00287       IF(IOPT1.LT.0.OR.IOPT1.GT.3) THEN
00288         IF(LNG.EQ.1) THEN
00289           WRITE(LU,*) 'CVTRVF_POS : OPTION IOPT1 INCONNUE : ',IOPT1
00290         ENDIF
00291         IF(LNG.EQ.2) THEN
00292           WRITE(LU,*) 'CVTRVF_POS: OPTION IOPT1 UNKNOWN: ',IOPT1
00293         ENDIF
00294         CALL PLANTE(1)
00295         STOP
00296       ENDIF
00297       IF(IOPT2.NE.0.AND.IOPT2.NE.1) THEN
00298         IF(LNG.EQ.1) THEN
00299           WRITE(LU,*) 'CVTRVF_POS : OPTION IOPT2 INCONNUE : ',IOPT2
00300         ENDIF
00301         IF(LNG.EQ.2) THEN
00302           WRITE(LU,*) 'CVTRVF_POS: OPTION IOPT2 UNKNOWN: ',IOPT2
00303         ENDIF
00304         CALL PLANTE(1)
00305         STOP
00306       ENDIF
00307 !
00308 !-----------------------------------------------------------------------
00309 !
00310 !     STARTING AGAIN FROM NON CORRECTED DEPTH
00311 !
00312       IF(TESTING) THEN
00313         C=1.D99
00314         CINIT=1.D99
00315         DO I=1,NPOIN
00316           C    =MIN(C    ,H%R(I))
00317           CINIT=MIN(CINIT,HN%R(I))
00318         ENDDO
00319         IF(NCSIZE.GT.1) THEN
00320           C=P_DMIN(C)
00321           CINIT=P_DMIN(CINIT)
00322         ENDIF
00323         WRITE(LU,*) 'AVANT TRAITEMENT HAUTEURS NEGATIVES, H MIN=',C
00324         WRITE(LU,*) 'AVANT TRAITEMENT HAUTEURS NEGATIVES, HN MIN=',CINIT
00325       ENDIF
00326 !
00327 !     CALCUL DES FLUX PAR NOEUDS
00328 !
00329       IF(.NOT.FLUX_GIVEN) THEN
00330 !
00331       FORMUL='HUGRADP         '
00332       IF(SOLSYS.EQ.2) FORMUL(8:8)='2'
00333       CALL VECTOR(T1,'=',FORMUL,H%ELM,-1.D0,
00334      &            HPROP,DM1,ZCONV,UDEL,VDEL,VDEL,MESH,MSK,MASKEL)
00335 !                 T1 AS HUGRADP IS NOT USED AS AN ASSEMBLED VECTOR
00336 !                 BUT TO GET THE NON ASSEMBLED FORM MESH%W
00337 !     CALCUL DES FLUX PAR SEGMENT (TE1 SUIVI DE FALSE NON UTILISE)
00338 !     FXMAT IS NOT ASSEMBLED IN //
00339 !
00340 !----------------------------------------
00341 ! DIFFERENT OPTIONS TO COMPUTE THE FLUXES
00342 !----------------------------------------
00343 !
00344       CALL FLUX_EF_VF(FXMAT,MESH%W%R,MESH%NSEG,MESH%NELEM,
00345      &                MESH%ELTSEG%I,MESH%ORISEG%I,
00346      &                MESH%IKLE%I,.TRUE.,IOPT1)
00347 !
00348 !-----------------------------------------------------
00349 ! LIMITATION OF FLUXES (IF REQUESTED, USED BY SISYPHE)
00350 !-----------------------------------------------------
00351 !
00352       IF(YAFLULIM) THEN
00353         DO I=1,MESH%NSEG
00354           FXMAT(I)=FXMAT(I)*FLULIM(I)
00355         ENDDO
00356       ENDIF
00357 !
00358       ELSE
00359         DO I=1,MESH%NSEG
00360           FXMAT(I)=GIVEN_FLUX%R(I)
00361         ENDDO
00362       ENDIF
00363 !
00364 !----------------------------------------
00365 !
00366 !     AVERAGING FLUXES ON INTERFACE SEGMENTS BY ASSEMBLING AND
00367 !     DIVIDING BY 2. THIS WILL GIVE THE UPWINDING INFORMATION
00368 !
00369       IF(NCSIZE.GT.1) THEN
00370         CALL PARCOM2_SEG(FXMAT,FXMAT,FXMAT,MESH%NSEG,1,2,1,MESH,
00371      &                   1,11)
00372         CALL MULT_INTERFACE_SEG(FXMAT,MESH%NH_COM_SEG%I,
00373      &                          MESH%NH_COM_SEG%DIM1,
00374      &                          MESH%NB_NEIGHB_SEG,
00375      &                          MESH%NB_NEIGHB_PT_SEG%I,
00376      &                          0.5D0,MESH%NSEG)
00377       ENDIF
00378 !
00379 !----------------------------------------
00380 ! END OF THE OPTIONS
00381 !----------------------------------------
00382 !
00383       CALL CPSTVC(H,T2)
00384 !     T2 WILL BE THE ASSEMBLED FLBOR, INITIALISATION HERE
00385 !     IS USELESS EXCEPT THAT PARCOM MAY ADD UNDEFINED
00386 !     NUMBERS (THAT WILL NOT BE USED BUT THAT WILL STOP
00387 !     A COMPILER... TOO BAD!)
00388       IF(NCSIZE.GT.1) CALL OS('X=0     ',X=T2)
00389 !
00390 !     INITIALIZING F AT THE OLD VALUE
00391 !
00392       CALL OS('X=Y     ',X=F,Y=FN)
00393 !
00394       CPREV=0.D0
00395       DO I=1,MESH%NSEG
00396         CPREV=CPREV+ABS(FXMAT(I))
00397       ENDDO
00398       IF(NCSIZE.GT.1) CPREV=P_DSUM(CPREV)
00399       CINIT=CPREV
00400       IF(TESTING) WRITE(LU,*) 'SOMME INITIALE DES FLUX=',CPREV
00401 !
00402 !     BOUCLE SUR LES SEGMENTS, POUR PRENDRE EN COMPTE LES FLUX
00403 !     ADMISSIBLES
00404 !
00405 !     ADDING THE SOURCES (SMH IS NATURALLY ASSEMBLED IN //)
00406       IF(YASMH) THEN
00407         IF(OPTSOU.EQ.1) THEN
00408           DO I=1,NPOIN
00409             HT%R(I)=HN%R(I)+DT*SMH%R(I)
00410             F%R(I)=FN%R(I)+DT/MAX(HT%R(I),1.D-4)*SMH%R(I)*FSCEXP%R(I)
00411           ENDDO
00412         ELSEIF(OPTSOU.EQ.2) THEN
00413           DO I=1,NPOIN
00414             HT%R(I)=HN%R(I)+DT*SMH%R(I)*UNSV2D%R(I)
00415             F%R(I)=FN%R(I)+DT/MAX(HT%R(I),1.D-4)*
00416      &                       UNSV2D%R(I)*SMH%R(I)*FSCEXP%R(I)
00417           ENDDO
00418         ENDIF
00419       ELSE
00420         DO I=1,NPOIN
00421           HT%R(I)=HN%R(I)
00422         ENDDO
00423       ENDIF
00424 !
00425 !     RAIN-EVAPORATION: RAIN FIRST, EVAPORATION IN THE END
00426 !
00427       IF(RAIN) THEN
00428         DO I=1,NPOIN
00429           C=MAX(PLUIE%R(I),0.D0)
00430           HT%R(I)=HT%R(I)+DT*C
00431 !                                                VALUE IN RAIN
00432           F%R(I)=F%R(I)+DT/MAX(HT%R(I),1.D-4)*C*(TRAIN-F%R(I))
00433         ENDDO
00434       ENDIF
00435 !
00436       IF(.NOT.YAFLBOR) THEN
00437 !       MASK=8 FOR LIQUID BOUNDARIES
00438         CALL VECTOR(FLBOR,'=','FLUBDF          ',1,1.D0,
00439      &              HPROP,HPROP,HPROP,
00440      &              UDEL , VDEL, VDEL,MESH,.TRUE.,MASKTR%ADR(8)%P)
00441       ENDIF
00442 !
00443 !     BOUNDARY FLUXES : ADDING THE ENTERING (NEGATIVE) FLUXES
00444 !     FIRST PUTTING FLBOR (BOUNDARY) IN T2 (DOMAIN)
00445       CALL OSDB( 'X=Y     ' ,T2,FLBOR,FLBOR,0.D0,MESH)
00446 !     ASSEMBLING T2 (FLBOR IS NOT ASSEMBLED)
00447       IF(NCSIZE.GT.1) CALL PARCOM(T2,2,MESH)
00448 !     POSSIBLE CORRECTION OF LIMTRA AND FBOR (LIMTRA HAS BEEN DONE BY
00449 !     DIFFIN WITH U.N)
00450       DO I=1,MESH%NPTFR
00451         N=MESH%NBOR%I(I)
00452         IF(LIMTRA(I).EQ.KDIR.AND.T2%R(N).GT.0.D0) THEN
00453           LIMTRA(I)=KDDL
00454         ELSEIF(LIMTRA(I).EQ.KDDL.AND.T2%R(N).LT.0.D0) THEN
00455           LIMTRA(I)=KDIR
00456           FBOR%R(I)=FN%R(N)
00457         ENDIF
00458       ENDDO
00459 !
00460       DO IPTFR=1,NPTFR
00461         I=NBOR(IPTFR)
00462         HT%R(I)=HT%R(I)-DT*UNSV2D%R(I)*MIN(T2%R(I),0.D0)
00463 !       ENTERING FLUXES OF TRACERS
00464 !       THE FINAL DEPTH IS TAKEN
00465         IF(LIMTRA(IPTFR).EQ.KDIR) THEN
00466           F%R(I)=FN%R(I)-DT/MAX(HT%R(I),1.D-4)*
00467      &       UNSV2D%R(I)*T2%R(I)*(FBOR%R(IPTFR)-FN%R(I))
00468         ELSEIF(LIMTRA(IPTFR).EQ.KDDL) THEN
00469 !         FLBORTRA IS NOT ASSEMBLED
00470           FLBORTRA%R(IPTFR)=FLBOR%R(IPTFR)*FN%R(I)
00471         ENDIF
00472       ENDDO
00473 !
00474 !     FOR OPTIMIZING THE LOOP ON SEGMENTS, ONLY SEGMENTS
00475 !     WITH NON ZERO FLUXES WILL BE CONSIDERED, THIS LIST
00476 !     WILL BE UPDATED. TO START WITH, ALL FLUXES ASSUMED NON ZERO
00477 !
00478       REMAIN_SEG=MESH%NSEG
00479       DO I=1,REMAIN_SEG
00480         INDIC(I)=I
00481       ENDDO
00482 !
00483       NITER = 0
00484 777   CONTINUE
00485       NITER = NITER + 1
00486 !
00487 !
00488 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00489 !     FOR DISTRIBUTING THE DEPTHS BETWEEN SEGMENTS
00490 !
00491       IF(OPTION.EQ.2) THEN
00492 !
00493 !       T1 : TOTAL FLUX REMOVED OF EACH POINT
00494 !       T4 : DEPTH H SAVED
00495 !       T6 : F SAVED
00496 !
00497         CALL CPSTVC(H,T1)
00498         IF(NITER.EQ.1) THEN
00499           DO I=1,NPOIN
00500             T1%R(I)=0.D0
00501             T4%R(I)=HT%R(I)
00502             T6%R(I)=F%R(I)
00503             T5%R(I)=HT%R(I)*F%R(I)
00504           ENDDO
00505           IF(NCSIZE.GT.1) THEN
00506             DO IPTFR=1,NPTIR
00507               I=MESH%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00508 !             AVAILABLE DEPTH IS SHARED BETWEEN PROCESSORS
00509               HT%R(I)=HT%R(I)*MESH%FAC%R(I)
00510               T5%R(I)=T5%R(I)*MESH%FAC%R(I)
00511             ENDDO
00512           ENDIF
00513         ELSE
00514 !         NOT ALL THE POINTS NEED TO BE INITIALISED NOW
00515           DO IR=1,REMAIN_SEG
00516             I=INDIC(IR)
00517             I1=GLOSEG1(I)
00518             I2=GLOSEG2(I)
00519             T1%R(I1)=0.D0
00520             T1%R(I2)=0.D0
00521 !           SAVING THE DEPTH AND TRACER
00522             T4%R(I1)=HT%R(I1)
00523             T4%R(I2)=HT%R(I2)
00524             T6%R(I1)=F%R(I1)
00525             T6%R(I2)=F%R(I2)
00526             T5%R(I1)=HT%R(I1)*F%R(I1)
00527             T5%R(I2)=HT%R(I2)*F%R(I2)
00528           ENDDO
00529 !         CANCELLING INTERFACE POINTS (SOME MAY BE ISOLATED IN A SUBDOMAIN
00530 !         AT THE TIP OF AN ACTIVE SEGMENT WHICH IS ELSEWHERE)
00531           IF(NCSIZE.GT.1) THEN
00532             DO IPTFR=1,NPTIR
00533               I=MESH%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00534               T1%R(I)=0.D0
00535 !             SAVING THE DEPTH AND TRACER
00536               T4%R(I)=HT%R(I)
00537               T6%R(I)=F%R(I)
00538 !             AVAILABLE DEPTH IS SHARED BETWEEN PROCESSORS
00539               HT%R(I)=HT%R(I)*MESH%FAC%R(I)
00540               T5%R(I)=T5%R(I)*MESH%FAC%R(I)
00541             ENDDO
00542           ENDIF
00543         ENDIF
00544         DO IR=1,REMAIN_SEG
00545           I=INDIC(IR)
00546           I1=GLOSEG1(I)
00547           I2=GLOSEG2(I)
00548           IF(FXMAT(I).GT.EPS_FLUX) THEN
00549             T1%R(I1)=T1%R(I1)+FXMAT(I)
00550             HT%R(I1)=0.D0
00551             T5%R(I1)=0.D0
00552           ELSEIF(FXMAT(I).LT.-EPS_FLUX) THEN
00553             T1%R(I2)=T1%R(I2)-FXMAT(I)
00554             HT%R(I2)=0.D0
00555             T5%R(I2)=0.D0
00556           ENDIF
00557         ENDDO
00558 !
00559         IF(NCSIZE.GT.1) CALL PARCOM(T1,2,MESH)
00560 !
00561 !       FOR ISOLATED POINTS CONNECTED TO AN ACTIVE SEGMENT
00562 !       THAT IS IN ANOTHER SUBDOMAIN
00563         IF(NCSIZE.GT.1) THEN
00564           DO IPTFR=1,NPTIR
00565             I=MESH%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00566             IF(T1%R(I).GT.EPS_FLUX) THEN
00567               HT%R(I)=0.D0
00568               T5%R(I)=0.D0
00569             ENDIF
00570           ENDDO
00571         ENDIF
00572 !
00573       ELSEIF(OPTION.EQ.1) THEN
00574 !
00575 !       AT THIS LEVEL H THE SAME AT INTERFACE POINTS
00576 !       THIS IS DONE EVEN FOR OPTION 2, TO ANTICIPATE THE FINAL PARCOM
00577         IF(NCSIZE.GT.1) THEN
00578           DO IPTFR=1,NPTIR
00579 !           AVAILABLE DEPTH IS SHARED BETWEEN PROCESSORS
00580 !           NACHB(1,IPTFR) WITH DIMENSION NACHB(NBMAXNSHARE,NPTIR)
00581             I=MESH%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00582             HT%R(I)=HT%R(I)*MESH%FAC%R(I)
00583           ENDDO
00584         ENDIF
00585 !
00586       ENDIF
00587 !
00588 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00589 !
00590       C=0.D0
00591       NEWREMAIN=0
00592 !
00593       IF(OPTION.EQ.1) THEN
00594 !
00595       DO IR=1,REMAIN_SEG
00596         I=INDIC(IR)
00597         IF(FXMAT(I).GT.EPS_FLUX) THEN
00598           I1=GLOSEG1(I)
00599           I2=GLOSEG2(I)
00600           HFL1= DT*UNSV2D%R(I1)*FXMAT(I)
00601           HFL2=-DT*UNSV2D%R(I2)*FXMAT(I)
00602 !         POUR TRACEURS
00603           H1N=HT%R(I1)
00604           H2N=HT%R(I2)
00605 !         FIN POUR TRACEURS
00606           IF(HFL1.GT.HT%R(I1)) THEN
00607             TET=HT%R(I1)/HFL1
00608             HT%R(I1)=0.D0
00609             HT%R(I2)=HT%R(I2)-HFL2*TET
00610             FXMAT(I)=FXMAT(I)*(1.D0-TET)
00611             C=C+FXMAT(I)
00612             NEWREMAIN=NEWREMAIN+1
00613             INDIC(NEWREMAIN)=I
00614           ELSE
00615             HT%R(I1)=HT%R(I1)-HFL1
00616             HT%R(I2)=HT%R(I2)-HFL2
00617           ENDIF
00618 !         TRACER (WITH TEST HT%R(I2) CANNOT BE 0.D0)
00619           IF(H2N.LT.HT%R(I2)) THEN
00620             F%R(I2)=F%R(I2)+(1.D0-H2N/HT%R(I2))*(F%R(I1)-F%R(I2))
00621           ENDIF
00622 !         END TRACER
00623         ELSEIF(FXMAT(I).LT.-EPS_FLUX) THEN
00624           I1=GLOSEG1(I)
00625           I2=GLOSEG2(I)
00626           HFL1= DT*UNSV2D%R(I1)*FXMAT(I)
00627           HFL2=-DT*UNSV2D%R(I2)*FXMAT(I)
00628 !         POUR TRACEURS
00629           H1N=HT%R(I1)
00630           H2N=HT%R(I2)
00631 !         FIN POUR TRACEURS
00632           IF(HFL2.GT.HT%R(I2)) THEN
00633             TET=HT%R(I2)/HFL2
00634             HT%R(I1)=HT%R(I1)-HFL1*TET
00635             HT%R(I2)=0.D0
00636             FXMAT(I)=FXMAT(I)*(1.D0-TET)
00637             C=C-FXMAT(I)
00638             NEWREMAIN=NEWREMAIN+1
00639             INDIC(NEWREMAIN)=I
00640           ELSE
00641             HT%R(I1)=HT%R(I1)-HFL1
00642             HT%R(I2)=HT%R(I2)-HFL2
00643           ENDIF
00644 !         TRACER (WITH TEST HT%R(I1) CANNOT BE 0.D0)
00645           IF(H1N.LT.HT%R(I1)) THEN
00646             F%R(I1)=F%R(I1)+(1.D0-H1N/HT%R(I1))*(F%R(I2)-F%R(I1))
00647           ENDIF
00648 !         FIN TRACEUR
00649         ENDIF
00650       ENDDO
00651 !
00652       ELSEIF(OPTION.EQ.2) THEN
00653 !
00654       DO IR=1,REMAIN_SEG
00655         I=INDIC(IR)
00656         I1=GLOSEG1(I)
00657         I2=GLOSEG2(I)
00658         IF(FXMAT(I).GT.EPS_FLUX) THEN
00659 !         SHARING ON DEMAND: FRACTION OF DEPTH TAKEN
00660 !         T4 IS THE STORED DEPTH
00661 !         1.D-20 ADDED HERE OTHERWISE IT OCCURED THAT THE PRODUCT OF 2
00662 !         STRICTLY POSITIVE NUMBERS GAVE 0.D0
00663           IF(T4%R(I1).GT.1.D-20) THEN
00664             HSEG1=T4%R(I1)*FXMAT(I)/T1%R(I1)
00665 !           END OF SHARING ON DEMAND
00666             HFL1= DT*UNSV2D%R(I1)*FXMAT(I)
00667             IF(HFL1.GT.HSEG1) THEN
00668               TET=HSEG1/HFL1
00669 !             HSEG2 AND THUS HT WILL BE STRICTLY POSITIVE
00670               HSEG2=DT*UNSV2D%R(I2)*FXMAT(I)*TET
00671               HT%R(I2)=HT%R(I2)+HSEG2
00672 !             GROUPING H*F
00673               T5%R(I2)=T5%R(I2)+HSEG2*T6%R(I1)
00674 !             RECOMPUTING F (AS WEIGHTED AVERAGE)
00675 !             THIS MAY BE DONE SEVERAL TIMES FOR THE SAME POINT
00676 !             BUT THE LAST ONE WILL BE THE GOOD ONE
00677               F%R(I2)=T5%R(I2)/HT%R(I2)
00678               FXMAT(I)=FXMAT(I)*(1.D0-TET)
00679               C=C+FXMAT(I)
00680               NEWREMAIN=NEWREMAIN+1
00681               INDIC(NEWREMAIN)=I
00682             ELSE
00683               HSEG1=HSEG1-HFL1
00684               HSEG2=DT*UNSV2D%R(I2)*FXMAT(I)
00685               HT%R(I2)=HT%R(I2)+HSEG2
00686               T5%R(I2)=T5%R(I2)+HSEG2*T6%R(I1)
00687 !             THE LAST ONE WILL BE THE GOOD ONE
00688               F%R(I2)=T5%R(I2)/HT%R(I2)
00689               IF(HSEG1.GT.0.D0) THEN
00690                 HT%R(I1)=HT%R(I1)+HSEG1
00691                 T5%R(I1)=T5%R(I1)+HSEG1*T6%R(I1)
00692 !               THE LAST ONE WILL BE THE GOOD ONE
00693                 F%R(I1)=T5%R(I1)/HT%R(I1)
00694               ENDIF
00695             ENDIF
00696           ELSE
00697 !           NO WATER NO FLUX TRANSMITTED, NOTHING CHANGED
00698             C=C+FXMAT(I)
00699             NEWREMAIN=NEWREMAIN+1
00700             INDIC(NEWREMAIN)=I
00701           ENDIF
00702         ELSEIF(FXMAT(I).LT.-EPS_FLUX) THEN
00703 !         SHARING ON DEMAND
00704           IF(T4%R(I2).GT.1.D-20) THEN
00705 !         1.D-20 ADDED HERE OTHERWISE IT OCCURED THAT THE PRODUCT OF 2
00706 !         STRICTLY POSITIVE NUMBERS GAVE 0.D0
00707             HSEG2=-T4%R(I2)*FXMAT(I)/T1%R(I2)
00708 !           END OF SHARING ON DEMAND
00709             HFL2=-DT*UNSV2D%R(I2)*FXMAT(I)
00710             IF(HFL2.GT.HSEG2) THEN
00711               TET=HSEG2/HFL2
00712 !             HSEG1 AND THUS HT WILL BE STRICTLY POSITIVE
00713               HSEG1=-DT*UNSV2D%R(I1)*FXMAT(I)*TET
00714               HT%R(I1)=HT%R(I1)+HSEG1
00715               T5%R(I1)=T5%R(I1)+HSEG1*T6%R(I2)
00716 !             THE LAST ONE WILL BE THE GOOD ONE
00717               F%R(I1)=T5%R(I1)/HT%R(I1)
00718               FXMAT(I)=FXMAT(I)*(1.D0-TET)
00719               C=C-FXMAT(I)
00720               NEWREMAIN=NEWREMAIN+1
00721               INDIC(NEWREMAIN)=I
00722             ELSE
00723               HSEG1=-DT*UNSV2D%R(I1)*FXMAT(I)
00724               HSEG2=HSEG2-HFL2
00725               HT%R(I1)=HT%R(I1)+HSEG1
00726               T5%R(I1)=T5%R(I1)+HSEG1*T6%R(I2)
00727               F%R(I1)=T5%R(I1)/HT%R(I1)
00728               IF(HSEG2.GT.0.D0) THEN
00729                 HT%R(I2)=HT%R(I2)+HSEG2
00730                 T5%R(I2)=T5%R(I2)+HSEG2*T6%R(I2)
00731 !               THE LAST ONE WILL BE THE GOOD ONE
00732                 F%R(I2)=T5%R(I2)/HT%R(I2)
00733               ENDIF
00734             ENDIF
00735           ELSE
00736 !           NO WATER NO FLUX TRANSMITTED, NOTHING CHANGED
00737             C=C-FXMAT(I)
00738             NEWREMAIN=NEWREMAIN+1
00739             INDIC(NEWREMAIN)=I
00740           ENDIF
00741         ENDIF
00742       ENDDO
00743 !
00744 !     ELSE
00745 !       UNKNOWN OPTION
00746       ENDIF
00747 !
00748       REMAIN_SEG=NEWREMAIN
00749 !
00750 !     MERGING DEPTHS AND F AT INTERFACE POINTS
00751 !
00752       IF(NCSIZE.GT.1) THEN
00753         DO IPTFR=1,NPTIR
00754 !         ARRAY WITH HT*F AT INTERFACE POINTS
00755           I=MESH%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00756           T1%R(I)=HT%R(I)*F%R(I)
00757         ENDDO
00758 !       SUMMING HT*F AT INTERFACE POINTS
00759         CALL PARCOM(T1,2,MESH)
00760 !       SUMMING THE NEW POSITIVE PARTIAL DEPTHS OF INTERFACE POINTS
00761         CALL PARCOM(HT,2,MESH)
00762 !       AVERAGE F AT INTERFACE POINTS
00763         DO IPTFR=1,NPTIR
00764           I=MESH%NACHB%I(NBMAXNSHARE*(IPTFR-1)+1)
00765           IF(HT%R(I).GT.0.D0) F%R(I)=T1%R(I)/HT%R(I)
00766         ENDDO
00767       ENDIF
00768 !
00769       IF(NCSIZE.GT.1) C=P_DSUM(C)
00770       IF(TESTING) WRITE(LU,*) 'FLUX NON PRIS EN COMPTE=',C
00771       IF(C.NE.CPREV.AND.ABS(C-CPREV).GT.CINIT*1.D-9
00772      &             .AND.C.NE.0.D0) THEN
00773         CPREV=C
00774         IF(NITER.LT.NITMAX) GO TO 777
00775       ENDIF
00776 !
00777 !     RAIN-EVAPORATION: RAIN DONE ABOVE, NOW EVAPORATION
00778 !
00779       IF(RAIN) THEN
00780         DO I=1,NPOIN
00781           C=MIN(PLUIE%R(I),0.D0)
00782 !         POSITIVITY NOT TESTED HERE, WOULD REQUIRE C=MAX(C,-HT%R(I)/DT)
00783 !         BUT THEN MASS-BALANCE WOULD NOT BE CORRECT,
00784           HT%R(I)=HT%R(I)+DT*C
00785 !                                                VALUE IN VAPOR
00786 !         F%R(I)=F%R(I)+DT/MAX(HT%R(I),1.D-4)*C*(0.D0-F%R(I))
00787           F%R(I)=F%R(I)-DT/MAX(HT%R(I),1.D-4)*C*F%R(I)
00788         ENDDO
00789       ENDIF
00790 !
00791 !     BOUNDARY FLUXES : ADDING THE EXITING (POSITIVE) FLUXES
00792 !                       WITH A POSSIBLE LIMITATION
00793 !
00794       DO IPTFR=1,NPTFR
00795         I=NBOR(IPTFR)
00796 !                               T2 = // ASSEMBLED FLBOR
00797         HFL1=DT*UNSV2D%R(I)*MAX(T2%R(I),0.D0)
00798         TET=1.D0
00799 !       NEXT LINE SHOULD NEVER HAPPEN (DONE IN POSITIVE_DEPTHS)
00800         IF(HFL1.GT.HT%R(I)) TET=HT%R(I)/HFL1
00801 !       MAX IS ONLY TO PREVENT TRUNCATION ERROR
00802         HT%R(I)=MAX(HT%R(I)-HFL1*TET,0.D0)
00803 !       LIMITATION OF FLBOR (MUST HAVE BEEN DONE ALREADY
00804 !                            IN POSITIVE_DEPTHS)
00805 !       FLBOR%R(IPTFR)=FLBOR%R(IPTFR)*TET
00806         IF(LIMTRA(IPTFR).EQ.KDIR) THEN
00807           F%R(I)=F%R(I)-HFL1*TET/MAX(HT%R(I),1.D-4)*
00808      &           (FBOR%R(IPTFR)-F%R(I))
00809           FLBORTRA%R(IPTFR)=FLBOR%R(IPTFR)*FBOR%R(IPTFR)
00810         ELSEIF(LIMTRA(IPTFR).EQ.KDDL) THEN
00811           FLBORTRA%R(IPTFR)=FLBOR%R(IPTFR)*F%R(I)
00812         ELSE
00813           FLBORTRA%R(IPTFR)=0.D0
00814         ENDIF
00815       ENDDO
00816 !
00817       IF(TESTING) THEN
00818         C=0.D0
00819         DO I=1,NPOIN
00820           C=C+(HT%R(I)-H%R(I))**2
00821         ENDDO
00822 !                       FAUX MAIS PAS GRAVE SI 0.
00823         IF(NCSIZE.GT.1) C=P_DSUM(C)
00824         WRITE(LU,*) 'DIFFERENCE ENTRE H ET HT =',C
00825 !
00826         C=1.D99
00827         DO I=1,NPOIN
00828           C=MIN(C,F%R(I))
00829         ENDDO
00830         IF(NCSIZE.GT.1) C=P_DMIN(C)
00831         WRITE(LU,*) 'APRES TRAITEMENT TRACEUR MIN=',C
00832         C=-1.D99
00833         DO I=1,NPOIN
00834           C=MAX(C,F%R(I))
00835         ENDDO
00836         IF(NCSIZE.GT.1) C=P_DMAX(C)
00837         WRITE(LU,*) 'APRES TRAITEMENT TRACEUR MAX=',C
00838       ENDIF
00839 !
00840 !-----------------------------------------------------------------------
00841 !
00842 !     SOURCE TERMS
00843 !
00844       IF(YASMI) THEN
00845 !
00846 !       IMPLICIT AND EXPLICIT SOURCE TERM
00847 !
00848         IF(IOPT2.EQ.0) THEN
00849           DO I = 1,MESH%NPOIN
00850             F%R(I)=(F%R(I)+DT*SM%R(I))/
00851      &           (1.D0-DT*SMI%R(I)/MAX(H%R(I),1.D-15))
00852 !           COULD BE DONE LIKE THIS...
00853 !           F%R(I)=H%R(I)*(F%R(I)+DT*SM%R(I))/(H%R(I)-DT*SMI%R(I))
00854           ENDDO
00855         ELSEIF(IOPT2.EQ.1) THEN
00856 !         HERE WE ASSUME THAT SMI WILL PREVENT A DIVISION BY ZERO
00857 !         THIS IS THE CASE WITH SETTLING VELOCITY IN SISYPHE
00858           DO I = 1,MESH%NPOIN
00859           F%R(I)=(F%R(I)*HT%R(I)+DT*SM%R(I)*H%R(I))/(H%R(I)-DT*SMI%R(I))
00860           ENDDO
00861         ENDIF
00862 !
00863       ELSE
00864 !
00865 !       EXPLICIT SOURCE TERM ONLY (AND IOPT2=1 NOT TREATED !!!)
00866 !
00867         DO I = 1,MESH%NPOIN
00868           F%R(I) = F%R(I)+DT*SM%R(I)
00869         ENDDO
00870 !
00871       ENDIF
00872 !
00873 !-----------------------------------------------------------------------
00874 !
00875       IF(ENTET) THEN
00876         IF(LNG.EQ.1) WRITE(LU,101) NITER
00877         IF(LNG.EQ.2) WRITE(LU,102) NITER
00878       ENDIF
00879 !
00880 101   FORMAT(' CVTRVF_POS (SCHEMA 13 OU 14) : ',1I3,' ITERATIONS')
00881 102   FORMAT(' CVTRVF_POS (SCHEME 13 OR 14): ',1I3,' ITERATIONS')
00882 !
00883 !-----------------------------------------------------------------------
00884 !
00885       RETURN
00886       END

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