difsou.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\difsou.f
00002 !
00133                      SUBROUTINE DIFSOU
00134 !                    *****************
00135 !
00136      &(TEXP,TIMP,YASMI,TSCEXP,HPROP,TN,TETAT,NREJTR,ISCE,DSCE,TSCE,
00137      & MAXSCE,MAXTRA,AT,DT,MASSOU,NTRAC,FAC,NSIPH,ENTSIP,SORSIP,
00138      & DSIP,TSIP,NBUSE,ENTBUS,SORBUS,DBUS,TBUS,NWEIRS,TYPSEUIL,
00139      & NPSING,NDGA1,NDGB1,TWEIRA,TWEIRB)
00140 !
00141 !***********************************************************************
00142 ! TELEMAC2D   V7P0
00143 !***********************************************************************
00144 !
00145 !
00146 !
00147 !
00148 !
00149 !
00150 !
00151 !
00152 !
00153 !
00154 !
00155 !
00156 !
00157 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00158 !| AT             |-->| TIME IN SECONDS
00159 !| DBUS           |-->| DISCHARGE OF TUBES.
00160 !| DSCE           |-->| DISCHARGE OF POINT SOURCES
00161 !| DSIP           |-->| DISCHARGE OF CULVERT.
00162 !| DT             |-->| TIME STEP
00163 !| ENTBUS         |-->| INDICES OF ENTRY OF TUBES IN GLOBAL NUMBERING
00164 !| ENTSIP         |-->| INDICES OF ENTRY OF PIPE IN GLOBAL NUMBERING
00165 !| FAC            |-->| IN PARALLEL :
00166 !|                |   | 1/(NUMBER OF SUB-DOMAINS OF THE POINT)
00167 !| HPROP          |-->| PROPAGATION DEPTH
00168 !| ISCE           |-->| NEAREST POINTS OF DISCHARGES
00169 !| MASSOU         |<--| MASS OF TRACER ADDED BY SOURCE TERM
00170 !| MAXSCE         |-->| MAXIMUM NUMBER OF SOURCES
00171 !| MAXTRA         |-->| MAXIMUM NUMBER OF TRACERS
00172 !| NBUSE          |-->| NUMBER OF TUBES
00173 !| NREJTR         |-->| NUMBER OF POINT SOURCES AS GIVEN BY TRACERS KEYWORDS
00174 !| NSIPH          |-->| NUMBER OF CULVERTS
00175 !| NTRAC          |-->| NUMBER OF TRACERS
00176 !| NWEIRS         |-->| NUMBER OF WEIRS
00177 !| SORBUS         |-->| INDICES OF TUBES EXITS IN GLOBAL NUMBERING
00178 !| SORSIP         |-->| INDICES OF PIPES EXITS IN GLOBAL NUMBERING
00179 !| TBUS           |-->| VALUES OF TRACERS AT TUBES EXTREMITY
00180 !| TETAT          |-->| COEFFICIENT OF IMPLICITATION FOR TRACERS.
00181 !| TEXP           |-->| EXPLICIT SOURCE TERM.
00182 !| TIMP           |-->| IMPLICIT SOURCE TERM.
00183 !| TN             |-->| TRACERS AT TIME N
00184 !| TSCE           |-->| PRESCRIBED VALUES OF TRACERS AT POINT SOURCES
00185 !| TSCEXP         |<--| EXPLICIT SOURCE TERM OF POINT SOURCES
00186 !|                |   | IN TRACER EQUATION, EQUAL TO:
00187 !|                |   | TSCE - ( 1 - TETAT ) TN
00188 !| TSIP           |-->| VALUES OF TRACERS AT CULVERT EXTREMITY
00189 !| TWEIRA         |-->| VALUES OF TRACERS ON SIDE A OF WEIR
00190 !| TWEIRB         |-->| VALUES OF TRACERS ON SIDE B OF WEIR
00191 !| TYPSEUIL       |-->| TYPE OF WEIRS (IF = 2, WEIRS TREATED AS SOURCES POINTS)
00192 !| YASMI          |<--| IF YES, THERE ARE IMPLICIT SOURCE TERMS
00193 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00194 !
00195       USE BIEF
00196       USE INTERFACE_PARALLEL
00197       USE DECLARATIONS_TELEMAC2D, ONLY : LOITRAC, COEF1TRAC, QWA, QWB,
00198      &   MAXNPS,U,V,UNSV2D,V2DPAR,VOLU2D,T1,T2,T3,T4,MESH,MSK,MASKEL,
00199      &   IELMU,S,NPOIN,CF,H,SECCURRENTS,SEC_AS,SEC_DS,SEC_R,WATQUA,IND_T
00200       USE DECLARATIONS_WAQTEL,ONLY: WAQPROCESS,FORMRS,O2SATU,ADDTR,
00201      &                              WATTEMP,RSW,ABRS
00202 !
00203       IMPLICIT NONE
00204       INTEGER LNG,LU
00205       COMMON/INFO/LNG,LU
00206 !
00207 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00208 !
00209       INTEGER          , INTENT(IN)    :: ISCE(*),NREJTR,NTRAC
00210       INTEGER          , INTENT(IN)    :: NSIPH,NBUSE,NWEIRS
00211       INTEGER          , INTENT(IN)    :: ENTSIP(NSIPH),SORSIP(NSIPH)
00212       INTEGER          , INTENT(IN)    :: ENTBUS(NBUSE),SORBUS(NBUSE)
00213       INTEGER          , INTENT(IN)    :: MAXSCE,MAXTRA,TYPSEUIL
00214       LOGICAL          , INTENT(INOUT) :: YASMI(*)
00215       DOUBLE PRECISION , INTENT(IN)    :: AT,DT,TETAT,DSCE(*)
00216       DOUBLE PRECISION , INTENT(IN)    :: DSIP(NSIPH),DBUS(NBUSE)
00217       DOUBLE PRECISION , INTENT(IN)    :: TSCE(MAXSCE,MAXTRA),FAC(*)
00218       DOUBLE PRECISION , INTENT(INOUT) :: MASSOU(*)
00219       TYPE(BIEF_OBJ)   , INTENT(IN)    :: TN,HPROP,TSIP,TBUS
00220       TYPE(BIEF_OBJ)   , INTENT(IN)    :: TWEIRA,TWEIRB
00221       TYPE(BIEF_OBJ)   , INTENT(IN)    :: NPSING,NDGA1,NDGB1
00222       TYPE(BIEF_OBJ)   , INTENT(INOUT) :: TSCEXP,TEXP,TIMP
00223 !
00224 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00225 !
00226       INTEGER I,K,IR,ITRAC,N,INDIC,NTRA
00227 !
00228       DOUBLE PRECISION DEBIT,TRASCE
00229       DOUBLE PRECISION DENOM,NUMER,NORM2,SEC_RMAX,RMAX
00230 !
00231       DOUBLE PRECISION H1,H2,TRUP,TRDO,AB,DZ
00232       DOUBLE PRECISION, PARAMETER :: EPS=1.D-6
00233 !
00234 !      DOUBLE PRECISION P_DSUM,P_DMIN,P_DMAX
00235 !      EXTERNAL         P_DSUM,P_DMIN,P_DMAX
00236 !
00237       INTRINSIC SQRT
00238 !
00239 !-----------------------------------------------------------------------
00240 !
00241 !     SECONDARY CURRENTS WILL BE TREATED APART
00242 !
00243       NTRA=NTRAC
00244       IF(SECCURRENTS) NTRA=NTRA-1
00245 !
00246 !-----------------------------------------------------------------------
00247 !
00248 !     EXPLICIT SOURCE TERMS
00249 !
00250       DO ITRAC=1,NTRA
00251         CALL OS('X=0     ',X=TSCEXP%ADR(ITRAC)%P)
00252         CALL OS('X=0     ',X=TEXP%ADR(ITRAC)%P)
00253         MASSOU(ITRAC) = 0.D0
00254       ENDDO
00255 !
00256 !-----------------------------------------------------------------------
00257 !
00258 !     IMPLICIT SOURCE TERMS (DEPENDING ON THE LAW CHOSEN)
00259 !
00260       DO ITRAC=1,NTRA
00261         IF(LOITRAC(ITRAC).EQ.0) THEN
00262           YASMI(ITRAC)=.FALSE.
00263         ELSEIF(LOITRAC(ITRAC).EQ.1) THEN
00264           YASMI(ITRAC)=.TRUE.
00265           CALL OS('X=CY    ',X=TIMP%ADR(ITRAC)%P,Y=HPROP,
00266      &            C=-2.3D0/COEF1TRAC(ITRAC)/3600.D0)
00267         ELSE
00268           IF(LNG.EQ.1) WRITE(LU,*) 'DIFSOU : LOI NON PROGRAMMEE'
00269           IF(LNG.EQ.2) WRITE(LU,*) 'DIFSOU : LAW NOT IMPLEMENTED'
00270         ENDIF
00271       ENDDO
00272 !
00273 !                                   N+1
00274 !     EXAMPLE WHERE WE ADD -0.0001 F      IN THE RIGHT HAND-SIDE
00275 !     OF THE TRACER EQUATION THAT BEGINS WITH DF/DT=...
00276 !     (T12=SMI WILL BE DIVIDED BY HPROP IN CVDFTR, THE EQUATION IS:
00277 !     DT/DT=...+SMI*T(N+1)/H
00278 !
00279 !     HERE THIS IS DONE FOR TRACER 3 ONLY IN A RECTANGULAR ZONE
00280 !
00281 !     CALL OS('X=0     ',X=TIMP%ADR(3)%P)
00282 !     DO I=1,HPROP%DIM1
00283 !       IF(X(I).GE.263277.D0.AND.X(I).LE.265037.D0) THEN
00284 !       IF(Y(I).GE.379007.D0.AND.Y(I).LE.380326.D0) THEN
00285 !         TIMP%ADR(3)%P%R(I)=-0.00001D0*HPROP%R(I)
00286 !       ENDIF
00287 !       ENDIF
00288 !     ENDDO
00289 !     YASMI(3)=.TRUE.
00290 !
00291 !-----------------------------------------------------------------------
00292 !
00293 !  TAKES THE SOURCES OF TRACER INTO ACCOUNT
00294 !
00295 !-----------------------------------------------------------------------
00296 !
00297       DO ITRAC=1,NTRA
00298 !
00299         IF(NREJTR.GT.0) THEN
00300 !
00301           DO I = 1 , NREJTR
00302 !
00303             IR = ISCE(I)
00304 !           TEST IR.GT.0 FOR THE PARALLELISM
00305             IF(IR.GT.0) THEN
00306               DEBIT=DSCE(I)
00307               IF(DEBIT.GT.0.D0) THEN
00308                 TRASCE = TSCE(I,ITRAC)
00309               ELSE
00310 !               THE VALUE AT THE SOURCE IS TN IF THE FLOW IS OUTGOING
00311                 TRASCE = TN%ADR(ITRAC)%P%R(IR)
00312               ENDIF
00313 !             SOURCE TERM ADDED TO THE MASS OF TRACER
00314               IF(NCSIZE.GT.1) THEN
00315 !               FAC TO AVOID COUNTING THE POINT SEVERAL TIMES
00316 !               (SEE CALL TO P_DSUM BELOW)
00317                 MASSOU(ITRAC)=MASSOU(ITRAC)+DT*DEBIT*TRASCE*FAC(IR)
00318               ELSE
00319                 MASSOU(ITRAC)=MASSOU(ITRAC)+DT*DEBIT*TRASCE
00320               ENDIF
00321               TRASCE = TRASCE - (1.D0 - TETAT) * TN%ADR(ITRAC)%P%R(IR)
00322               TSCEXP%ADR(ITRAC)%P%R(IR)=TSCEXP%ADR(ITRAC)%P%R(IR)+TRASCE
00323 !
00324 !             THE IMPLICIT PART OF THE TERM - T * SCE
00325 !             IS DEALT WITH IN CVDFTR.
00326 !
00327             ENDIF
00328 !
00329           ENDDO
00330 !
00331         ENDIF
00332 !
00333         IF(NSIPH.GT.0) THEN
00334           DO I = 1 , NSIPH
00335             IR = ENTSIP(I)
00336             IF(IR.GT.0) THEN
00337               IF(NCSIZE.GT.1) THEN
00338 !               FAC TO AVOID COUNTING THE POINT SEVERAL TIMES
00339 !               (SEE CALL TO P_DSUM BELOW)
00340                 MASSOU(ITRAC)=MASSOU(ITRAC)-DT*DSIP(I)*
00341      &                        TSIP%ADR(ITRAC)%P%R(I)*FAC(IR)
00342               ELSE
00343                 MASSOU(ITRAC)=MASSOU(ITRAC)-DT*DSIP(I)*
00344      &                        TSIP%ADR(ITRAC)%P%R(I)
00345               ENDIF
00346               TSCEXP%ADR(ITRAC)%P%R(IR)=TSCEXP%ADR(ITRAC)%P%R(IR) +
00347      &           TSIP%ADR(ITRAC)%P%R(I) -
00348      &           (1.D0 - TETAT) * TN%ADR(ITRAC)%P%R(IR)
00349             ENDIF
00350             IR = SORSIP(I)
00351             IF(IR.GT.0) THEN
00352               IF(NCSIZE.GT.1) THEN
00353 !               FAC TO AVOID COUNTING THE POINT SEVERAL TIMES
00354 !               (SEE CALL TO P_DSUM BELOW)
00355                 MASSOU(ITRAC)=MASSOU(ITRAC)+DT*DSIP(I)*
00356      &                        TSIP%ADR(ITRAC)%P%R(NSIPH+I)*FAC(IR)
00357               ELSE
00358                 MASSOU(ITRAC)=MASSOU(ITRAC)+DT*DSIP(I)*
00359      &                        TSIP%ADR(ITRAC)%P%R(NSIPH+I)
00360               ENDIF
00361               TSCEXP%ADR(ITRAC)%P%R(IR)=TSCEXP%ADR(ITRAC)%P%R(IR) +
00362      &           TSIP%ADR(ITRAC)%P%R(NSIPH+I) -
00363      &           (1.D0 - TETAT) * TN%ADR(ITRAC)%P%R(IR)
00364             ENDIF
00365           ENDDO
00366         ENDIF
00367 !
00368         IF(NBUSE.GT.0) THEN
00369           DO I = 1 , NBUSE
00370             IR = ENTBUS(I)
00371             IF(IR.GT.0) THEN
00372               IF(NCSIZE.GT.1) THEN
00373 !               FAC TO AVOID COUNTING THE POINT SEVERAL TIMES
00374 !               (SEE CALL TO P_DSUM BELOW)
00375                 MASSOU(ITRAC)=MASSOU(ITRAC)-DT*DBUS(I)*
00376      &                        TBUS%ADR(ITRAC)%P%R(I)*FAC(IR)
00377               ELSE
00378                 MASSOU(ITRAC)=MASSOU(ITRAC)-DT*DBUS(I)*
00379      &                        TBUS%ADR(ITRAC)%P%R(I)
00380               ENDIF
00381               TSCEXP%ADR(ITRAC)%P%R(IR)=TSCEXP%ADR(ITRAC)%P%R(IR) +
00382      &           TBUS%ADR(ITRAC)%P%R(I) -
00383      &           (1.D0 - TETAT) * TN%ADR(ITRAC)%P%R(IR)
00384             ENDIF
00385             IR = SORBUS(I)
00386             IF(IR.GT.0) THEN
00387               IF(NCSIZE.GT.1) THEN
00388 !               FAC TO AVOID COUNTING THE POINT SEVERAL TIMES
00389 !               (SEE CALL TO P_DSUM BELOW)
00390                 MASSOU(ITRAC)=MASSOU(ITRAC)+DT*DBUS(I)*
00391      &                        TBUS%ADR(ITRAC)%P%R(NBUSE+I)*FAC(IR)
00392               ELSE
00393                 MASSOU(ITRAC)=MASSOU(ITRAC)+DT*DBUS(I)*
00394      &                        TBUS%ADR(ITRAC)%P%R(NBUSE+I)
00395               ENDIF
00396               TSCEXP%ADR(ITRAC)%P%R(IR)=TSCEXP%ADR(ITRAC)%P%R(IR) +
00397      &           TBUS%ADR(ITRAC)%P%R(NBUSE+I) -
00398      &           (1.D0 - TETAT) * TN%ADR(ITRAC)%P%R(IR)
00399             ENDIF
00400           ENDDO
00401         ENDIF
00402 !
00403         IF(NWEIRS.GT.0.AND.TYPSEUIL.EQ.2) THEN
00404           DO N = 1 , NWEIRS
00405             DO I = 1 , NPSING%I(N)
00406               INDIC = (N-1)*MAXNPS + I
00407               IR = NDGA1%ADR(N)%P%I(I)
00408               H1 = 0.D0
00409               TRUP = 0.D0
00410               IF(IR.GT.0) THEN
00411                 IF(NCSIZE.GT.1) THEN
00412 !                 FAC TO AVOID COUNTING THE POINT SEVERAL TIMES
00413 !                 (SEE CALL TO P_DSUM BELOW)
00414                   MASSOU(ITRAC)=MASSOU(ITRAC)-DT*QWA%ADR(N)%P%R(I)*
00415      &                          TWEIRA%ADR(ITRAC)%P%R(INDIC)*FAC(IR)
00416                 ELSE
00417                   MASSOU(ITRAC)=MASSOU(ITRAC)-DT*QWA%ADR(N)%P%R(I)*
00418      &                          TWEIRA%ADR(ITRAC)%P%R(INDIC)
00419                 ENDIF
00420                 TSCEXP%ADR(ITRAC)%P%R(IR)=TSCEXP%ADR(ITRAC)%P%R(IR) +
00421      &             TWEIRA%ADR(ITRAC)%P%R(INDIC) -
00422      &             (1.D0 - TETAT) * TN%ADR(ITRAC)%P%R(IR)
00423 !               RECUPERATE H FOR WAQ
00424                 IF(WATQUA.AND.WAQPROCESS.EQ.1)THEN
00425                   H1   = HPROP%R(IR)
00426                   TRUP = TN%ADR(NTRAC-ADDTR+1)%P%R(IR)
00427                   IF(NCSIZE.GT.1)THEN
00428                     H1   = P_DMIN(H1  )+P_DMAX(H1  )
00429                     TRUP = P_DMIN(TRUP)+P_DMAX(TRUP)
00430                   ENDIF
00431                 ENDIF
00432               ENDIF
00433               IR = NDGB1%ADR(N)%P%I(I)
00434               H2   = 0.D0
00435               TRDO = 0.D0
00436               IF(IR.GT.0) THEN
00437                 IF(NCSIZE.GT.1) THEN
00438 !                 FAC TO AVOID COUNTING THE POINT SEVERAL TIMES
00439 !                 (SEE CALL TO P_DSUM BELOW)
00440                   MASSOU(ITRAC)=MASSOU(ITRAC)+DT*QWB%ADR(N)%P%R(I)*
00441      &                          TWEIRB%ADR(ITRAC)%P%R(INDIC)*FAC(IR)
00442                 ELSE
00443                   MASSOU(ITRAC)=MASSOU(ITRAC)+DT*QWB%ADR(N)%P%R(I)*
00444      &                          TWEIRB%ADR(ITRAC)%P%R(INDIC)
00445                 ENDIF
00446                 TSCEXP%ADR(ITRAC)%P%R(IR)=TSCEXP%ADR(ITRAC)%P%R(IR) +
00447      &             TWEIRB%ADR(ITRAC)%P%R(INDIC) -
00448      &             (1.D0 - TETAT) * TN%ADR(ITRAC)%P%R(IR)
00449 !               RECUPERATE H FOR WAQ
00450                 IF(WATQUA.AND.WAQPROCESS.EQ.1)THEN
00451                   H2  = HPROP%R(IR)
00452                   IF(NCSIZE.GT.1)THEN
00453                     H2   = P_DMIN(H2  )+P_DMAX(H2  )
00454                   ENDIF
00455                 ENDIF
00456 !               CONTRIBUTION TO WAQ
00457                 IF(WATQUA.AND.WAQPROCESS.EQ.1)THEN
00458                   DZ= ABS(H2-H1)
00459                   RSW = 0.D0
00460 !                 LETS'S COMPUTE RS IF IT IS NOT TAKEN CONSTANT
00461                   IF(FORMRS.NE.0) THEN
00462                     AB=ABRS(1)*ABRS(2)
00463 !                   GAMESON FORMULA 1
00464                     IF(FORMRS.EQ.1)     THEN
00465                       RSW=1.D0+0.5D0*AB*DZ
00466 !                   GAMESON FORMULA2
00467                     ELSEIF(FORMRS.EQ.2) THEN
00468                       RSW = 0.11D0*AB*(1.D0+0.046D0*WATTEMP)*DZ
00469 !                   WRL FORMULA 1 (NO NEED TO AB ? )
00470                     ELSEIF(FORMRS.EQ.3 )THEN
00471                       RSW = 1.D0+0.69D0*DZ*(1.D0-0.11D0*DZ )
00472      &                      *( 1.D0+0.046D0*WATTEMP)
00473 !                   WRL FORMULA 2
00474                     ELSEIF (FORMRS.EQ.4)THEN
00475                       RSW = 1.D0+0.38D0*AB*DZ*(1.D0-0.11D0*DZ)
00476      &                      * (1.D0+0.046D0*WATTEMP )
00477                     ELSE
00478                       IF(LNG.EQ.1)THEN
00479                         WRITE(LU,*)'FORMULE DE RS (REAERATION AU SEUIL)'
00480                         WRITE(LU,*)'INCONNUE  :',FORMRS
00481                         WRITE(LU,*)'LES CHOIX POSSIBLES SONT DE 1 A 4'
00482                       ELSE
00483                         WRITE(LU,*)'FORMULA FOR RS (WEIR REAERATION) '
00484                         WRITE(LU,*)' NOT VALID  :',FORMRS
00485                         WRITE(LU,*)'POSSIBLE CHOICES ARE FROM 1 TO 4'
00486                       ENDIF
00487                       CALL PLANTE(1)
00488                       STOP
00489                     ENDIF
00490 !
00491 !                   FORCING O2 DENSITY DOWNSTREAM THE WEIR
00492 !
00493                     IF(ABS(RSW).GT.EPS)THEN
00494                       TRDO = O2SATU + (TRUP-O2SATU)/RSW
00495                     ELSE
00496                       WRITE(LU,*)'DIFSOU:RSW VERY SMALL',RSW
00497                       CALL PLANTE(1)
00498                       STOP
00499                     ENDIF
00500                     IF(NCSIZE.GT.1)TRDO = P_DMIN(TRDO)+P_DMAX(TRDO)
00501                     TN%ADR(NTRAC-ADDTR+1)%P%R(IR)=TRDO
00502                   ENDIF
00503                 ENDIF
00504               ENDIF
00505             ENDDO
00506           ENDDO
00507         ENDIF
00508 !
00509         IF(NCSIZE.GT.1.AND.
00510      &     (NREJTR.GT.0.OR.NSIPH.GT.0.OR.NBUSE.GT.0.OR.
00511      &      (NWEIRS.GT.0.AND.TYPSEUIL.EQ.2))) THEN
00512           MASSOU(ITRAC)=P_DSUM(MASSOU(ITRAC))
00513         ENDIF
00514 !
00515       ENDDO
00516 !
00517 !     WATER QUALITY CONTRIBUTION TO TRACER SOURCES
00518       IF(WATQUA)THEN
00519         CALL SOURCE_WAQ(NPOIN,TEXP,TIMP,YASMI,TSCEXP,HPROP,TN,TETAT,
00520      &                  AT,DT,NTRAC,WAQPROCESS)
00521       ENDIF
00522 !
00523 !-----------------------------------------------------------------------
00524 !
00525 !     SECONDARY CURRENTS (OMEGA IS THE TRACER OF RANK NTRAC)
00526 !
00527       IF(SECCURRENTS) THEN
00528 !
00529         CALL VECTOR(T1,'=','GRADF          Y',IELMU,
00530      &              1.D0,V,S,S,S,S,S,MESH,MSK,MASKEL)
00531         CALL VECTOR(T2,'=','GRADF          X',IELMU,
00532      &              1.D0,U,S,S,S,S,S,MESH,MSK,MASKEL)
00533         CALL VECTOR(T3,'=','GRADF          X',IELMU,
00534      &              1.D0,V,S,S,S,S,S,MESH,MSK,MASKEL)
00535         CALL VECTOR(T4,'=','GRADF          Y',IELMU,
00536      &              1.D0,U,S,S,S,S,S,MESH,MSK,MASKEL)
00537         IF(NCSIZE.GT.1) THEN
00538           CALL PARCOM (T1, 2, MESH)
00539           CALL PARCOM (T2, 2, MESH)
00540           CALL PARCOM (T3, 2, MESH)
00541           CALL PARCOM (T4, 2, MESH)
00542         ENDIF
00543 !
00544 !       INITIALISATIONS
00545 !
00546         CALL OS('X=0     ',X=TSCEXP%ADR(NTRAC)%P)
00547         YASMI(NTRAC)=.TRUE.
00548 !
00549 !       SOURCE TERMS
00550 !
00551         DO K=1,NPOIN
00552           NORM2=U%R(K)**2+V%R(K)**2
00553           NUMER = (U%R(K)*V%R(K)*(T1%R(K)-T2%R(K))+U%R(K)**2*(T3%R(K))
00554      &           -V%R(K)**2*(T4%R(K)))*UNSV2D%R(K)
00555           SEC_R%R(K)=NUMER/MAX(SQRT(NORM2)**3,1.D-9)
00556 !         GEOMETRY: R OBVIOUSLY LARGER THAN 0.5 LOCAL MESH SIZE
00557 !         THEORY ALSO SAYS R > 2H
00558 !         LOCAL MESH SIZE HERE ASSUMED TO BE SQRT(V2DPAR)
00559           RMAX=MAX(2.D0*H%R(K),0.5D0*SQRT(V2DPAR%R(K)))
00560 !         RMAX=0.5D0*SQRT(V2DPAR%R(K))
00561           SEC_RMAX=1.D0/RMAX
00562           SEC_R%R(K)=MAX(-SEC_RMAX,MIN(SEC_RMAX,SEC_R%R(K)))
00563 !         EXPLICIT SOURCE TERMS (CREATION OF OMEGA)
00564 !         CLIPPING OF H AT 1.D-2
00565 !         NOTE: IMPLICIT TERMS (DESTRUCTION) IN CVDFTR CLIPPED AT 1.D-4
00566           DENOM=MAX(H%R(K),1.D-2)*(9.D0*(H%R(K)*SEC_R%R(K))**2+1.D0)
00567           TEXP%ADR(NTRAC)%P%R(K)=
00568      &                 SEC_AS*SQRT(0.5D0*CF%R(K))*NORM2*SEC_R%R(K)/DENOM
00569 !         IMPLICIT SOURCE TERMS (DEPENDING ON THE LAW CHOSEN)
00570           TIMP%ADR(NTRAC)%P%R(K)=-SEC_DS*SQRT(0.5D0*CF%R(K)*NORM2)
00571         ENDDO
00572 !
00573 !       MASS ADDED BY EXPLICIT TERMS
00574 !       THE MASS ADDED BY IMPLICIT TERMS IS COMPUTED IN CVDFTR
00575 !
00576         MASSOU(NTRAC) = 0.D0
00577         DO K=1,NPOIN
00578           MASSOU(NTRAC)=MASSOU(NTRAC)
00579      &                 +H%R(K)*TEXP%ADR(NTRAC)%P%R(K)*VOLU2D%R(K)
00580         ENDDO
00581         MASSOU(NTRAC)=MASSOU(NTRAC)*DT
00582         IF(NCSIZE.GT.1) MASSOU(NTRAC)=P_DSUM(MASSOU(NTRAC))
00583 !
00584       ENDIF
00585 !
00586 !-----------------------------------------------------------------------
00587 !
00588       RETURN
00589       END

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