majtrac.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\majtrac.f
00002 !
00069                      SUBROUTINE MAJTRAC
00070 !                    ******************
00071 !
00072      &(NS,NT,DIMT,DLIMT,NSEG,NPTFR,NUBO,
00073      & X,Y,AIRS,NU,AIRE,HT,HTN,TN,ZF,NBOR,
00074      & TBOR,FLUTENT,FLUTSOR,SMTR,NORDRE,CMI,JMI,
00075      & DJXT,DJYT,DXT,DYT,DPX,DPY,DIFT,CVIST,BETA,DSZ,AIRST,HSTOK,
00076      & HCSTOK,FLUXT,FLUHBOR,MASSOU,DTT,MESH,
00077      & ELTSEG,IFABOR,VNOCL)
00078 !
00079 !***********************************************************************
00080 ! TELEMAC2D   V6P3                                         27/07/2013
00081 !***********************************************************************
00082 !
00083 !
00084 !
00085 !
00086 !
00087 !
00088 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00089 !| AIRE           |-->| ELEMENT AREA
00090 !| AIRS           |-->| CELL AREA
00091 !| AIRST          |-->| AREA OF SUB-TRIANGLES (SECOND ORDER)
00092 !| BETA           |---| COEFFICIENT OF EXTRAPOLATION FOR ORDRE 2
00093 !| CMI            |-->| COORDINATES OF MIDDLE PONTS OF EDGES
00094 !| CVIST          |-->| COEFFICIENT OF DIFFUSION FOR TRACER
00095 !| DIFT           |-->| LOGICAL TELLING IF THERE IS DIFFUSION FOR TRACER OR NOT
00096 !| DIMT           |-->| DIMENSION OF TRACER
00097 !| DJXT,DJYT      |---| GRADIENTS PER TRIANGLES
00098 !| DLIMT          |-->| DIMENSION OF TRACER AT THE BOUNDARY
00099 !| DSZ            |-->| VARIATION OF Z FOR ORDRE 2
00100 !| DTT            |-->| TIME STEP FOR TRACER
00101 !| DXT,DYT        |---| GRADIENTS AT THE NODES
00102 !| FLUHBOR        |-->| MASS FLUX AT THE BOUNDARY
00103 !| FLUTENT,FLUTSOR|<--| TRACER FLUX AT THE INLET AND OUTLET
00104 !| FLUXT          |-->| MASS FLUX OF TRACER
00105 !| HCSTOK         |-->| STOCKED H RECONSTRUCTED FOR ORDRE 2
00106 !| HSTOK          |-->| STOCKED WATER DEPTH
00107 !| HT             |<--| HT AT TIME N+1
00108 !| HTN,TN         |-->| HT AT TIME N
00109 !| JMI            |-->| NUMBER OF THE TRIANGLE IN WHICH IS LOCATED
00110 !|                |   | THE MIDPOINT OF THE INTERFACE
00111 !| MASSOU         |<--| MASS OF TRACER ADDED BY SOURCE TERM
00112 !| NBOR           |-->| GLOBAL NUMBERING OF BOUNDARY POINTS
00113 !| NORDRE         |-->| ORDRE OF THE SCHEME
00114 !| NPTFR          |-->| TOTAL NUMBER OF BOUNDARY NODES
00115 !| NS             |-->| TOTAL NUMBER OF NODES IN THE MESH
00116 !| NSEG           |-->| TOTAL NUMBER OF EDGES
00117 !| NT             |-->| TOTAL NUMBER OF ELEMENTS
00118 !| NU             |-->| NUMBERING OF NODES IN THE TRIANGLES
00119 !| NUBO           |-->| GLOBAL INDICES OF EDGE EXTREMITIES
00120 !| SMTR           |-->| TRACER SOURCE TERMS
00121 !| TBOR           |-->| BOUNDARY CONDITIONS FOR T
00122 !| X,Y            |-->| COORDINATES IF THE NODES
00123 !| ZF             |-->| BATHYMETRY
00124 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00125 !
00126       USE INTERFACE_TELEMAC2D, EX_MAJTRAC => MAJTRAC
00127       USE BIEF_DEF
00128       USE DECLARATIONS_TELEMAC2D,ONLY: DEBUG
00129 !
00130       IMPLICIT NONE
00131       INTEGER LNG,LU
00132       COMMON/INFO/LNG,LU
00133 !
00134 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00135 !
00136       LOGICAL, INTENT(IN) :: DIFT
00137       INTEGER, INTENT(IN) :: NSEG,NPTFR,NORDRE,DIMT,DLIMT,NS,NT
00138       INTEGER, INTENT(IN) :: NUBO(2,NSEG),NU(NT,3)
00139       INTEGER, INTENT(IN) :: NBOR(NPTFR),JMI(*)
00140       INTEGER, INTENT(IN)             :: ELTSEG(NT,3)
00141       INTEGER, INTENT(IN)             :: IFABOR(NT,3)
00142       DOUBLE PRECISION, INTENT(INOUT) :: HT(DIMT),FLUTENT,FLUTSOR
00143       DOUBLE PRECISION, INTENT(INOUT) :: MASSOU
00144       DOUBLE PRECISION, INTENT(IN)    :: TBOR(DLIMT),DSZ(2,*)
00145       DOUBLE PRECISION, INTENT(IN)    :: X(NS),Y(NS),AIRS(NS),AIRE(NT)
00146       DOUBLE PRECISION, INTENT(IN)    :: HTN(DIMT),TN(DIMT),ZF(*)
00147       DOUBLE PRECISION, INTENT(IN)    :: SMTR(DIMT),DPX(3,NT),DPY(3,NT)
00148       DOUBLE PRECISION, INTENT(IN)    :: CMI(2,*),AIRST(2,*),CVIST
00149       DOUBLE PRECISION, INTENT(INOUT) :: DJXT(*),DJYT(*),DXT(*),DYT(*)
00150       DOUBLE PRECISION, INTENT(INOUT) :: BETA
00151       DOUBLE PRECISION, INTENT(IN)    :: HSTOK(*),VNOCL(3,NSEG)
00152       DOUBLE PRECISION, INTENT(IN)    :: HCSTOK(2,*),FLUXT(*)
00153       DOUBLE PRECISION, INTENT(IN)    :: FLUHBOR(*),DTT
00154       TYPE(BIEF_MESH), INTENT(INOUT)  :: MESH
00155 !
00156 !
00157 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00158 !
00159       INTEGER IS,K,NSG,NUBO1,NUBO2,J,ILIM,ERR,I,IEL
00160 !
00161       DOUBLE PRECISION ZF1,ZF2,FLUH,FLUT,HI0,HJ0,AIX,AIY,AJX,AJY,AMDS
00162       DOUBLE PRECISION FLU11,FLU41,UAS41,UAS42,DSZ1,DSZ2
00163       DOUBLE PRECISION DEMI,PROD_SCAL
00164 !
00165 !     DYNAMIC ARRAY ALLOCATION !!!!!!!!
00166 !
00167       DOUBLE PRECISION, ALLOCATABLE, SAVE :: CET(:),DST(:,:)
00168       DOUBLE PRECISION, ALLOCATABLE, SAVE :: DSP(:),DSM(:),CORRT(:)
00169 !
00170       DOUBLE PRECISION,ALLOCATABLE,SAVE  :: GRADI(:),GRADJ(:)
00171       DOUBLE PRECISION,ALLOCATABLE,SAVE  :: GRADIJ(:),GRADJI(:)
00172       LOGICAL, ALLOCATABLE ::   YESNO(:)
00173 !
00174       LOGICAL DEJA
00175       DATA DEJA/.FALSE./
00176 !
00177 !-----------------------------------------------------------------------
00178 !
00179       IF(.NOT.DEJA) THEN
00180         ALLOCATE(CET(NS),STAT=ERR)
00181         IF(ERR.NE.0) GO TO 1001
00182         ALLOCATE(DST(2,NSEG),STAT=ERR)
00183         IF(ERR.NE.0) GO TO 1001
00184         ALLOCATE(DSP(NS),STAT=ERR)
00185         IF(ERR.NE.0) GO TO 1001
00186         ALLOCATE(DSM(NS)    ,STAT=ERR)
00187         IF(ERR.NE.0) GO TO 1001
00188         ALLOCATE(CORRT(NS)   ,STAT=ERR)
00189         IF(ERR.NE.0) GO TO 1001
00190         ALLOCATE(GRADI(NSEG),GRADJ(NSEG),STAT=ERR)
00191         IF(ERR.NE.0) GO TO 1001
00192         ALLOCATE(GRADIJ(NSEG),GRADJI(NSEG),STAT=ERR)
00193         IF(ERR.NE.0) GO TO 1001
00194         GO TO 1002
00195 1001    CONTINUE
00196         IF(LNG.EQ.1) WRITE(LU,1000) ERR
00197         IF(LNG.EQ.2) WRITE(LU,2000) ERR
00198 1000    FORMAT(1X,'MAJTRAC : ERREUR A L''ALLOCATION DE MEMOIRE : ',/,1X,
00199      &         'CODE D''ERREUR : ',1I6)
00200 2000    FORMAT(1X,'MAFTRAC: ERROR DURING ALLOCATION OF MEMORY: ',/,1X,
00201      &        'ERROR CODE: ',1I6)
00202         CALL PLANTE(1)
00203         STOP
00204 1002    CONTINUE
00205         DEJA=.TRUE.
00206       ENDIF
00207       DEMI = 0.5D0
00208 !
00209 !-----------------------------------------------------------------------
00210 !
00211 !  INITIALISES
00212 !
00213       CET(:)=(/(0.D0,IS=1,NS)/)
00214       ALLOCATE(YESNO(NSEG),STAT=ERR)
00215       IF(ERR.GT.0)THEN
00216         IF(LNG.EQ.1) WRITE(LU,1000) ERR
00217         IF(LNG.EQ.2) WRITE(LU,2000) ERR
00218         CALL PLANTE(1)
00219         STOP
00220       ENDIF
00221 !
00222 ! INITIALIZATION OF YESNO
00223       DO I=1,NSEG
00224         YESNO(I)=.FALSE.
00225       ENDDO
00226 !
00227 !   COMPUTES THE TRACER GRADIENTS BY TRIANGLE AND BY NODE
00228 !   COMPUTES THE DIFFUSION TERM
00229 !
00230       IF(DIFT.OR.NORDRE.EQ.2) CALL GRADNODT(NS,NT,NU,AIRE,AIRS,
00231      &HSTOK,TN,DPX,DPY,DJXT,DJYT,DXT,DYT,DIFT,CVIST,CET,DTT,MESH)
00232 !
00233       IF(NORDRE.EQ.2) THEN
00234 !
00235 !  REBUILDS 2ND ORDER FOR TRACER
00236 !  *************************************
00237 !
00238 !    INITIALIZATION
00239       DSP(:)  =(/(0.D0,IS=1,NS)/)
00240       DSM(:)  =(/(0.D0,IS=1,NS)/)
00241       DST(1,:)=(/(0.D0,IS=1,NSEG)/)
00242       DST(2,:)=(/(0.D0,IS=1,NSEG)/)
00243 !    INITIALIZATION  OF GRADIENTS
00244       GRADI(:) =(/(0.D0,IS=1,NSEG)/)
00245       GRADJ(:) =(/(0.D0,IS=1,NSEG)/)
00246       GRADIJ(:)=(/(0.D0,IS=1,NSEG)/)
00247       GRADJI(:)=(/(0.D0,IS=1,NSEG)/)
00248 !
00249 !
00250       DO IEL=1, NT
00251         DO I = 1,3
00252           IF(.NOT.YESNO(ELTSEG(IEL,I)))THEN
00253             NSG = ELTSEG(IEL,I)
00254 !     RECUPERATE JMI
00255             J   = JMI(NSG) ! THIS THE TRIANGLE IN WHICH IS LOCATED CMI
00256             IF(NCSIZE.GT.1.AND.J.EQ.0) CYCLE  ! THAT MEANS CMI IS NOT LOCATED IN TRIANGLE J
00257 !
00258 !    RECUPERATE NODES OF THE EDGE WITH THE GOOD ORIENTATION
00259 !     WITH RESPECT TO THE NORMAL
00260             NUBO1 = NUBO(1,NSG)
00261             NUBO2 = NUBO(2,NSG)
00262             PROD_SCAL= ((X(NUBO2)-X(NUBO1))*VNOCL(1,NSG)+
00263      &                  (Y(NUBO2)-Y(NUBO1))*VNOCL(2,NSG))
00264             IF(PROD_SCAL.LT.0.D0)THEN
00265               NUBO1 = NUBO(2,NSG)
00266               NUBO2 = NUBO(1,NSG)
00267             ENDIF
00268 !
00269             ZF1 = ZF(NUBO1)
00270             ZF2 = ZF(NUBO2)
00271 !
00272             IF(PROD_SCAL.LT.0.D0)THEN
00273               DSZ1 = DSZ(2,NSG)
00274               DSZ2 = DSZ(1,NSG)
00275             ELSE
00276               DSZ1 = DSZ(1,NSG)
00277               DSZ2 = DSZ(2,NSG)
00278             ENDIF
00279 !
00280             HI0 = HSTOK(NUBO1)
00281             HJ0 = HSTOK(NUBO2)
00282 !
00283 !           STICKS TO 1ST ORDER FOR A COVERED EDGE
00284 !
00285             IF(ZF1.GE. (HJ0+ZF2) .OR. ZF2.GE. (HI0+ZF1)
00286      &         .OR. 2.*ABS(DSZ1).GE.HI0
00287      &         .OR. 2.*ABS(DSZ1).GE.HJ0
00288      &         .OR. 2.*ABS(DSZ2).GE.HI0
00289      &         .OR. 2.*ABS(DSZ2).GE.HJ0)  THEN
00290 !             DST(1,NSG) =0.D0
00291 !             DST(2,NSG) =0.D0
00292               CYCLE
00293             ELSE
00294 !
00295               AIX         = CMI(1,NSG)-X(NUBO1)
00296               AIY         = CMI(2,NSG)-Y(NUBO1)
00297               AJX         = CMI(1,NSG)-X(NUBO2)
00298               AJY         = CMI(2,NSG)-Y(NUBO2)
00299 !
00300               GRADI(NSG)  = AIX*DXT(NUBO1) + AIY*DYT(NUBO1)
00301               GRADJ(NSG)  = AJX*DXT(NUBO2) + AJY*DYT(NUBO2)
00302               GRADIJ(NSG) = AIX*DJXT(J) + AIY*DJYT(J)
00303               GRADJI(NSG) = AJX*DJXT(J) + AJY*DJYT(J)
00304             ENDIF
00305             YESNO(NSG)=.TRUE.
00306           ENDIF
00307         ENDDO
00308       ENDDO
00309       IF(NCSIZE.GT.1)THEN      ! NPON,NPLAN,ICOM,IAN , HERE ICOM=1 VALUE WITH MAX | |
00310         CALL PARCOM2_SEG(GRADI,GRADJ,GRADI,
00311      &              NSEG,1,2,2,MESH,1,11)
00312         CALL PARCOM2_SEG(GRADIJ,GRADJI,GRADJI,
00313      &              NSEG,1,2,2,MESH,1,11)
00314       ENDIF
00315 !
00316 !    EXTRAPOLATES THE GRADIENTS AND SLOPE LIMITOR
00317 !
00318       ILIM=2
00319       BETA=0.3333D0
00320       DO NSG=1,NSEG
00321         DST(1,NSG)  = EXLIM (ILIM,BETA,GRADI(NSG),GRADIJ(NSG))
00322         DST(2,NSG)  = EXLIM (ILIM,BETA,GRADJ(NSG),GRADJI(NSG))
00323       ENDDO
00324 !
00325 ! INITIALIZATION OF YESNO
00326       DO I=1,NSEG
00327         YESNO(I)=.FALSE.
00328       ENDDO
00329       DO IEL=1, NT
00330         DO I = 1,3
00331           IF(.NOT.YESNO(ELTSEG(IEL,I)))THEN
00332             NSG = ELTSEG(IEL,I)
00333 !    RECUPERATE NODES OF THE EDGE WITH THE GOOD ORIENTATION
00334 !     WITH RESPECT TO THE NORMAL
00335             NUBO1 = NUBO(1,NSG)
00336             NUBO2 = NUBO(2,NSG)
00337             PROD_SCAL= ((X(NUBO2)-X(NUBO1))*VNOCL(1,NSG)+
00338      &                  (Y(NUBO2)-Y(NUBO1))*VNOCL(2,NSG))
00339             IF(PROD_SCAL.LT.0.D0)THEN
00340               NUBO1 = NUBO(2,NSG)
00341               NUBO2 = NUBO(1,NSG)
00342             ENDIF
00343 !
00344 !           FOR PARALLELILSM
00345             IF(NCSIZE.GT.1.AND.IFABOR(IEL,I).EQ.-2)THEN ! THIS IS AN INTERFACE EDGE
00346               IF(DST(1,NSG).GE.0.D0) THEN
00347                 DSP(NUBO1) = DSP(NUBO1) +
00348      &                      DEMI*AIRST(1,NSG)*HCSTOK(1,NSG)*DST(1,NSG) ! WE CONSIDER ONLY
00349               ELSE                                                      ! 0.5 AIRST
00350                 DSM(NUBO1) = DSM(NUBO1) -
00351      &                      DEMI*AIRST(1,NSG)*HCSTOK(1,NSG)*DST(1,NSG) ! PARCOM2 WILL ADD
00352               ENDIF                                                     ! CONTRIBUTIONS
00353               IF(DST(2,NSG).GE.0.D0) THEN
00354                 DSP(NUBO2) = DSP(NUBO2) +
00355      &                      DEMI*AIRST(2,NSG)*HCSTOK(2,NSG)*DST(2,NSG)
00356               ELSE
00357                 DSM(NUBO2) = DSM(NUBO2) -
00358      &                      DEMI*AIRST(2,NSG)*HCSTOK(2,NSG)*DST(2,NSG)
00359               ENDIF
00360             ELSE ! NO PARALLELILSM OR NO INTERFACE EDGE
00361               IF(DST(1,NSG).GE.0.D0) THEN
00362                 DSP(NUBO1) = DSP(NUBO1) +
00363      &          AIRST(1,NSG)* HCSTOK(1,NSG)*DST(1,NSG)
00364               ELSE
00365                 DSM(NUBO1) = DSM(NUBO1) -
00366      &          AIRST(1,NSG)* HCSTOK(1,NSG)*DST(1,NSG)
00367               ENDIF
00368               IF(DST(2,NSG).GE.0.) THEN
00369                 DSP(NUBO2) = DSP(NUBO2) +
00370      &          AIRST(2,NSG)* HCSTOK(2,NSG)*DST(2,NSG)
00371               ELSE
00372                 DSM(NUBO2) = DSM(NUBO2) -
00373      &          AIRST(2,NSG)* HCSTOK(2,NSG)*DST(2,NSG)
00374               ENDIF
00375             ENDIF
00376 !
00377             YESNO(NSG)=.TRUE.
00378           ENDIF
00379         ENDDO
00380       ENDDO
00381       !  FOR PARALLELILSM
00382       IF(NCSIZE.GT.1)THEN
00383         CALL PARCOM2(DSP,DSM,DSM,NS,1,2,2,MESH)
00384       ENDIF
00385 !
00386 !     COMPUTES THE CORRECTIONS TO ENSURE CONSERVATION OF HT
00387 !                  ***********           ******************
00388 !
00389       DO IS=1,NS
00390         CORRT(IS) =  DSM(IS) - DSP(IS)
00391         AMDS =MAX(DSP(IS),DSM(IS))
00392         IF(AMDS.GT.0.D0) THEN
00393           CORRT(IS) = CORRT(IS)/AMDS
00394         ENDIF
00395       ENDDO
00396 !
00397       ENDIF ! ENDIF OF NORDRE.EQ.2
00398 !
00399 !     COMPUTES FLUXES FOR THE INTERNAL INTERFACES
00400 !
00401 !  REINITIALIZATION OF YESNO
00402       DO I=1,NSEG
00403         YESNO(I)=.FALSE.
00404       ENDDO
00405 !
00406 !     LOOP ON GLOBAL LIST OF EDGES
00407 !    ******************************
00408 !
00409       DO IEL=1, NT
00410         DO I = 1,3
00411           IF(.NOT.YESNO(ELTSEG(IEL,I)))THEN
00412             NSG = ELTSEG(IEL,I)
00413 !
00414 !    RECUPERATE NODES OF THE EDGE WITH THE GOOD ORIENTATION
00415 !     WITH RESPECT TO THE NORMAL
00416             NUBO1 = NUBO(1,NSG)
00417             NUBO2 = NUBO(2,NSG)
00418             PROD_SCAL= ((X(NUBO2)-X(NUBO1))*VNOCL(1,NSG)+
00419      &                  (Y(NUBO2)-Y(NUBO1))*VNOCL(2,NSG))
00420             IF(PROD_SCAL.LT.0.D0)THEN
00421               NUBO1 = NUBO(2,NSG)
00422               NUBO2 = NUBO(1,NSG)
00423             ENDIF
00424 !
00425             UAS41     = TN(NUBO1)
00426             UAS42     = TN(NUBO2)
00427 !
00428             FLU11=FLUXT(NSG)
00429 !
00430             IF (FLU11.GE.0.) THEN
00431               IF(NORDRE.GE.2) THEN
00432                 UAS41 = UAS41  + DST(1,NSG) +
00433      &          MIN(0.D0,CORRT(NUBO1))*MAX(0.D0,DST(1,NSG))+
00434      &          MAX(0.D0,CORRT(NUBO1))*MAX(0.D0,-DST(1,NSG))
00435               ENDIF
00436               FLU41 =  UAS41 * FLU11
00437             ELSE
00438               IF(NORDRE.GE.2) THEN
00439                 UAS42 = UAS42 + DST(2,NSG) +
00440      &          MIN(0.D0,CORRT(NUBO2))*MAX(0.D0,DST(2,NSG))+
00441      &          MAX(0.D0,CORRT(NUBO2))*MAX(0.D0,-DST(2,NSG))
00442               ENDIF
00443               FLU41 =  UAS42 * FLU11
00444             ENDIF
00445 !
00446             CET(NUBO1) = CET(NUBO1) - FLU41
00447             CET(NUBO2) = CET(NUBO2) + FLU41
00448 !
00449             YESNO(NSG)=.TRUE.
00450           ENDIF
00451         ENDDO
00452       ENDDO
00453       !  FOR PARALLELILSM
00454       IF(NCSIZE.GT.1)THEN
00455         CALL PARCOM2(CET,CET,CET,NS,1,2,1,MESH)
00456       ENDIF
00457 !
00458 !     BOUNDARY FLUX
00459 !
00460       IF(NPTFR.GT.0)THEN  ! USEFUL FOR PARALLEL CASE
00461         DO K=1,NPTFR
00462           IS=NBOR(K)
00463 !
00464           FLUH =FLUHBOR(K)
00465 !
00466           IF(FLUH.GE.0.D0) THEN
00467             FLUT= TN(IS)*FLUH
00468             FLUTSOR = FLUTSOR +FLUT
00469           ELSE
00470             FLUT= TBOR(K)*FLUH
00471             FLUTENT = FLUTENT +FLUT
00472           ENDIF
00473 !
00474           CET(IS)  = CET(IS) - FLUT
00475 !
00476         ENDDO
00477       ENDIF
00478 !
00479 !     UPDATES HT
00480 !
00481       DO IS =1,NS
00482 !
00483         HT(IS)  = HTN(IS) +  (CET(IS)+SMTR(IS))/AIRS(IS)
00484         MASSOU = MASSOU + SMTR(IS)
00485 !
00486         IF(HT(IS).LE.1.D-15) HT(IS)=0.D0
00487 !
00488       ENDDO
00489 !
00490       DEALLOCATE(YESNO)
00491 !
00492 !-----------------------------------------------------------------------
00493 !
00494       RETURN
00495       END

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