flucint.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\flucint.f
00002 !
00051                      SUBROUTINE FLUCINT
00052 !                    ******************
00053 !
00054      &(NS,NSEG,DIMT,NUBO,G,X,Y,UA,TN,ZF,VNOCL,CE,
00055      & NORDRE,CMI,JMI,DJX,DJY,DX,DY,DJXT,DJYT,DXT,DYT,EPSWL)
00056 !
00057 !***********************************************************************
00058 ! TELEMAC2D   V6P2                                   21/08/2010
00059 !***********************************************************************
00060 !
00061 !
00062 !
00063 !
00064 !
00065 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00066 !| CE             |<->| FLUX INCREMENTS AT INTERNAL INTERFACES
00067 !| CMI            |-->| COORDINATES OF MIDDLE PONTS OF EDGES
00068 !| DIMT           |-->| DIMENSION OF TN
00069 !| DJX            |-->| GRADIENTS PER TRIANGLE
00070 !| DJXT           |-->| WORKING TABLES FOR TRACER
00071 !| DJY            |-->| GRADIENTS PER TRIANGLE
00072 !| DJYT           |---| WORKING TABLES FOR TRACER
00073 !| DX             |-->| GRADIENTS PER NODE
00074 !| DXT            |-->| WORKING TABLES FOR TRACER
00075 !| DY             |-->| GRADIENTS PER NODE
00076 !| DYT            |-->| WORKING TABLES FOR TRACER
00077 !| EPSWL          |-->| THRESHOLD DEPTH
00078 !| G              |-->| GRAVITY
00079 !| JMI            |-->| NUMBER OF THE TRIANGLE IN WHICH IS LOCATED
00080 !|                |   | THE MIDDLE POINT OF THE INTERFACE
00081 !| NORDRE         |-->| ORDER OF THE SCHEME
00082 !| NS             |-->| TOTAL NUMBER OF NODES IN THE MESH
00083 !| NSEG           |-->| TOTAL NUMBER OF SGMENTS IN THE MESH
00084 !| NUBO           |-->| GLOBAL NUMBERS OF THE NODES FORMING THE EDGE
00085 !| TN             |-->| CURRENT TIME
00086 !| UA             |-->| UA(1,IS) = H,  UA(2,IS)=U  ,UA(3,IS)=V
00087 !| VNOCL          |-->| NORMAL VECTOR TO THE INTERFACE
00088 !|                |   | (2 FIRST COMPONENTS) AND
00089 !|                |   | LENGTH OF THE SEGMENT (3RD COMPONENT)
00090 !| X              |-->| X COORDINATES OF THE NODES
00091 !| Y              |-->| Y COORDINATES OF THE NODES
00092 !| ZF             |-->| BATHYMETRY
00093 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00094 !
00095       IMPLICIT NONE
00096       INTEGER LNG,LU
00097       COMMON/INFO/LNG,LU
00098 !
00099 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00100 !
00101       INTEGER, INTENT(IN) :: NS,NSEG,DIMT,NORDRE
00102       INTEGER, INTENT(IN) :: NUBO(2,NSEG)
00103       DOUBLE PRECISION, INTENT(IN) :: G,X(NS),Y(NS),VNOCL(3,NSEG),ZF(NS)
00104       DOUBLE PRECISION, INTENT(INOUT) :: CE(NS,3)
00105       DOUBLE PRECISION, INTENT(IN)    :: UA(3,NS),TN(DIMT),CMI(2,*)
00106       DOUBLE PRECISION, INTENT(IN)    :: JMI(*),EPSWL
00107       DOUBLE PRECISION, INTENT(IN)    :: DJX(3,*),DJY(3,*)
00108       DOUBLE PRECISION, INTENT(IN)    :: DX(3,*),DY(3,*)
00109       DOUBLE PRECISION, INTENT(IN)    :: DJXT(*),DJYT(*),DXT(*),DYT(*)
00110 !
00111 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00112 !
00113       INTEGER NSG,NUBO1,NUBO2,J,IVAR
00114 !
00115       DOUBLE PRECISION E2,RA2,RA3,ALP,ZF1,ZF2,XNN,YNN,RNN
00116       DOUBLE PRECISION UAS11,UAS12,UAS21,UAS22,UAS31,UAS32
00117       DOUBLE PRECISION GRADI(4),GRADJ(4),GRADIJ(4),GRADJI(4),GRIJ,GRJI
00118       DOUBLE PRECISION GRI,GRJ,AUX1,AUX2,GRI2,GRIJ2,GRJ2,GRJI2
00119       DOUBLE PRECISION AIX,AIY,AJX,AJY,HI,HI0,HIJ,CIJ,ETAI,UAS110,UAS120
00120       DOUBLE PRECISION HJ,HJ0,HJI,CJI,ETAJ,EXT1,EXT2,FLU11,FLU12,FLU41
00121       DOUBLE PRECISION A01,A11,A02,A12,UAS41,UAS42
00122       DOUBLE PRECISION UAS410,UAS420,BETA,BETA1
00123 !
00124 !-----------------------------------------------------------------------
00125 !
00126       RA2  = SQRT(2.D0)
00127       RA3  = SQRT(1.5D0*G)
00128       ALP  = 0.5D0/RA3
00129 !
00130       BETA = 0.333D0
00131       BETA1 = 1.D0+BETA
00132       E2    = 1.E-6
00133 !
00134 !     LOOP ON GLOBAL LIST OF EDGES
00135 !
00136       DO NSG=1,NSEG
00137         J         = INT(JMI(NSG))
00138 !
00139         NUBO1     = NUBO(1,NSG)
00140         NUBO2     = NUBO(2,NSG)
00141 !
00142         ZF1   =    ZF(NUBO1)
00143         ZF2   =    ZF(NUBO2)
00144 !
00145         XNN       = VNOCL(1,NSG)
00146         YNN       = VNOCL(2,NSG)
00147         RNN       = VNOCL(3,NSG)
00148 !
00149         UAS11     = UA(1,NUBO1)
00150         UAS12     = UA(1,NUBO2)
00151         UAS21     = UA(2,NUBO1)
00152         UAS22     = UA(2,NUBO2)
00153         UAS31     = UA(3,NUBO1)
00154         UAS32     = UA(3,NUBO2)
00155         UAS41     = TN(NUBO1)
00156         UAS42     = TN(NUBO2)
00157 !
00158         HI0=UAS11*UAS11
00159         HJ0=UAS12*UAS12
00160 !
00161 ! ROTATION
00162 !
00163         UAS21  = XNN*UAS21+YNN*UAS31
00164 !
00165         UAS22  = XNN*UAS22+YNN*UAS32
00166 !
00167 !
00168         IF(NORDRE.EQ.1)  GOTO 1234
00169 !
00170         AIX       = CMI(1,NSG)-X(NUBO1)
00171         AIY       = CMI(2,NSG)-Y(NUBO1)
00172         AJX       = CMI(1,NSG)-X(NUBO2)
00173         AJY       = CMI(2,NSG)-Y(NUBO2)
00174 !
00175         DO IVAR=1,3
00176 !
00177           GRADI(IVAR)  = AIX*DX(IVAR,NUBO1) + AIY*DY(IVAR,NUBO1)
00178 !
00179           GRADJ(IVAR)  = AJX*DX(IVAR,NUBO2) + AJY*DY(IVAR,NUBO2)
00180 !
00181           GRADIJ(IVAR)  = AIX*DJX(IVAR,J) + AIY*DJY(IVAR,J)
00182 !
00183           GRADJI(IVAR)  = AJX*DJX(IVAR,J) + AJY*DJY(IVAR,J)
00184 !
00185         ENDDO
00186 !
00187         GRADI(4)  = AIX*DXT(NUBO1) + AIY*DYT(NUBO1)
00188 !
00189         GRADJ(4)  = AJX*DXT(NUBO2) + AJY*DYT(NUBO2)
00190 !
00191         GRADIJ(4)  = AIX*DJXT(J) + AIY*DJYT(J)
00192 !
00193         GRADJI(4)  = AJX*DJXT(J) + AJY*DJYT(J)
00194 !
00195 ! GRADIENTS ROTATION
00196 !
00197 !
00198         GRADI(2)  = XNN*GRADI(2)+YNN*GRADI(3)
00199 !
00200         GRADIJ(2)  = XNN*GRADIJ(2)+YNN*GRADIJ(3)
00201 !
00202         GRADJ(2)  = XNN*GRADJ(2)+YNN*GRADJ(3)
00203 !
00204         GRADJI(2)  = XNN*GRADJI(2)+YNN*GRADJI(3)
00205 !
00206         DO IVAR=1,4
00207           IF(IVAR.NE.3) THEN
00208 !
00209             GRIJ = GRADIJ(IVAR)
00210             GRJI = GRADJI(IVAR)
00211 !
00212             GRI = BETA1*GRADI(IVAR) - BETA*GRIJ
00213             GRJ = BETA1*GRADJ(IVAR) - BETA*GRJI
00214 !
00215             AUX1 = 0.5*(1.0+  DSIGN(1.0D0, GRI*GRIJ))
00216             AUX2 = 0.5*(1.0 + DSIGN(1.0D0, GRJ*GRJI))
00217 !
00218 !    VAN ALBADA
00219 !
00220             GRI2  = GRI*GRI    + E2
00221             GRIJ2 = GRIJ*GRIJ  + E2
00222             GRJ2  = GRJ*GRJ    + E2
00223             GRJI2 = GRJI*GRJI  + E2
00224 !
00225             GRADI(IVAR)  = AUX1*
00226      &         (GRI2  *GRIJ  + GRIJ2  *GRI )/(GRI2 + GRIJ2)
00227 !
00228             GRADJ(IVAR)  = AUX2*
00229      &         (GRJ2  *GRJI + GRJI2 *GRJ )/(GRJ2 + GRJI2)
00230 !
00231           ENDIF
00232         ENDDO
00233 !
00234 !
00235         UAS110 =   UAS11
00236         UAS11 = MIN(MAX(UAS110/RA2,UAS11 + GRADI(1)),RA2*UAS110)
00237         UAS21     = UAS21 + GRADI(2)
00238 !
00239         UAS120 =   UAS12
00240         UAS12 = MIN(MAX(UAS120/RA2,UAS12 + GRADJ(1)),RA2*UAS120)
00241         UAS22     = UAS22 + GRADJ(2)
00242 !
00243         IF(UAS11.LE.0.D0) THEN
00244           UAS11 =0.D0
00245           UAS21 =0.D0
00246         ENDIF
00247         IF(UAS12.LE.0.D0)  THEN
00248            UAS12 =0.D0
00249            UAS22 =0.D0
00250         ENDIF
00251 !
00252 !
00253 1234    CONTINUE
00254 !
00255         HI   = UAS11*UAS11
00256 !
00257 !   ETAI = FREE SURFACE ELEVATION
00258 !   IF HI < EPSWL KEEP FLAT BED
00259 !
00260         ETAI = HI0+ZF1
00261         IF(HI0.LE.EPSWL) ETAI = MIN(ETAI,HJ0+ZF2)
00262 !
00263 !       HIJ= ETAI-MIN(ETAI,MAX(ZF1,ZF2))
00264 !
00265         IF(ZF1.GE.ZF2) THEN
00266           HIJ=HI0
00267         ELSE
00268           HIJ=MAX(0.D0,ETAI-ZF2)
00269         ENDIF
00270 !
00271         IF(HI.LE.0.D0) THEN
00272           CIJ=0.D0
00273           FLU11=0.D0
00274         ELSE
00275 !
00276           IF(HIJ.LE.0.D0) THEN
00277             CIJ=0.D0
00278 !
00279             IF(UAS21.GE.0.D0) THEN
00280               EXT1=-RA3
00281             ELSE
00282               EXT1=RA3
00283             ENDIF
00284 !
00285           ELSE
00286             CIJ  =HIJ*SQRT(HIJ)/HI
00287             EXT1 = MIN(RA3,MAX(-RA3,-UAS21/CIJ))
00288           ENDIF
00289 !
00290         A01  = ALP*(RA3-EXT1)
00291         A11  = ALP*(RA3**2-EXT1**2)/2.D0
00292 !
00293         FLU11= HI*(UAS21*A01+CIJ*A11)
00294 !
00295         ENDIF
00296 !
00297 !
00298         HJ   = UAS12 *UAS12
00299 !
00300         ETAJ = HJ0+ZF2
00301         IF(HJ0.LE.EPSWL) ETAJ = MIN(ETAJ,HI0+ZF1)
00302 !
00303 !       HJI= ETAJ-MIN(ETAJ,MAX(ZF1,ZF2))
00304         IF(ZF2.GE.ZF1) THEN
00305           HJI=HJ0
00306         ELSE
00307           HJI=MAX(0.D0,ETAJ-ZF1)
00308         ENDIF
00309 !
00310         IF(HJ.LE.0.D0) THEN
00311           CJI=0.D0
00312           FLU12=0.D0
00313         ELSE
00314 !
00315           IF(HJI.LE.0.D0) THEN
00316             CJI=0.D0
00317 !
00318             IF(UAS22.GE.0.D0) THEN
00319               EXT2=-RA3
00320             ELSE
00321               EXT2=RA3
00322             ENDIF
00323 !
00324           ELSE
00325             CJI  =HJI*SQRT(HJI)/HJ
00326             EXT2 = MIN(RA3,MAX(-RA3,-UAS22/CJI))
00327           ENDIF
00328 !
00329           A02  = ALP*(RA3+EXT2)
00330           A12  = ALP*(EXT2**2-RA3**2)/2.D0
00331 !
00332           FLU12= HJ*(UAS22*A02+CJI*A12)
00333         ENDIF
00334 !
00335 !
00336         FLU11=(FLU11 + FLU12)*RNN
00337 !
00338         IF(FLU11.GE.0.D0) THEN
00339           IF(NORDRE.GE.2) THEN
00340             UAS410 = UAS41
00341             UAS41 = MIN(MAX(0.5D0*UAS410,UAS41 + GRADI(4)),2.D0*UAS410)
00342           ENDIF
00343           FLU41 =  UAS41 * FLU11
00344         ELSE
00345           IF(NORDRE.GE.2) THEN
00346             UAS420 = UAS42
00347             UAS42 = MIN(MAX(0.5D0*UAS420,UAS42 + GRADJ(4)),2.D0*UAS420)
00348           ENDIF
00349           FLU41 =  UAS42 * FLU11
00350         ENDIF
00351 !
00352         CE(NUBO1,1) = CE(NUBO1,1) - FLU41
00353         CE(NUBO2,1) = CE(NUBO2,1) + FLU41
00354 !
00355       ENDDO ! NSG
00356 !
00357 !-----------------------------------------------------------------------
00358 !
00359       RETURN
00360       END

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