cvtrvf_pos_2.f

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

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