cdl.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\cdl.f
00002 !
00056                         SUBROUTINE CDL
00057 !                       **************
00058 !
00059      &(NS,NPTFR,NBOR,LIMPRO,XNEBOR,YNEBOR,KDIR,KNEU,G,HBOR,
00060      & UBOR,VBOR,UA,CE,FLUENT,FLUSORT,FLBOR,
00061      & DTHAUT,DT,CFL,FLUHBTEMP,NTRAC)
00062 !
00063 !***********************************************************************
00064 ! TELEMAC 2D VERSION 6.1                                          INRIA
00065 !***********************************************************************
00066 !
00067 !
00068 !    UA(1,IS) = H,  UA(2,IS)=U  ,UA(3,IS)=V
00069 !
00070 !
00071 !
00072 !
00073 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00074 !|  NS            |-->|  TOTAL NUMNER OF NODES
00075 !|  NPTFR         |-->|  TOTAL NUMBER OF BOUNDARY NODES
00076 !|  NBOR          |-->|  GLOBAL NUMBERS OF BOUNDARY POINTS
00077 !|  LIMPRO        |-->|  TYPES OF BOUNDARY CONDITION
00078 !|  XNEBOR        |-->|  UNIT OUTWARD NORMAL COMPONENT AT BOUNDARY POINTS
00079 !|  YNEBOR        |-->|  UNIT OUTWARD NORMAL COMPONENT AT BOUNDARY POINTS
00080 !|  KDIR          |-->|  CONVENTION FOR DIRICHLET POINTS
00081 !|  KNEU          |-->|  CONVENTION FOR NEUMANN POINTS
00082 !|  G             |-->|  GRAVITY CONSTANT
00083 !|  HBOR          |-->|  IMPOSED VALUES FOR H
00084 !|  UBOR          |-->|  IMPOSED VALUES FOR U
00085 !|  VBOR          |-->|  IMPOSED VALUES FOR V
00086 !|  UA            |-->|  UA(1,IS) = H,  UA(2,IS)=U  ,UA(3,IS)=V
00087 !|  CE            |<->|  FLUX
00088 !|  FLUENT        |<--|  ENTERING MASS FLUX
00089 !|  FLUSORT       |<--|  EXITING MASS FLUX
00090 !|  DTHAUT        |-->|  CHARACTERISTIC LENGTH (DX) FOR CFL
00091 !|  DT            |<->|  TIME STEP
00092 !|  CFL           |-->|  CFL NUMBER
00093 !|  FLUHBTEMP     |<--|  BOUNDARY FLUX FOR THE TRACER
00094 !|  TRAC          |-->|  LOGICAL: TO INDICATE THE PRESENCE OF A TRACER
00095 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00096 !
00097       USE BIEF
00098       IMPLICIT NONE
00099       INTEGER LNG,LU
00100       COMMON/INFO/LNG,LU
00101 !
00102 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00103 !
00104       INTEGER, INTENT(IN)             :: NS,NPTFR,KDIR,KNEU,NTRAC
00105       INTEGER, INTENT(IN)             :: NBOR(NPTFR),LIMPRO(NPTFR,6)
00106       DOUBLE PRECISION, INTENT(IN)    :: XNEBOR(2*NPTFR),YNEBOR(2*NPTFR)
00107       DOUBLE PRECISION, INTENT(IN)    :: HBOR(NPTFR),UA(3,NS),DTHAUT(*)
00108       DOUBLE PRECISION, INTENT(IN)    :: UBOR(NPTFR),VBOR(NPTFR)
00109       DOUBLE PRECISION, INTENT(IN)    :: G,CFL
00110       DOUBLE PRECISION, INTENT(INOUT) :: DT
00111       DOUBLE PRECISION, INTENT(INOUT) :: CE(NS,3),FLUENT,FLUSORT
00112       TYPE(BIEF_OBJ), INTENT(INOUT)   :: FLUHBTEMP
00113       TYPE(BIEF_OBJ), INTENT(INOUT)   :: FLBOR
00114 !
00115 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00116 !
00117       INTEGER IS,K,NIT,ITRAC
00118 !
00119       DOUBLE PRECISION RA3,RA32,RA33, ALP,ALP2,ALP3,SG,SQ2
00120       DOUBLE PRECISION VNX,VNY,VNX1,VNY1,VNL,H,U,V,RUN
00121       DOUBLE PRECISION FLUH,FLUU,FLUV,AUX,FLUTMP,RH,HRH,UNN,VNN
00122       DOUBLE PRECISION FHPLUS,FUPLUS,FHMOINS,FUMOINS
00123       DOUBLE PRECISION A,A1,A2,A3,ALPHA0,ALPHA1,ALPHA2,C,VP1,VP2 ,VP3
00124       DOUBLE PRECISION HG ,RHG,HRHG,UG,VG,DEST,RVG,CA1,AM
00125       DOUBLE PRECISION UIN,VIN,HUIN,HVIN,SIGMAX,DTL,UNORM
00126       DOUBLE PRECISION P_DMIN
00127       EXTERNAL P_DMIN
00128 !
00129       SQ2   = SQRT(2.D0)
00130       SG    = SQRT(G)
00131       RA3   = SQRT(1.5D0*G)
00132       RA32  = RA3**2
00133       RA33  = RA3*RA32
00134       ALP   = 0.5D0/RA3
00135       ALP2  = 0.5D0 *ALP
00136       ALP3  = ALP/3.D0
00137 !
00138       FLUENT=0.D0
00139       FLUSORT=0.D0
00140 !
00141       IF(NPTFR.GT.0)THEN ! USEFUL FOR PARALLEL
00142         DO K=1,NPTFR
00143           IS=NBOR(K)
00144           VNX1=XNEBOR(K)
00145           VNY1=YNEBOR(K)
00146           VNX=XNEBOR(K+NPTFR)
00147           VNY=YNEBOR(K+NPTFR)
00148           VNL=SQRT(VNX**2+VNY**2)
00149 !
00150           H   = UA(1,IS)
00151           RH  = SQRT(H)
00152           U   = UA(2,IS)
00153           V   = UA(3,IS)
00154 !
00155 !         SOLID WALLS
00156 !         **************
00157 
00158 !     PERFECT SLIPPING CONDITION
00159 !     **************************
00160 !
00161           IF(LIMPRO(K,1).EQ.KNEU) THEN
00162 !
00163             AUX=0.5D0*G*H**2
00164             FLUH = 0.D0
00165             FLUU = AUX*VNX
00166             FLUV = AUX*VNY
00167           ELSE
00168 !
00169 !        LIQUID BOUNDARY
00170 !        *******************
00171 !
00172 !    CALCULATION OF F+(H,U,V)
00173 !
00174             HRH = RH * H
00175 !
00176             IF(H.LE.0.D0) THEN
00177               U=0.D0
00178               V=0.D0
00179               UNN=0.D0
00180               VNN=0.D0
00181               FHPLUS = 0.D0
00182               FUPLUS = 0.D0
00183             ELSE
00184               UNN= +VNX1*U+VNY1*V
00185               VNN= -VNY1*U+VNX1*V
00186 !
00187               A=MIN(RA3,MAX(-RA3,-UNN/RH))
00188               A2 =A * A
00189               A3 =A2 * A
00190               ALPHA0=ALP*(RA3-A)
00191               ALPHA1=ALP2*(RA32-A2)
00192               ALPHA2=ALP3*(RA33-A3)
00193 !
00194               FHPLUS = H*UNN*ALPHA0 + HRH*ALPHA1
00195               FUPLUS = UNN*(FHPLUS+HRH*ALPHA1) + H*H*ALPHA2
00196             ENDIF
00197 !
00198 !
00199 !    CALCULATION OF FICTIVE STATE (HG,UG,VG)
00200 !
00201 !
00202 !    H GIVEN
00203 !    ########
00204 !
00205             IF(LIMPRO(K,1).EQ.KDIR) THEN
00206 !
00207               C   = SG*RH
00208               VP1 = UNN
00209               VP2 = VP1  + C
00210               VP3 = VP1  - C
00211 !
00212               HG     =HBOR(K)
00213               RHG    =SQRT (HG)
00214               HRHG   =RHG*HG
00215 !
00216               IF (VP2*VP3.LE.0.D0.OR. VP1.LE.0.D0) THEN
00217 !
00218                 IF(HG.EQ.0.D0) THEN
00219                   UG=0.D0
00220                   VG=0.D0
00221                   FHMOINS = 0.D0
00222                   FUMOINS = 0.D0
00223                   SIGMAX=1.D-2
00224                 ELSE
00225 !
00226 !   FLUVIAL REGIME
00227 !   --------------
00228 !
00229                   IF (VP2*VP3.LE.0.D0) THEN
00230 !
00231                     UG=UNN+2.D0*SG*(RH-RHG)
00232                     VG=VNN
00233 !
00234 !   TORRENTIAL REGIME
00235 !   -----------------
00236 !
00237                   ELSE
00238 !
00239 !  IMPOSED FLUX
00240 !  -----------
00241                     IF(LIMPRO(K,2).EQ.KDIR) THEN
00242 !
00243                       UIN = UBOR(K)
00244                       VIN = VBOR(K)
00245                       HUIN = H*UIN
00246                       HVIN = H*VIN
00247 !
00248                       DEST=HUIN*VNX1+HVIN*VNY1
00249                       RVG =-HUIN*VNY1+HVIN*VNX1
00250 !
00251                       A1 = DEST-FHPLUS
00252                       CA1= SQ2*A1/(SG*HG*RHG)
00253                       CALL ZEROPHI(-1.D0,AM,NIT,CA1)
00254 !
00255                       UG= AM*SG*RHG
00256                       VG=RVG/HG
00257 !
00258                     ELSE
00259 !
00260 !  ONE DATUM IS MISSING, WE SUPPOSE "THE LAKE AT REST"
00261 
00262                       UG= 0.D0
00263                       VG= 0.D0
00264 !
00265                     ENDIF
00266 !
00267                   ENDIF
00268 !
00269                   GOTO 220
00270                 ENDIF
00271                 GOTO 200
00272 !
00273 ! THE OUTFLOW IS TORRENTIAL SO WE HAVE NO NEED FOR THE GIVEN H
00274 !
00275               ELSE
00276                 GOTO 100
00277               ENDIF
00278 !
00279 !
00280 !    GIVEN VELOCITY
00281 !    ################
00282 !
00283             ELSE IF(LIMPRO(K,2).EQ.KDIR) THEN
00284 !
00285               UIN = UBOR(K)
00286               VIN = VBOR(K)
00287               HUIN = H*UIN
00288               HVIN = H*VIN
00289 !
00290               DEST=HUIN*VNX1+HVIN*VNY1
00291               RVG =-HUIN*VNY1+HVIN*VNX1
00292 !             WARNING: SIGN CHANGE / INRIA REPORT
00293               A1 = -DEST+FHPLUS
00294               A2 = -UNN - 2.D0*SG*RH
00295 !
00296               IF (A1.LE.0.D0) THEN
00297 !
00298 !               FH- =-A1 CANNOT BE SATISFIED
00299 !
00300                 FHMOINS = 0.D0
00301                 FUMOINS = 0.D0
00302                 VG=0.D0
00303                 SIGMAX=1.E-2
00304               ELSE
00305                 CA1= 1.D0/(G*SQ2*A1)**(1.D0/3.D0)
00306                 CALL ZEROPSI(-0.5D0,AM,NIT,CA1,A2)
00307 !
00308                 RHG =A2/(SG*(AM-2.D0))
00309                 HG= RHG * RHG
00310                 HRHG= RHG * HG
00311 !
00312                 IF (HG.EQ.0.D0) THEN
00313                   UG=0.D0
00314                   VG=0.D0
00315                   FHMOINS = 0.D0
00316                   FUMOINS = 0.D0
00317                   SIGMAX=1.D-2
00318                 ELSE
00319                   UG=-AM*A2/(AM-2.D0)
00320                   VG=RVG/HG
00321                   GOTO 220
00322                 ENDIF
00323               ENDIF
00324               GOTO 200
00325 !
00326 ! NO CONDITION
00327 !
00328             ELSE
00329 !
00330               GOTO 100
00331 !
00332             ENDIF
00333             GOTO 1000
00334 !
00335 !
00336 !   CALCULATION OF F-(HG,UG,VG)
00337 !
00338  220        CONTINUE
00339 !
00340             A=MIN(RA3,MAX(-RA3,-UG/RHG))
00341             A2 =A * A
00342             A3 =A2 * A
00343             ALPHA0=ALP*(A+RA3)
00344             ALPHA1=ALP2*(A2-RA32)
00345             ALPHA2=ALP3*(A3+RA33)
00346 !
00347             FHMOINS = HG*UG*ALPHA0 + HRHG*ALPHA1
00348             FUMOINS = UG*(FHMOINS + HRHG*ALPHA1)
00349      &      + HG*HG*ALPHA2
00350 !
00351             SIGMAX= RHG
00352             UNORM=SQRT(UG *UG + VG*VG)
00353             SIGMAX=MAX( 1.D-2, RA3 *SIGMAX +UNORM )
00354 !
00355 !           CALCUL DES FLUX ET ROTATION INVERSE
00356 !
00357  200        CONTINUE
00358             FLUH=(FHPLUS +FHMOINS)*VNL
00359             FLUU=(FUPLUS +FUMOINS)*VNL
00360 !
00361             IF (FLUH.GE.0.D0) THEN
00362               FLUV= VNN*FLUH
00363             ELSE
00364               FLUV= VG*FLUH
00365             ENDIF
00366 !
00367             FLUTMP=FLUU
00368             FLUU = +VNX1*FLUTMP-VNY1*FLUV
00369             FLUV = +VNY1*FLUTMP+VNX1*FLUV
00370 !
00371 !           CORRECTION OF THE TIME STEP
00372 !
00373             DTL = CFL*DTHAUT(IS)/SIGMAX
00374             DT  = MIN(DT, DTL)
00375 !
00376             GOTO 1000
00377 100         CONTINUE
00378             RUN     = H*UNN
00379 !
00380             FLUH =  RUN* VNL
00381             FLUU =  (U *RUN + 0.5D0*G*H**2* VNX)*VNL
00382             FLUV =  (V *RUN + 0.5D0*G*H**2* VNY)*VNL
00383 !
00384 1000        CONTINUE
00385           ENDIF
00386 !
00387           IF(LIMPRO(K,1).EQ.KDIR)  FLUSORT = FLUSORT + FLUH
00388           IF(LIMPRO(K,2).EQ.KDIR)  FLUENT = FLUENT +FLUH
00389 !RA
00390           FLBOR%R(K)=FLUH
00391 !
00392           CE(IS,1)  = CE(IS,1) - FLUH
00393           CE(IS,2)  = CE(IS,2) - FLUU
00394           CE(IS,3)  = CE(IS,3) - FLUV
00395 !
00396           IF(NTRAC.GT.0) THEN
00397             DO ITRAC=1,NTRAC
00398               FLUHBTEMP%ADR(ITRAC)%P%R(K)=FLUH
00399             ENDDO
00400           ENDIF
00401 !
00402         ENDDO
00403       ENDIF
00404       IF(NCSIZE.GT.1)DT=P_DMIN(DT)
00405 !
00406 !-----------------------------------------------------------------------
00407 !
00408       RETURN
00409       END

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