flucin.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\flucin.f
00002 !
00066                      SUBROUTINE FLUCIN
00067 !                    *****************
00068 !
00069      &(NS,NELEM,NSEG,NUBO,G,X,Y,CFL,DT,UA,ZF,VNOCL,CE,NORDRE,CMI,JMI,
00070      & DJX,DJY,DX,DY,BETA,DSZ0,AIRS,AIRST,HC,FLUXTEMP,NPTFR,NBOR,
00071      & XNEBOR,YNEBOR,NTRAC,ELTSEG,IFABOR,MESH)
00072 !
00073 !***********************************************************************
00074 ! TELEMAC2D   V6P3                                   21/07/2013
00075 !***********************************************************************
00076 !
00077 !
00078 !
00079 !
00080 !
00081 !
00082 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00083 !| AIRS           |-->| CELL'S AREA
00084 !| AIRST          |-->| AREA OF SUB-TRIANGLES IN THE CELLS
00085 !| BETA           |-->| EXTRAPOLATION COEFFICIENT FOR ORDRE 2
00086 !| CE             |<->| FLUX INCREMENTS AT INTERNAL INTERFACES
00087 !| CFL            |-->| CFL NUMBER
00088 !| CMI            |-->| COORDINATES OF MIDDLE POINTS OF INTERFACES
00089 !| DJX,DJY        |-->| GRADIENTS PER TRIANGLE
00090 !| DSZ0           |-->| VARIATIONS Z (BATHY) FOR ORDRE 2
00091 !| DT             |<->| TIME STEP (CAN CHANGE IF ORDRE 2)
00092 !| DX,DY          |-->| GRADIENTS PER NODE
00093 !| FLUXTEMP       |<--| MASS FLUX OF TRACER
00094 !| G              |-->| GRAVITY
00095 !| HC             |<--| REBUILT H (ORDRE 2)
00096 !| JMI            |-->| NUMBER OF THE TRIANGLE IN WHICH IS LOCATED
00097 !|                |   | THE MIDDLE POINT OF THE INTERFACE
00098 !| NBOR           |-->| GLOBAL NUMBER OF BOUNDARY NODES
00099 !| NORDRE         |-->| ORDRE OF THE SCHEME
00100 !| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
00101 !| NS             |-->| TOTAL NUMER OF POINTS IN THE MESH
00102 !| NSEG           |-->| NUMBER OF EDGES IN THE MESH
00103 !| NTRAC          |-->| NUMBER OF TRACERS
00104 !| NUBO           |-->| GLOBAL NUMBERS OF THE NODES FORMING THE EDGE
00105 !| UA             |-->| UA(1,IS) = H,  UA(2,IS)=U  ,UA(3,IS)=V
00106 !| VNOCL          |-->| NORMAL VECTOR TO THE INTERFACE
00107 !|                |   | (2 FIRST COMPONENTS) AND
00108 !|                |   | LENGTH OF THE SEGMENT (3RD COMPONENT)
00109 !| X,Y            |-->| COORDINATES IF THE NODES
00110 !| XNEBOR,YNEBOR  |-->| NORMAL VECTOR TO BOUNDARY NODES
00111 !| ZF             |-->| BATHYMETRY
00112 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00113 !
00114       USE BIEF
00115       USE DECLARATIONS_TELEMAC2D , ONLY:DEBUG
00116 !
00117       IMPLICIT NONE
00118       INTEGER LNG,LU
00119       COMMON/INFO/LNG,LU
00120 !
00121 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00122 !
00123       INTEGER, INTENT(IN) :: NS,NSEG,NPTFR,NORDRE,NTRAC,NELEM
00124       INTEGER, INTENT(IN) :: NBOR(*),NUBO(2,NSEG),JMI(*)
00125       INTEGER, INTENT(IN)             :: ELTSEG(NELEM,3)
00126       DOUBLE PRECISION, INTENT(IN)    :: XNEBOR(*),YNEBOR(*),X(NS),Y(NS)
00127       DOUBLE PRECISION, INTENT(IN)    :: ZF(NS),VNOCL(3,NSEG),AIRS(*)
00128       DOUBLE PRECISION, INTENT(IN)    :: G,CFL,UA(3,NS),AIRST(2,*)
00129       DOUBLE PRECISION, INTENT(IN)    :: DSZ0(2,*),CMI(2,*)
00130       DOUBLE PRECISION, INTENT(INOUT) :: BETA,DT,CE(NS,3),HC(2,*)
00131       DOUBLE PRECISION, INTENT(IN)    :: DJX(3,*),DJY(3,*)
00132       DOUBLE PRECISION, INTENT(IN)    :: DX(3,*),DY(3,*)
00133       TYPE(BIEF_OBJ), INTENT(INOUT)   :: FLUXTEMP
00134       INTEGER, INTENT(IN)             :: IFABOR(NELEM,3)
00135       TYPE(BIEF_MESH),INTENT(INOUT)   :: MESH
00136 !
00137 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00138 !
00139 !     AUTOMATIC EX ARRAYS!!!!!!
00140 !
00141       DOUBLE PRECISION, ALLOCATABLE,SAVE :: DSH(:,:),DSU(:,:)
00142       DOUBLE PRECISION, ALLOCATABLE,SAVE :: DSV(:,:)
00143       DOUBLE PRECISION, ALLOCATABLE,SAVE :: DSP(:),DSM(:),DSZ(:,:)
00144       DOUBLE PRECISION, ALLOCATABLE,SAVE :: CORR(:),DTLL(:)
00145 ! ra29/042013      DOUBLE PRECISION GRADI(3),GRADJ(3),GRADIJ(3),GRADJI(3)
00146       DOUBLE PRECISION,ALLOCATABLE,SAVE  :: GRADI(:,:),GRADJ(:,:)
00147       DOUBLE PRECISION,ALLOCATABLE,SAVE  :: GRADIJ(:,:),GRADJI(:,:)
00148 !end ra
00149 !
00150       LOGICAL DEJA
00151       DATA DEJA/.FALSE./
00152 !
00153 !-----------------------------------------------------------------------
00154 !
00155       INTEGER NSG,NUBO1,NUBO2,J,IVAR,IS,K,ILIM,ERR,ITRAC,I,IEL
00156 !
00157       DOUBLE PRECISION VNX,VNY,VNL,RA2,RA3,ALP,ZF1,ZF2,XNN,YNN,RNN
00158       DOUBLE PRECISION UAS11,UAS12,UAS21,UAS22,UAS31,UAS32,AMDS
00159       DOUBLE PRECISION GRADI2,GRADIJ2,GRADJ2,GRADJI2,DEMI
00160       DOUBLE PRECISION AIX,AIY,AJX,AJY,HI,HI0,HIJ,CIJ,UAS210,UAS220
00161       DOUBLE PRECISION HJ,HJ0,HJI,CJI,DZIJ,DZJI
00162       DOUBLE PRECISION EXT1,EXT2,FLU11,FLU21,FLU31,FLU12,FLU22
00163       DOUBLE PRECISION FLU210,A01,A11,A21,A02,A12,A22
00164       DOUBLE PRECISION HGZI,HGZJ,HDXZ1,HDYZ1,HDXZ2,HDYZ2
00165       DOUBLE PRECISION SIGMAX,DTL,UNORM,DSZ1,DSZ2,AUX,AUX1
00166       DOUBLE PRECISION PROD_SCAL
00167       LOGICAL, ALLOCATABLE ::   YESNO(:)
00168 !
00169       DOUBLE PRECISION EXLIM,P_DMIN
00170       EXTERNAL         EXLIM,P_DMIN
00171 !
00172 !-----------------------------------------------------------------------
00173 !
00174       IF(.NOT.DEJA) THEN
00175         ALLOCATE(DSH(2,NSEG),STAT=ERR)
00176         IF(ERR.NE.0) GO TO 1001
00177         ALLOCATE(DSU(2,NSEG),STAT=ERR)
00178         IF(ERR.NE.0) GO TO 1001
00179         ALLOCATE(DSV(2,NSEG),STAT=ERR)
00180         IF(ERR.NE.0) GO TO 1001
00181         ALLOCATE(DSP(NS)    ,STAT=ERR)
00182         IF(ERR.NE.0) GO TO 1001
00183         ALLOCATE(DSM(NS)    ,STAT=ERR)
00184         IF(ERR.NE.0) GO TO 1001
00185         ALLOCATE(DSZ(2,NSEG),STAT=ERR)
00186         IF(ERR.NE.0) GO TO 1001
00187         ALLOCATE(CORR(NS)   ,STAT=ERR)
00188         IF(ERR.NE.0) GO TO 1001
00189         ALLOCATE(DTLL(NS)   ,STAT=ERR)
00190         IF(ERR.NE.0) GO TO 1001
00191         ALLOCATE(GRADI(3,NSEG),GRADJ(3,NSEG),STAT=ERR)
00192         IF(ERR.NE.0) GO TO 1001
00193         ALLOCATE(GRADIJ(3,NSEG),GRADJI(3,NSEG),STAT=ERR)
00194         IF(ERR.NE.0) GO TO 1001
00195         GO TO 1002
00196 1001    CONTINUE
00197         IF(LNG.EQ.1) WRITE(LU,1000) ERR
00198         IF(LNG.EQ.2) WRITE(LU,2000) ERR
00199 1000    FORMAT(1X,'FLUCIN : ERREUR A L''ALLOCATION DE MEMOIRE : ',/,1X,
00200      &         'CODE D''ERREUR : ',1I6)
00201 2000    FORMAT(1X,'FLUCIN: ERROR DURING ALLOCATION OF MEMORY: ',/,1X,
00202      &        'ERROR CODE: ',1I6)
00203         CALL PLANTE(1)
00204         STOP
00205 1002    CONTINUE
00206         DEJA=.TRUE.
00207       ENDIF
00208       DEMI = 0.5D0
00209 !
00210 !-----------------------------------------------------------------------
00211 !
00212       RA2 = SQRT(2.D0)
00213       RA3 = SQRT(1.5D0*G)
00214       ALP = 0.5D0/RA3
00215       ALLOCATE(YESNO(NSEG),STAT=ERR)
00216       IF(ERR.NE.0)THEN
00217         IF(LNG.EQ.1) WRITE(LU,1000) ERR
00218         IF(LNG.EQ.2) WRITE(LU,2000) ERR
00219         CALL PLANTE(1)
00220         STOP
00221       ENDIF
00222 !
00223 ! INITIALIZATION OF YESNO
00224       DO I=1,NSEG
00225         YESNO(I)=.FALSE.
00226       ENDDO
00227 !
00228       IF(NORDRE.EQ.1) GOTO 12
00229 !
00230 !  2ND ORDER RECONSTRUCTION
00231 !  ************************
00232 !
00233 !    INITIALIZATION
00234       DTLL(:)=(/(1.E10,I=1,NS)/)
00235       DSP (:)=(/(0.D0,I=1,NS)/)
00236       DSM (:)=(/(0.D0,I=1,NS)/)
00237       DSH(1,:)=(/(0.D0,I=1,NSEG)/)
00238       DSH(2,:)=(/(0.D0,I=1,NSEG)/)
00239       DSU(1,:)=(/(0.D0,I=1,NSEG)/)
00240       DSU(2,:)=(/(0.D0,I=1,NSEG)/)
00241       DSV(1,:)=(/(0.D0,I=1,NSEG)/)
00242       DSV(2,:)=(/(0.D0,I=1,NSEG)/)
00243       DSZ(1,:)=(/(0.D0,I=1,NSEG)/)
00244       DSZ(2,:)=(/(0.D0,I=1,NSEG)/)
00245 !GIVES ERROR WITH INTEL COMPILER IF WRITTEN LIKE THIS
00246 !       CALL OV( 'X=C     ' ,DSH(1,1:NSEG),DSM ,DSM ,0.D0,NSEG)
00247 !
00248 !    INITIALIZATION  OF GRADIENTS
00249       GRADI(1,:)=(/(0.D0,I=1,NSEG)/)
00250       GRADI(2,:)=(/(0.D0,I=1,NSEG)/)
00251       GRADI(3,:)=(/(0.D0,I=1,NSEG)/)
00252 !
00253       GRADJ(1,:)=(/(0.D0,I=1,NSEG)/)
00254       GRADJ(2,:)=(/(0.D0,I=1,NSEG)/)
00255       GRADJ(3,:)=(/(0.D0,I=1,NSEG)/)
00256 !
00257       GRADIJ(1,:)=(/(0.D0,I=1,NSEG)/)
00258       GRADIJ(2,:)=(/(0.D0,I=1,NSEG)/)
00259       GRADIJ(3,:)=(/(0.D0,I=1,NSEG)/)
00260 !
00261       GRADJI(1,:)=(/(0.D0,I=1,NSEG)/)
00262       GRADJI(2,:)=(/(0.D0,I=1,NSEG)/)
00263       GRADJI(3,:)=(/(0.D0,I=1,NSEG)/)
00264 !
00265       DO IEL=1, NELEM
00266         DO I = 1,3
00267           IF(.NOT.YESNO(ELTSEG(IEL,I)))THEN
00268             NSG = ELTSEG(IEL,I)
00269 !     RECUPERATE JMI
00270             J   = JMI(NSG) ! THIS THE TRIANGLE IN WHICH IS LOCATED CMI
00271             IF(NCSIZE.GT.1.AND.J.EQ.0)CYCLE  ! THAT MEANS CMI IS NOT LOCATED IN TRIANGLE J
00272 !
00273 !    RECUPERATE NODES OF THE EDGE WITH THE GOOD ORIENTATION
00274 !     WITH RESPECT TO THE NORMAL
00275             NUBO1 = NUBO(1,NSG)
00276             NUBO2 = NUBO(2,NSG)
00277             PROD_SCAL= ((X(NUBO2)-X(NUBO1))*VNOCL(1,NSG)+
00278      &                  (Y(NUBO2)-Y(NUBO1))*VNOCL(2,NSG))
00279             IF(PROD_SCAL.LT.0.D0)THEN
00280               NUBO1 = NUBO(2,NSG)
00281               NUBO2 = NUBO(1,NSG)
00282             ENDIF
00283 !
00284             ZF1        = ZF(NUBO1)
00285             ZF2        = ZF(NUBO2)
00286             IF(PROD_SCAL.LT.0.D0)THEN
00287               DSZ(1,NSG) = DSZ0(2,NSG)
00288               DSZ(2,NSG) = DSZ0(1,NSG)
00289             ELSE
00290               DSZ(1,NSG) = DSZ0(1,NSG)
00291               DSZ(2,NSG) = DSZ0(2,NSG)
00292             ENDIF
00293 !
00294             HI0=UA(1,NUBO1)
00295             HJ0=UA(1,NUBO2)
00296 !
00297 !   FOR AN EDGE BEING RECOVERED, ONE WILL REMAIN 1ST ORDER
00298 !
00299             IF(ZF1.GE. (HJ0+ZF2) .OR. ZF2.GE. (HI0+ZF1)
00300      &     .OR. 2.*ABS(DSZ(1,NSG)).GE.HI0
00301      &     .OR. 2.*ABS(DSZ(1,NSG)).GE.HJ0
00302      &     .OR. 2.*ABS(DSZ(2,NSG)).GE.HI0
00303      &     .OR. 2.*ABS(DSZ(2,NSG)).GE.HJ0)  THEN
00304 !ra02/05/2013 FOR OPTIMIZATION
00305               CYCLE
00306 !           DSH(1,NSG) =0.D0
00307 !           DSH(2,NSG) =0.D0
00308 !           DSU(1,NSG) =0.D0
00309 !           DSU(2,NSG) =0.D0
00310 !           DSV(1,NSG) =0.D0
00311 !           DSV(2,NSG) =0.D0
00312 !           DSZ(1,NSG) =0.D0
00313 !           DSZ(2,NSG) =0.D0
00314             ELSE
00315 !
00316 !     NORMALIZED UNIT NORMAL (VNOCL), RNN LENGTH OF LIJ
00317 !
00318               XNN = VNOCL(1,NSG)
00319               YNN = VNOCL(2,NSG)
00320               RNN = VNOCL(3,NSG)
00321 !
00322               AIX = CMI(1,NSG)-X(NUBO1) ! THESE ARE COORDINATES OF
00323               AIY = CMI(2,NSG)-Y(NUBO1) !  VECTOR PM (EQ 5.1)
00324               AJX = CMI(1,NSG)-X(NUBO2) ! P: NUBO1 OR NUBO2
00325               AJY = CMI(2,NSG)-Y(NUBO2) ! M: CMI(NSG)
00326 !
00327               DO IVAR=1,3
00328                 GRADI(IVAR,NSG) = AIX*DX(IVAR,NUBO1)+AIY*DY(IVAR,NUBO1)!NODE GRADIENT (PM.GRADZ)
00329                 GRADJ(IVAR,NSG) = AJX*DX(IVAR,NUBO2)+AJY*DY(IVAR,NUBO2)!eq 5.1 of audusse paper)
00330 !
00331                 GRADIJ(IVAR,NSG) = AIX*DJX(IVAR,J) + AIY*DJY(IVAR,J)
00332                 GRADJI(IVAR,NSG) = AJX*DJX(IVAR,J) + AJY*DJY(IVAR,J)
00333               ENDDO
00334 ! ROTATION OF THE GRADIENTS
00335 !
00336               GRADI2       = GRADI(2,NSG)
00337               GRADI(2,NSG) = XNN*GRADI2+YNN*GRADI(3,NSG)
00338               GRADI(3,NSG) =-YNN*GRADI2+XNN*GRADI(3,NSG)
00339 !
00340               GRADIJ2      = GRADIJ(2,NSG)
00341               GRADIJ(2,NSG)= XNN*GRADIJ2+YNN*GRADIJ(3,NSG)
00342               GRADIJ(3,NSG)=-YNN*GRADIJ2+XNN*GRADIJ(3,NSG)
00343 !
00344               GRADJ2       = GRADJ(2,NSG)
00345               GRADJ(2,NSG) = XNN*GRADJ2+YNN*GRADJ(3,NSG)
00346               GRADJ(3,NSG) =-YNN*GRADJ2+XNN*GRADJ(3,NSG)
00347 !
00348               GRADJI2      = GRADJI(2,NSG)
00349               GRADJI(2,NSG)= XNN*GRADJI2+YNN*GRADJI(3,NSG)
00350               GRADJI(3,NSG)=-YNN*GRADJI2+XNN*GRADJI(3,NSG)
00351 !
00352             ENDIF
00353             YESNO(NSG)=.TRUE.
00354           ENDIF
00355         ENDDO
00356       ENDDO
00357       IF(NCSIZE.GT.1)THEN      ! NPON,NPLAN,ICOM,IAN , HERE ICOM=1 VALUE WITH MAX | |
00358         CALL PARCOM2_SEG(GRADI(1,1:NSEG),
00359      &                   GRADI(2,1:NSEG),
00360      &                   GRADI(3,1:NSEG),
00361      &              NSEG,1,2,3,MESH,1,11)
00362         CALL PARCOM2_SEG(GRADJ(1,1:NSEG),
00363      &                   GRADJ(2,1:NSEG),
00364      &                   GRADJ(3,1:NSEG),
00365      &              NSEG,1,2,3,MESH,1,11)
00366         CALL PARCOM2_SEG(GRADIJ(1,1:NSEG),
00367      &                   GRADIJ(2,1:NSEG),
00368      &                   GRADIJ(3,1:NSEG),
00369      &              NSEG,1,2,3,MESH,1,11)
00370         CALL PARCOM2_SEG(GRADJI(1,1:NSEG),
00371      &                   GRADJI(2,1:NSEG),
00372      &                   GRADJI(3,1:NSEG),
00373      &              NSEG,1,2,3,MESH,1,11)
00374       ENDIF
00375 !
00376 !    EXTRAPOLATES THE GRADIENTS AND USES OF SLOPE LIMITERS
00377 !
00378 !
00379 ! INITIALIZATION OF YESNO
00380       DO I=1,NSEG
00381         YESNO(I)=.FALSE.
00382       ENDDO
00383       DO IEL=1, NELEM
00384         DO I = 1,3
00385           IF(.NOT.YESNO(ELTSEG(IEL,I)))THEN
00386             NSG = ELTSEG(IEL,I)
00387 !    RECUPERATE NODES OF THE EDGE WITH THE GOOD ORIENTATION
00388 !     WITH RESPECT TO THE NORMAL
00389             NUBO1 = NUBO(1,NSG)
00390             NUBO2 = NUBO(2,NSG)
00391             PROD_SCAL= ((X(NUBO2)-X(NUBO1))*VNOCL(1,NSG)+
00392      &                  (Y(NUBO2)-Y(NUBO1))*VNOCL(2,NSG))
00393             IF(PROD_SCAL.LT.0.D0)THEN
00394               NUBO1 = NUBO(2,NSG)
00395               NUBO2 = NUBO(1,NSG)
00396             ENDIF
00397 !
00398             ZF1        = ZF(NUBO1)
00399             ZF2        = ZF(NUBO2)
00400             IF(PROD_SCAL.LT.0.D0)THEN
00401               DSZ(1,NSG) = DSZ0(2,NSG)
00402               DSZ(2,NSG) = DSZ0(1,NSG)
00403             ELSE
00404               DSZ(1,NSG) = DSZ0(1,NSG)
00405               DSZ(2,NSG) = DSZ0(2,NSG)
00406             ENDIF
00407 !
00408             HI0=UA(1,NUBO1)
00409             HJ0=UA(1,NUBO2)
00410 !
00411 !   FOR AN EDGE BEING RECOVERED, ONE WILL REMAIN 1ST ORDER
00412 !
00413             IF(ZF1.GE. (HJ0+ZF2) .OR. ZF2.GE. (HI0+ZF1)
00414      &     .OR. 2.*ABS(DSZ(1,NSG)).GE.HI0
00415      &     .OR. 2.*ABS(DSZ(1,NSG)).GE.HJ0
00416      &     .OR. 2.*ABS(DSZ(2,NSG)).GE.HI0
00417      &     .OR. 2.*ABS(DSZ(2,NSG)).GE.HJ0)  THEN
00418 !ra02/05/2013 FOR APTIMIZATION
00419               CYCLE
00420             ELSE
00421 !
00422 !   ONE REBUILDS H+Z, DSH = VARIATION OF H+Z
00423 !
00424               ILIM=1
00425               BETA=1.D0
00426 !
00427               DSH(1,NSG) = EXLIM(ILIM,BETA,GRADI(1,NSG),GRADIJ(1,NSG))
00428               DSH(2,NSG) = EXLIM(ILIM,BETA,GRADJ(1,NSG),GRADJI(1,NSG))
00429               !   FOR PARALLELILSM
00430               IF(NCSIZE.GT.1.AND.IFABOR(IEL,I).EQ.-2)THEN ! THIS IS AN INTERFACE EDGE
00431                 IF(DSH(1,NSG).GE.0.D0) THEN
00432                   DSP(NUBO1) = DSP(NUBO1)+DEMI*AIRST(1,NSG)*DSH(1,NSG) ! WE CONSIDER ONLY
00433                 ELSE                                                  ! 0.5 AIRST
00434                   DSM(NUBO1) = DSM(NUBO1)-DEMI*AIRST(1,NSG)*DSH(1,NSG) ! PARCOM2 WILL ADD
00435                 ENDIF                                                 ! CONTRIBUTIONS
00436                 IF(DSH(2,NSG).GE.0.D0) THEN
00437                   DSP(NUBO2) = DSP(NUBO2)+DEMI*AIRST(2,NSG)*DSH(2,NSG)
00438                 ELSE
00439                   DSM(NUBO2) = DSM(NUBO2)-DEMI*AIRST(2,NSG)*DSH(2,NSG)
00440                 ENDIF
00441               ELSE ! NO PARALLELILSM OR NO INTERFACE EDGE
00442                 IF(DSH(1,NSG).GE.0.D0) THEN
00443                   DSP(NUBO1) = DSP(NUBO1) + AIRST(1,NSG)*DSH(1,NSG)
00444                 ELSE
00445                   DSM(NUBO1) = DSM(NUBO1) - AIRST(1,NSG)*DSH(1,NSG)
00446                 ENDIF
00447                 IF(DSH(2,NSG).GE.0.D0) THEN
00448                   DSP(NUBO2) = DSP(NUBO2) + AIRST(2,NSG)*DSH(2,NSG)
00449                 ELSE
00450                   DSM(NUBO2) = DSM(NUBO2) - AIRST(2,NSG)*DSH(2,NSG)
00451                 ENDIF
00452               ENDIF
00453 !
00454               ILIM=2
00455               BETA=0.3333D0 ! THESE ARE CHOICES OF INRIA 1/3 FOR
00456                             ! VELOCITIES AND 1 FOR H
00457 !
00458               DSU(1,NSG) = EXLIM(ILIM,BETA,GRADI(2,NSG),GRADIJ(2,NSG))
00459               DSU(2,NSG) = EXLIM(ILIM,BETA,GRADJ(2,NSG),GRADJI(2,NSG))
00460 !
00461               DSV(1,NSG) = EXLIM(ILIM,BETA,GRADI(3,NSG),GRADIJ(3,NSG))
00462               DSV(2,NSG) = EXLIM(ILIM,BETA,GRADJ(3,NSG),GRADJI(3,NSG))
00463 !
00464             ENDIF
00465             YESNO(NSG)=.TRUE.
00466           ENDIF
00467         ENDDO
00468       ENDDO
00469       !  FOR PARALLELILSM
00470       IF(NCSIZE.GT.1)THEN
00471         CALL PARCOM2(DSP,DSM,DSM,NS,1,2,2,MESH)
00472       ENDIF
00473 !
00474 !  ONE CALCULATES THE CORRECTIONS TO ENSURE THE CONSERVATION OF H
00475 !
00476       DO IS=1,NS
00477         CORR(IS) =  DSM(IS) - DSP(IS)
00478         AMDS =MAX(DSP(IS),DSM(IS))
00479         IF(AMDS.GT.0.D0) THEN
00480           CORR(IS) = CORR(IS)/AMDS
00481         ENDIF
00482       ENDDO
00483 ! IF ORDER 2 REINITIALIZATION OF YESNO
00484       DO I=1,NSEG
00485         YESNO(I)=.FALSE.
00486       ENDDO
00487 ! FIRST ORDER
00488  12   CONTINUE
00489 !
00490 !     LOOP ON GLOBAL LIST OF EDGES
00491 !    ******************************
00492 !
00493       DO IEL=1, NELEM
00494         DO I = 1,3
00495           IF(.NOT.YESNO(ELTSEG(IEL,I)))THEN
00496             NSG = ELTSEG(IEL,I)
00497 !
00498 !    RECUPERATE NODES OF THE EDGE WITH THE GOOD ORIENTATION
00499 !     WITH RESPECT TO THE NORMAL
00500             NUBO1 = NUBO(1,NSG)
00501             NUBO2 = NUBO(2,NSG)
00502             PROD_SCAL= ((X(NUBO2)-X(NUBO1))*VNOCL(1,NSG)+
00503      &                  (Y(NUBO2)-Y(NUBO1))*VNOCL(2,NSG))
00504             IF(PROD_SCAL.LT.0.D0)THEN
00505              NUBO1 = NUBO(2,NSG)
00506              NUBO2 = NUBO(1,NSG)
00507             ENDIF
00508 !
00509             ZF1  = ZF(NUBO1)
00510             ZF2  = ZF(NUBO2)
00511             DSZ1 = 0.D0
00512             DSZ2 = 0.D0
00513 !
00514             XNN  = VNOCL(1,NSG)
00515             YNN  = VNOCL(2,NSG)
00516             RNN  = VNOCL(3,NSG)
00517 !
00518             UAS11 = UA(1,NUBO1)
00519             UAS12 = UA(1,NUBO2)
00520             UAS21 = UA(2,NUBO1)
00521             UAS22 = UA(2,NUBO2)
00522             UAS31 = UA(3,NUBO1)
00523             UAS32 = UA(3,NUBO2)
00524 !
00525             HI0 = UAS11
00526             HJ0 = UAS12
00527 !
00528 ! ROTATION
00529 !
00530             UAS210 = UAS21
00531             UAS21  = XNN*UAS210+YNN*UAS31
00532             UAS31  =-YNN*UAS210+XNN*UAS31
00533 !
00534             UAS220 = UAS22
00535             UAS22  = XNN*UAS220+YNN*UAS32
00536             UAS32  =-YNN*UAS220+XNN*UAS32
00537 !
00538             IF(NORDRE.EQ.1)  GOTO 1234
00539 !
00540 !    REBUILDING FOR 2ND ORDER
00541 !    ***************************
00542 !
00543 !   FOR AN EDGE BEING RECOVERED, ONE REMAINS 1ST ORDER
00544 !
00545             IF(ZF1.GE. (HJ0+ZF2) .OR. ZF2.GE. (HI0+ZF1)
00546      &     .OR. 2.*ABS(DSZ(1,NSG)).GE.HI0
00547      &     .OR. 2.*ABS(DSZ(1,NSG)).GE.HJ0
00548      &     .OR. 2.*ABS(DSZ(2,NSG)).GE.HI0
00549      &     .OR. 2.*ABS(DSZ(2,NSG)).GE.HJ0)  GOTO 1234
00550 !
00551 !
00552             UAS11  = UAS11  + DSH(1,NSG)  -  DSZ(1,NSG) +
00553      &        MIN(0.D0,CORR(NUBO1))*MAX(0.D0,DSH(1,NSG))+
00554      &        MAX(0.D0,CORR(NUBO1))*MAX(0.D0,-DSH(1,NSG))
00555 !
00556             UAS21  = UAS21 + DSU(1,NSG)
00557             UAS31  = UAS31 + DSV(1,NSG)
00558             DSZ1   = DSZ(1,NSG)
00559             ZF1    = ZF1 + DSZ1
00560 !
00561             UAS12  = UAS12  + DSH(2,NSG)  -  DSZ(2,NSG) +
00562      &        MIN(0.D0,CORR(NUBO2))*MAX(0.D0,DSH(2,NSG))+
00563      &        MAX(0.D0,CORR(NUBO2))*MAX(0.D0,-DSH(2,NSG))
00564 !
00565             UAS22  = UAS22 + DSU(2,NSG)
00566             UAS32  = UAS32 + DSV(2,NSG)
00567             DSZ2   = DSZ(2,NSG)
00568             ZF2    = ZF2 + DSZ2
00569 !
00570 !
00571             IF(UAS11.LE.0.D0) THEN
00572               UAS11 = 0.D0
00573               UAS21 = 0.D0
00574               UAS31 = 0.D0
00575             ENDIF
00576             IF(UAS12.LE.0.D0)  THEN
00577               UAS12 = 0.D0
00578               UAS22 = 0.D0
00579               UAS32 = 0.D0
00580             ENDIF
00581 !
00582 !
00583 !    LIMITATION OF THE TIME STEP
00584 !    ***************************
00585 !
00586 !
00587             SIGMAX=MAX( 1.D-2, RA3* SQRT(UAS11) + ABS(UAS21) )
00588             DTL    = CFL*AIRST(1,NSG)/(RNN*SIGMAX)
00589             DTLL(NUBO1) = MIN (DTL,DTLL(NUBO1))
00590             DT          = MIN(DT, DTL)
00591 !
00592             SIGMAX=MAX( 1.D-2, RA3* SQRT(UAS12) + ABS(UAS22) )
00593             DTL    = CFL*AIRST(2,NSG)/(RNN*SIGMAX)
00594             DTLL(NUBO2) = MIN (DTL,DTLL(NUBO2))
00595             DT          = MIN(DT, DTL)
00596 ! PARALLEL: TAKE MIN DT OF ALL SUBDOMAINS
00597 !          IF(NCSIZE.GT.1)DT = P_DMIN(DT) !WILL BE PLACED AT THE END (SEE BELOW)
00598 !
00599 1234        CONTINUE
00600 !
00601 !
00602 !   MAIN FLUX COMPUTATAION
00603 !
00604             HI        = UAS11
00605             HC(1,NSG) = UAS11
00606 !
00607             DZIJ = MAX(0.D0,ZF2-ZF1)
00608             HIJ  = MAX(0.D0, HI- DZIJ)
00609 !
00610 !
00611             IF(HIJ.LE.0.D0) THEN
00612               CIJ  = 0.D0
00613               FLU11= 0.D0
00614               FLU21= 0.D0
00615             ELSE
00616               CIJ  = SQRT(HIJ)
00617               EXT1 = MIN(RA3,MAX(-RA3,-UAS21/CIJ))
00618 !
00619               A01  = ALP*(RA3-EXT1)
00620               A11  = ALP*(RA3**2-EXT1**2)/2.D0
00621               A21  = ALP*(RA3**3-EXT1**3)/3.D0
00622 !
00623               FLU11= HIJ*(UAS21*A01+CIJ*A11)
00624               FLU21= UAS21*(FLU11+CIJ*HIJ*A11) +A21*HIJ*HIJ
00625             ENDIF
00626 !
00627 !
00628             HJ   = UAS12
00629             HC(2,NSG) = UAS12
00630 !
00631             DZJI = MAX(0.D0,ZF1-ZF2)
00632             HJI  = MAX(0.D0, HJ- DZJI)
00633 !
00634             IF(HJI.LE.0.D0) THEN
00635               CJI   = 0.D0
00636               FLU12 = 0.D0
00637               FLU22 = 0.D0
00638             ELSE
00639               CJI  = SQRT(HJI)
00640               EXT2 = MIN(RA3,MAX(-RA3,-UAS22/CJI))
00641 !
00642               A02  = ALP*(RA3+EXT2)
00643               A12  = ALP*(EXT2**2-RA3**2)/2.D0
00644               A22  = ALP*(RA3**3+EXT2**3)/3.D0
00645 !
00646               FLU12= HJI*(UAS22*A02+CJI*A12)
00647               FLU22= UAS22*(FLU12+CJI*HJI*A12) +A22*HJI*HJI
00648             ENDIF
00649 !
00650             HGZI = 0.5D0*RNN*(HIJ+HI)*(HIJ-HI)
00651             HGZJ = 0.5D0*RNN*(HJI+HJ)*(HJI-HJ)
00652 !
00653             IF(NORDRE.EQ.2) THEN
00654               HGZI = HGZI - 0.5D0*RNN*(HI0+HI)*DSZ1
00655               HGZJ = HGZJ - 0.5D0*RNN*(HJ0+HJ)*DSZ2
00656             ENDIF
00657 !
00658             FLU11=(FLU11+FLU12)*RNN
00659             FLU21=(FLU21+FLU22)*RNN
00660 !
00661             IF(NTRAC.GT.0) THEN
00662               DO ITRAC=1,NTRAC
00663                 FLUXTEMP%ADR(ITRAC)%P%R(NSG)=FLU11
00664               ENDDO
00665             ENDIF
00666 !
00667             IF(FLU11.GE.0.D0) THEN
00668               FLU31 =  UAS31 * FLU11
00669             ELSE
00670               FLU31 =  UAS32 * FLU11
00671             ENDIF
00672 !
00673 ! OPPOSITE ROTATION
00674 !
00675             FLU210 = FLU21
00676             FLU21  = XNN*FLU210-YNN*FLU31
00677             FLU31  = YNN*FLU210+XNN*FLU31
00678 !
00679 !       TERM DUE TO THE BOTTOM GRADIENT
00680 !
00681             HDXZ1  = G*XNN*HGZI
00682             HDYZ1  = G*YNN*HGZI
00683 !
00684             HDXZ2  = G*XNN*HGZJ
00685             HDYZ2  = G*YNN*HGZJ
00686 !
00687 !FOR PARALLELISM
00688 !
00689             IF(NCSIZE.GT.1)THEN
00690               IF(IFABOR(IEL,I).EQ.-2)THEN !THIS IS INTERFACE EDGE
00691                 FLU11= DEMI*FLU11
00692                 FLU21= DEMI*FLU21
00693                 FLU31= DEMI*FLU31
00694 !
00695                 HDXZ1 = DEMI*HDXZ1
00696                 HDYZ1 = DEMI*HDYZ1
00697                 HDXZ2 = DEMI*HDXZ2
00698                 HDYZ2 = DEMI*HDYZ2
00699                 IF(NTRAC.GT.0) THEN
00700                   DO ITRAC=1,NTRAC
00701                     FLUXTEMP%ADR(ITRAC)%P%R(NSG)=DEMI*FLU11
00702                   ENDDO
00703                 ENDIF
00704               ENDIF
00705             ENDIF
00706 !
00707 !***********************************************************
00708 !       FLUX INCREMENT
00709 !***********************************************************
00710 !
00711             CE(NUBO1,1) = CE(NUBO1,1) - FLU11
00712             CE(NUBO1,2) = CE(NUBO1,2) - FLU21 + HDXZ1
00713             CE(NUBO1,3) = CE(NUBO1,3) - FLU31 + HDYZ1
00714 !
00715             CE(NUBO2,1) = CE(NUBO2,1) + FLU11
00716             CE(NUBO2,2) = CE(NUBO2,2) + FLU21 - HDXZ2
00717             CE(NUBO2,3) = CE(NUBO2,3) + FLU31 - HDYZ2
00718 !
00719             YESNO(NSG)=.TRUE.
00720           ENDIF
00721         ENDDO
00722       ENDDO
00723 !
00724 !       LIMITATION OF THE TIME STEP FOR THE BOUNDARY NODES
00725 !
00726       IF(NORDRE.EQ.2) THEN
00727         IF(NPTFR.GT.0)THEN  ! USEFUL FOR PARALLEL CASE
00728           DO K=1,NPTFR
00729             IS = NBOR(K)
00730             VNX= XNEBOR(K+NPTFR)
00731             VNY= YNEBOR(K+NPTFR)
00732             VNL= SQRT(VNX**2+VNY**2)
00733             SIGMAX= SQRT(UA(1,IS))
00734             UNORM=SQRT(UA(2,IS)*UA(2,IS) + UA(3,IS)*UA(3,IS))
00735             SIGMAX=MAX( 1.D-2, RA3*SIGMAX +UNORM )
00736             DTL   = CFL*AIRS(IS)/(VNL*SIGMAX)
00737             AUX   = DTL/DTLL(IS)
00738             AUX1  =AUX/(1.D0+AUX)
00739             DT    =MIN(DT, AUX1*DTLL(IS))
00740           ENDDO
00741         ENDIF
00742 !       FOR PARALLELISME
00743         IF(NCSIZE.GT.1) DT=P_DMIN(DT)
00744       ENDIF
00745 !
00746       DEALLOCATE(YESNO)
00747 !
00748 !-----------------------------------------------------------------------
00749 !
00750       RETURN
00751       END

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