thomps.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\thomps.f
00002 !
00081                      SUBROUTINE THOMPS
00082 !                    *****************
00083 !
00084      &(HBOR,UBOR,VBOR,TBOR,U,V,H,T,ZF,X,Y,NBOR,FRTYPE,C,
00085      & UCONV,VCONV,T7,XCONV,YCONV,LIHBOR,LIUBOR,LIVBOR,
00086      & LITBOR,IT1,ITRAV2,W1R,W2R,W3R,
00087      & HBTIL,UBTIL,VBTIL,TBTIL,ZBTIL,SURDET,IKLE,IFABOR,NELEM,
00088      & MESH,XNEBOR,YNEBOR,NPOIN,NPTFR,DT,GRAV,
00089      & NTRAC,NFRLIQ,KSORT,KINC,KENT,KENTU,LV,MSK,MASKEL,
00090      & NELMAX,IELM,SHPP,NUMLIQ,SHP,
00091      & DX_T,DY_T,DZ_T,IT3,IT4,HFIELD,UFIELD,VFIELD,ZS,GZSX,GZSY,
00092      & SHPBUF)
00093 !
00094 !***********************************************************************
00095 ! TELEMAC2D   V6P2                                   21/08/2010
00096 !***********************************************************************
00097 !
00098 !
00099 !
00100 !
00101 !
00102 !
00103 !
00104 !
00105 !
00106 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00107 !| C              |-->| WORK ARRAY: CELERITY OF WAVES
00108 !| DT             |-->| TIME STEP
00109 !| DX_T           |<->| WORK ARRAY
00110 !| DY_T           |<->| WORK ARRAY
00111 !| DZ_T           |<->| WORK ARRAY
00112 !| ELT_T          |<->| WORK ARRAY
00113 !| FRTYPE         |-->| TYPE OF LIQUID BOUNDARIES
00114 !| GRAV           |-->| GRAVITY
00115 !| GZSX           |<--| FREE SURFACE GRADIENT ALONG X
00116 !| GZSY           |<--| FREE SURFACE GRADIENT ALONG Y
00117 !| H              |-->| WATER DEPTH AT TIME N
00118 !| HBOR           |<--| PRESCRIBED DEPTH AT BOUNDARIES
00119 !| HBTIL          |<->| WORK ARRAY, DEPTH AT BOUNDARIES AFTER ADVECTION
00120 !| HFIELD         |<->| WORK ARRAY, DEPTH WITH RELAXATION
00121 !| IELM           |-->| TYPE OF ELEMENT
00122 !| IFABOR         |-->| ELEMENTS BEHIND THE EDGES OF A TRIANGLE
00123 !|                |   | IF NEGATIVE OR ZERO, THE EDGE IS A LIQUID
00124 !|                |   | BOUNDARY
00125 !| IKLE           |-->| CONNECTIVITY TABLE
00126 !| ITRAV2         |-->| INTEGER WORK ARRAY
00127 !| IT3            |<->| INTEGER WORK ARRAY
00128 !| IT3            |<->| INTEGER WORK ARRAY
00129 !| KSORT          |-->| CONVENTION FOR FREE OUTPUT
00130 !| KINC           |-->| CONVENTION FOR INCIDENT WAVE BOUNDARY CONDITION
00131 !| KENT           |-->| CONVENTION FOR LIQUID INPUT WITH PRESCRIBED VALUE
00132 !| KENTU          |-->| CONVENTION FOR LIQUID INPUT WITH PRESCRIBED VELOCITY
00133 !| LIHBOR         |-->| TYPE OF BOUNDARY CONDITIONS ON DEPTH
00134 !| LISPFR         |-->| LIST OF BOUNDARY POINTS TO BE DEALT WITH
00135 !| LITBOR         |-->| TYPE OF BOUNDARY CONDITIONS ON TRACERS
00136 !| LIUBOR         |-->| TYPE OF BOUNDARY CONDITIONS ON U
00137 !| LIVBOR         |-->| TYPE OF BOUNDARY CONDITIONS ON V
00138 !| LV             |-->| VECTOR LENGTH (FOR VECTOR MACHINES)
00139 !| MASKEL         |-->| MASKING OF ELEMENTS
00140 !|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
00141 !| MESH           |-->| MESH STRUCTURE
00142 !| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
00143 !| NBOR           |-->| GLOBAL NUMBER OF BOUNDARY POINTS
00144 !| NELEM          |-->| NUMBER OF ELEMENTS
00145 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00146 !| NFRLIQ         |-->| NUMBER OF LIQUID BOUNDARIES
00147 !| NPOIN          |-->| NUMBER OF POINTS
00148 !| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
00149 !| NTRAC          |-->| NUMBER OF TRACERS
00150 !| NUMLIQ         |-->| LIQUID BOUNDARY NUMBER OF BOUNDARY POINTS
00151 !| SHP            |<--| BARYCENTRIC COORDINATES AT THE FOOT
00152 !|                |   | OF CHARACTERISTICS
00153 !| SHPP           |<--| BARYCENTRIC COORDINATES AT THE FOOT
00154 !|                |   | OF CHARACTERISTICS, FOR BOUNDARY POINTS
00155 !| SURDET         |-->| 1/(DETERMINANT OF ISOPARAMETRIC TRANSFORMATION)
00156 !| T              |-->| BLOCK OF TRACERS AT TIME N
00157 !| T7             |<->| WORK BIEF_OBJ STRUCTURE
00158 !| TBOR           |<--| PRESCRIBED BOUNDARY CONDITION ON TRACER
00159 !| TBTIL          |<--| BLOCK OF WORK ARRAYS, TRACERS AFTER ADVECTION
00160 !| U              |<->| X-COMPONENT OF VELOCITY
00161 !| UBOR           |<--| PRESCRIBED VELOCITY U.
00162 !| UBTIL          |<--| WORK ARRAY, U AT BOUNDARIES AFTER ADVECTION
00163 !| UCONV          |-->| WORK ARRAY: ADVECTION FIELDS
00164 !| UFIELD         |<->| WORK ARRAY, U WITH RELAXATION
00165 !| V              |<->| Y-COMPONENT OF VELOCITY
00166 !| UBOR           |<--| PRESCRIBED VELOCITY V.
00167 !| VBTIL          |<--| WORK ARRAY, V AT BOUNDARIES AFTER ADVECTION
00168 !| VCONV          |-->| WORK ARRAY: ADVECTION FIELDS
00169 !| VFIELD         |<->| WORK ARRAY, V WITH RELAXATION
00170 !| W1R            |<->| WORK ARRAY
00171 !| W2R            |<->| WORK ARRAY
00172 !| W3R            |<->| WORK ARRAY
00173 !| X              |-->| ABSCISSAE OF POINTS IN THE MESH
00174 !| XNEBOR         |-->| X-COMPONENT OF NORMAL AT NODES
00175 !| Y              |-->| ORDINATES OF POINTS IN THE MESH
00176 !| YNEBOR         |-->| Y-COMPONENT OF NORMAL AT NODES
00177 !| ZBTIL          |<--| WORK ARRAY, BOTTOM AT BOUNDARIES AFTER ADVECTION
00178 !| ZF             |-->| ELEVATION OF BOTTOM
00179 !| ZS             |<--| FREE SURFACE
00180 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00181 !
00182       USE BIEF
00183       USE INTERFACE_TELEMAC2D, EX_THOMPS => THOMPS
00184       USE STREAMLINE, ONLY : SCARACT
00185 !
00186       IMPLICIT NONE
00187       INTEGER LNG,LU
00188       COMMON/INFO/LNG,LU
00189 !
00190 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00191 !
00192       INTEGER, INTENT(IN) :: NPTFR,NPOIN,NELEM,NELMAX,NFRLIQ,LV
00193       INTEGER, INTENT(IN) :: IELM,NTRAC
00194       INTEGER, INTENT(IN) :: KSORT,KINC,KENT,KENTU
00195       INTEGER, INTENT(IN) :: NBOR(NPTFR)
00196       INTEGER, INTENT(IN) :: IKLE(NELMAX,3),IFABOR(*)
00197       INTEGER, INTENT(IN) :: LIHBOR(NPTFR),LIUBOR(NPTFR),LIVBOR(NPTFR)
00198       INTEGER, INTENT(IN) :: FRTYPE(NFRLIQ),NUMLIQ(NPTFR)
00199 !                                           2*NPTFR REALLY USED HERE
00200       INTEGER, INTENT(INOUT),TARGET  :: IT1(NPOIN)
00201 !     ITRAV2 : OF DIMENSION NPOIN
00202       INTEGER, INTENT(INOUT) :: ITRAV2(*)
00203       LOGICAL, INTENT(IN) :: MSK
00204       DOUBLE PRECISION, INTENT(INOUT) :: HBOR(NPTFR)
00205       DOUBLE PRECISION, INTENT(INOUT) :: UBOR(NPTFR),VBOR(NPTFR)
00206       DOUBLE PRECISION, INTENT(IN)    :: X(NPOIN),Y(NPOIN)
00207       DOUBLE PRECISION, INTENT(IN)    :: XNEBOR(NPTFR),YNEBOR(NPTFR)
00208       DOUBLE PRECISION, INTENT(IN)    :: SURDET(*)
00209       DOUBLE PRECISION, INTENT(IN)    :: GRAV,DT
00210       DOUBLE PRECISION, INTENT(INOUT) :: W1R(NPTFR),W2R(NPTFR)
00211       DOUBLE PRECISION, INTENT(INOUT) :: W3R(NPTFR)
00212       DOUBLE PRECISION, INTENT(INOUT) :: SHPP(3,NPTFR),SHP(3,*)
00213       DOUBLE PRECISION, INTENT(INOUT) :: DX_T(NPTFR),DY_T(NPTFR)
00214       DOUBLE PRECISION, INTENT(INOUT) :: DZ_T(NPTFR)
00215       TYPE(BIEF_OBJ), INTENT(IN)      :: MASKEL
00216       TYPE(BIEF_OBJ), INTENT(INOUT)   :: HBTIL,UBTIL,VBTIL,ZBTIL,T7
00217       TYPE(BIEF_OBJ), INTENT(INOUT)   :: XCONV,YCONV
00218       TYPE(BIEF_OBJ), INTENT(INOUT)   :: UCONV,VCONV,C,U,V
00219       TYPE(BIEF_OBJ), INTENT(INOUT)   :: H,T,TBOR,TBTIL,IT3,IT4
00220       TYPE(BIEF_OBJ), INTENT(INOUT)   :: HFIELD,UFIELD,VFIELD,ZS
00221       TYPE(BIEF_OBJ), INTENT(INOUT)   :: GZSX,GZSY,SHPBUF
00222       TYPE(BIEF_OBJ), INTENT(IN)      :: ZF,LITBOR
00223       TYPE(BIEF_MESH), INTENT(INOUT)  :: MESH
00224 !
00225 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00226 !
00227       INTEGER, DIMENSION(:), POINTER  :: LISPFR
00228       INTEGER, DIMENSION(:),POINTER :: ELT_T
00229       INTEGER K,NPT,J,ITRAC,N,NOMB,NDP,NPLAN,IELMU,ETA(1),I,IFR,IPT
00230       INTEGER FREBUF(1)
00231       DOUBLE PRECISION HMIN,HHBOR,DETADX,DETADY,TBAR(100),TT(100)
00232       DOUBLE PRECISION UCSI,UCSIBAR,UETA,UETABAR,CBAR,HH,TETA
00233       DOUBLE PRECISION ZSTAR(1),ZCONV(1,1),SHZ(1),Z(1,1),UNORM,NORMZS
00234       INTEGER IELEM,SIZEBUF
00235 !
00236       DATA HMIN  /2.D-2/
00237 !
00238       LOGICAL INIT
00239       DATA    INIT/.TRUE./
00240       TYPE(BIEF_OBJ) :: FNCAR1,FTILD1
00241       SAVE
00242 !
00243 !-----------------------------------------------------------------------
00244 !
00245 !     SLICING IT1
00246 !
00247       I=MAX(1,NPTFR)
00248       LISPFR=>IT1(1:I)
00249       ELT_T =>IT1(I+1:2*I)
00250 !
00251 !     SIZE OF BUFFER FOR SHP
00252 !
00253       SIZEBUF=SHPBUF%DIM1/3
00254 !
00255 !-----------------------------------------------------------------------
00256 !
00257       IF(INIT) THEN
00258         CALL ALLBLO(FNCAR1,'FNCAR1')
00259         CALL ALLBLO(FTILD1,'FTILD1')
00260         CALL ADDBLO(FNCAR1,UFIELD)
00261         CALL ADDBLO(FNCAR1,VFIELD)
00262         CALL ADDBLO(FNCAR1,HFIELD)
00263         CALL ADDBLO(FNCAR1,ZF)
00264         CALL ADDBLO(FTILD1,UBTIL)
00265         CALL ADDBLO(FTILD1,VBTIL)
00266         CALL ADDBLO(FTILD1,HBTIL)
00267         CALL ADDBLO(FTILD1,ZBTIL)
00268         IF(NTRAC.GT.0) THEN
00269           DO ITRAC=1,NTRAC
00270             CALL ADDBLO(FNCAR1,T%ADR(ITRAC)%P)
00271             CALL ADDBLO(FTILD1,TBTIL%ADR(ITRAC)%P)
00272           ENDDO
00273         ENDIF
00274         INIT=.FALSE.
00275       ENDIF
00276 !
00277       ETA(1)=1
00278 !     NDP ALSO FIRST DIMENSION OF SHP AND SHPP
00279       NDP=3
00280       NPLAN=1
00281       IELMU=IELM
00282 !
00283 !-----------------------------------------------------------------------
00284 !
00285 !     CREATING FIELDS OF H, U AND V, THAT CONTAIN THE PRESCRIBED DATA
00286 !     WITH RELAXATION TETA (IF 1.D0, THE ROUGH DATA ARE TAKEN ON THE
00287 !     THOMPSON BOUNDARIES FOR INTERPOLATION AT THE FOOT OF
00288 !     CHARACTERISTICS)
00289 !
00290       CALL OS('X=Y     ',X=HFIELD,Y=H)
00291       CALL OS('X=Y     ',X=UFIELD,Y=U)
00292       CALL OS('X=Y     ',X=VFIELD,Y=V)
00293 !
00294 !     FORCING DISCRETISATION TO IELM
00295 !
00296       UFIELD%ELM=IELM
00297       VFIELD%ELM=IELM
00298       HFIELD%ELM=IELM
00299 !
00300       UFIELD%DIM1=BIEF_NBPTS(IELM,MESH)
00301       VFIELD%DIM1=BIEF_NBPTS(IELM,MESH)
00302       HFIELD%DIM1=BIEF_NBPTS(IELM,MESH)
00303 !
00304       TETA=1.D0
00305 !     CAN BE RELAXED, TESTED UP TO 0.05 IN PALUEL BOX MODEL
00306 !     TETA=0.05D0
00307 !
00308       DO K=1,NPTFR
00309         IF(NUMLIQ(K).NE.0) THEN
00310           IF(FRTYPE(NUMLIQ(K)).EQ.2) THEN
00311             N=NBOR(K)
00312             IF(LIHBOR(K).EQ.KENT) THEN
00313               HFIELD%R(N)=TETA*HBOR(K)+(1.D0-TETA)*HFIELD%R(N)
00314             ENDIF
00315             IF(LIUBOR(K).EQ.KENT.OR.LIUBOR(K).EQ.KENTU) THEN
00316               UFIELD%R(N)=TETA*UBOR(K)+(1.D0-TETA)*UFIELD%R(N)
00317             ENDIF
00318             IF(LIVBOR(K).EQ.KENT.OR.LIVBOR(K).EQ.KENTU) THEN
00319               VFIELD%R(N)=TETA*VBOR(K)+(1.D0-TETA)*VFIELD%R(N)
00320             ENDIF
00321           ENDIF
00322         ENDIF
00323       ENDDO
00324 !
00325 !-----------------------------------------------------------------------
00326 !
00327 !     COMPUTES THE CELERITY
00328 !
00329       CALL OS('X=CY    ',C,H,H,GRAV)
00330       CALL CLIP(C,0.D0,.TRUE.,1.D6,.FALSE.,0)
00331       CALL OS('X=SQR(Y)',X=C,Y=C)
00332 !
00333 !     COMPUTES MINUS THE FREE SURFACE GRADIENT
00334 !
00335       CALL OS('X=Y+Z   ',X=ZS,Y=H,Z=ZF)
00336       CALL VECTOR(GZSX,'=','GRADF          X',IELM,
00337      &            -1.D0,ZS,ZS,ZS,ZS,ZS,ZS,MESH,MSK,MASKEL)
00338       CALL VECTOR(GZSY,'=','GRADF          Y',IELM,
00339      &            -1.D0,ZS,ZS,ZS,ZS,ZS,ZS,MESH,MSK,MASKEL)
00340       IF(NCSIZE.GT.1) THEN
00341         CALL PARCOM(GZSX,2,MESH)
00342         CALL PARCOM(GZSY,2,MESH)
00343       ENDIF
00344 !
00345 !     REGROUPS THE POINTS WITH CHARACTERISTICS TREATMENT BY THIS
00346 !     PROCESSOR AND BUILDS THEIR BARYCENTRIC COORDINATES
00347 !     NPT: NUMBER OF POINTS
00348 !     LISPFR: LIST OF THOSE POINTS (BOUNDARY NODE NUMBERS)
00349 !
00350       NPT=0
00351       DO K=1,NPTFR
00352         IF(NUMLIQ(K).NE.0) THEN
00353           N=NBOR(K)
00354           IF(FRTYPE(NUMLIQ(K)).EQ.2.AND.MESH%ELTCAR%I(N).NE.0) THEN
00355             NPT=NPT+1
00356             LISPFR(NPT)=K
00357             IELEM= MESH%ELTCAR%I(N)
00358             IF(IKLE(IELEM,1).EQ.N) THEN
00359               SHP(1,N)=1.D0
00360               SHP(2,N)=0.D0
00361               SHP(3,N)=0.D0
00362             ENDIF
00363             IF(IKLE(IELEM,2).EQ.N) THEN
00364               SHP(1,N)=0.D0
00365               SHP(2,N)=1.D0
00366               SHP(3,N)=0.D0
00367             ENDIF
00368             IF(IKLE(IELEM,3).EQ.N) THEN
00369               SHP(1,N)=0.D0
00370               SHP(2,N)=0.D0
00371               SHP(3,N)=1.D0
00372             ENDIF
00373           ENDIF
00374         ENDIF
00375       ENDDO
00376 !
00377 !--------------------------------------------------------------------
00378 !     CHARACTERISTICS WITH ADVECTION FIELD U
00379 !     MAY BE CALLED WITH NPT=0 IF OTHER PROCESSORS STILL AT WORK
00380 !--------------------------------------------------------------------
00381 !
00382 !     COMPUTES THE ADVECTION FIELD U
00383 !
00384       DO N=1,NPOIN
00385         UCONV%R(N)=U%R(N)
00386         VCONV%R(N)=V%R(N)
00387       ENDDO
00388 !
00389 !     PREPARING THE CALL TO SCARACT
00390 !
00391       DO IFR=1,NPT
00392         IPT=NBOR(LISPFR(IFR))
00393         XCONV%R(IFR) = X(IPT)
00394         YCONV%R(IFR) = Y(IPT)
00395         ELT_T(IFR)   = MESH%ELTCAR%I(IPT)
00396         SHPP(1,IFR)= SHP(1,IPT)
00397         SHPP(2,IFR)= SHP(2,IPT)
00398         SHPP(3,IFR)= SHP(3,IPT)
00399       ENDDO
00400 !
00401       NOMB=4+NTRAC
00402       CALL SCARACT(FNCAR1,FTILD1,UCONV%R,VCONV%R,VCONV%R,VCONV%R,
00403      &             X,Y,
00404      &             ZSTAR,ZSTAR,XCONV%R,YCONV%R,ZCONV,ZCONV,
00405      &             DX_T,DY_T,DZ_T,DZ_T,Z,
00406      &             SHPP,SHZ,SHZ,
00407      &             SURDET,DT,IKLE,IFABOR,ELT_T,
00408      &             ETA,ETA,IT3%I,IT4%I,IELM,IELMU,NELEM,NELMAX,
00409      &             NOMB,NPOIN,NPOIN,NDP,NPLAN,1,MESH,NPT,U%DIM1,-1,
00410      &             SHPBUF%R,SHPBUF%R,SHPBUF%R,FREBUF,SIZEBUF)
00411 !
00412 !----------------------------------------------------------------------
00413 !     UBTIL, VBTIL, HBTIL, TBTIL AT BOUNDARY NODES NUMBERING
00414 !----------------------------------------------------------------------
00415 !
00416 !     BACKWARD LOOP TO AVOID ERASING DATA, K ALWAYS GREATER THAN J
00417       DO J=NPT,1,-1
00418         K=LISPFR(J)
00419         UBTIL%R(K)=UBTIL%R(J)
00420         VBTIL%R(K)=VBTIL%R(J)
00421         HBTIL%R(K)=HBTIL%R(J)
00422 !       ZBTIL%R(K)=ZBTIL%R(J)
00423       ENDDO
00424       IF(NTRAC.GT.0) THEN
00425         DO ITRAC=1,NTRAC
00426           DO J=NPT,1,-1
00427             K=LISPFR(J)
00428             TBTIL%ADR(ITRAC)%P%R(K)=TBTIL%ADR(ITRAC)%P%R(J)
00429           ENDDO
00430         ENDDO
00431       ENDIF
00432 !
00433 !     BEFORE PARCOM_BORD ,CANCELLING VALUES OF POINTS TREATED
00434 !     IN ANOTHER SUB-DOMAIN
00435 !
00436       IF(NCSIZE.GT.1) THEN
00437         DO K=1,NPTFR
00438           IF(MESH%ELTCAR%I(NBOR(K)).EQ.0) THEN
00439             UBTIL%R(K)=0.D0
00440             VBTIL%R(K)=0.D0
00441             HBTIL%R(K)=0.D0
00442 !           ZBTIL%R(K)=0.D0
00443           ENDIF
00444         ENDDO
00445         IF(NTRAC.GT.0) THEN
00446           DO K=1,NPTFR
00447             IF(MESH%ELTCAR%I(NBOR(K)).EQ.0) THEN
00448               DO ITRAC=1,NTRAC
00449                 TBTIL%ADR(ITRAC)%P%R(K)=0.D0
00450               ENDDO
00451             ENDIF
00452           ENDDO
00453         ENDIF
00454         CALL PARCOM_BORD(UBTIL%R,1,MESH)
00455         CALL PARCOM_BORD(VBTIL%R,1,MESH)
00456         CALL PARCOM_BORD(HBTIL%R,1,MESH)
00457 !       CALL PARCOM_BORD(ZBTIL%R,1,MESH)
00458         IF(NTRAC.GT.0) THEN
00459           DO ITRAC=1,NTRAC
00460             CALL PARCOM_BORD(TBTIL%ADR(ITRAC)%P%R,1,MESH)
00461           ENDDO
00462         ENDIF
00463       ENDIF
00464 !
00465 !     COMPUTES THE RIEMANN INVARIANTS W1 AND W4 (SECOND DIMENSION OF TBOR)
00466 !     CARRIED BY THIS FIELD
00467 !
00468 !     IF W1=0 IS OK, THIS PART COULD BE DONE ONLY WHEN NTRAC.GT.0
00469 !
00470       DO K=1,NPTFR
00471         IF(NUMLIQ(K).NE.0) THEN
00472         IF(FRTYPE(NUMLIQ(K)).EQ.2) THEN
00473           N=NBOR(K)
00474           UNORM=SQRT(U%R(N)**2+V%R(N)**2)
00475           IF(UNORM.GT.1.D-12) THEN
00476             DETADX=U%R(N)/UNORM
00477             DETADY=V%R(N)/UNORM
00478           ELSE
00479             DETADX=0.D0
00480             DETADY=0.D0
00481           ENDIF
00482           UCSIBAR=-U%R(N)*DETADY+V%R(N)*DETADX
00483           HH=HBTIL%R(K)
00484           UCSI=-UBTIL%R(K)*DETADY+VBTIL%R(K)*DETADX
00485           W1R(K)=HH*(UCSIBAR-UCSI)
00486           IF(NTRAC.GT.0) THEN
00487             DO ITRAC=1,NTRAC
00488               TBAR(ITRAC)=T%ADR(ITRAC)%P%R(N)
00489             ENDDO
00490             IF(UCONV%R(N)*XNEBOR(K)+VCONV%R(N)*YNEBOR(K).GT.0.D0) THEN
00491 !             VELOCITY EXITING THE DOMAIN, THE CHARACTERISTICS ARE USED
00492               DO ITRAC=1,NTRAC
00493                 TT(ITRAC)=TBTIL%ADR(ITRAC)%P%R(K)
00494               ENDDO
00495             ELSE
00496 !             VELOCITY ENTERING THE DOMAIN, PRESCRIBED VALUES TAKEN
00497               DO ITRAC=1,NTRAC
00498                 TT(ITRAC)=TBOR%ADR(ITRAC)%P%R(K)
00499               ENDDO
00500             ENDIF
00501             DO ITRAC=1,NTRAC
00502 !             W4
00503               TBOR%ADR(ITRAC)%P%R(K+NPTFR)=HH*(TT(ITRAC)-TBAR(ITRAC))
00504             ENDDO
00505           ENDIF
00506         ENDIF
00507         ENDIF
00508       ENDDO
00509 !
00510 !----------------------------------------------------------
00511 !     COMPUTES THE ADVECTION FIELD U + C
00512 !----------------------------------------------------------
00513 !
00514       DO N=1,NPOIN
00515         UNORM=SQRT(U%R(N)**2+V%R(N)**2)
00516         IF(UNORM.GT.1.D-12) THEN
00517           DETADX=U%R(N)/UNORM
00518           DETADY=V%R(N)/UNORM
00519         ELSE
00520 !         VERY IMPORTANT
00521           NORMZS=SQRT(GZSX%R(N)**2+GZSY%R(N)**2)
00522           IF(NORMZS.GT.1.D-12) THEN
00523             DETADX=GZSX%R(N)/NORMZS
00524             DETADY=GZSY%R(N)/NORMZS
00525           ELSE
00526             DETADX=0.D0
00527             DETADY=0.D0
00528           ENDIF
00529         ENDIF
00530         UETABAR=U%R(N)*DETADX+V%R(N)*DETADY+C%R(N)
00531         UCONV%R(N)=UETABAR*DETADX
00532         VCONV%R(N)=UETABAR*DETADY
00533       ENDDO
00534 !
00535 !----------------------------------------------------------------------
00536 !     CHARACTERISTICS FOR THE GROUP OF POINTS, FIELD U + C
00537 !     MAY BE CALLED WITH NPT=0 IF OTHER PROCESSORS STILL AT WORK
00538 !----------------------------------------------------------------------
00539 !
00540 !     PREPARING THE CALL TO SCARACT
00541 !
00542       DO IFR=1,NPT
00543         IPT=NBOR(LISPFR(IFR))
00544         XCONV%R(IFR) = X(IPT)
00545         YCONV%R(IFR) = Y(IPT)
00546         ELT_T(IFR)   = MESH%ELTCAR%I(IPT)
00547         SHPP(1,IFR)  = SHP(1,IPT)
00548         SHPP(2,IFR)  = SHP(2,IPT)
00549         SHPP(3,IFR)  = SHP(3,IPT)
00550       ENDDO
00551 !
00552       NOMB=4
00553       CALL SCARACT(FNCAR1,FTILD1,UCONV%R,VCONV%R,VCONV%R,VCONV%R,
00554      &             X,Y,
00555      &             ZSTAR,ZSTAR,XCONV%R,YCONV%R,ZCONV,ZCONV,
00556      &             DX_T,DY_T,DZ_T,DZ_T,Z,
00557      &             SHPP,SHZ,SHZ,
00558      &             SURDET,DT,IKLE,IFABOR,ELT_T,ETA,ETA,
00559      &             IT3%I,IT4%I,IELM,IELMU,NELEM,NELMAX,NOMB,NPOIN,
00560      &             NPOIN,NDP,NPLAN,1,MESH,NPT,U%DIM1,-1,
00561      &             SHPBUF%R,SHPBUF%R,SHPBUF%R,FREBUF,SIZEBUF)
00562 !
00563 !----------------------------------------------------------------------
00564 !     UBTIL, VBTIL, HBTIL AT BOUNDARY NODES NUMBERING
00565 !----------------------------------------------------------------------
00566 !
00567 !     BACKWARD LOOP TO AVOID ERASING DATA, K ALWAYS GREATER THAN J
00568       DO J=NPT,1,-1
00569         K=LISPFR(J)
00570         UBTIL%R(K)=UBTIL%R(J)
00571         VBTIL%R(K)=VBTIL%R(J)
00572         HBTIL%R(K)=HBTIL%R(J)
00573         ZBTIL%R(K)=ZBTIL%R(J)
00574       ENDDO
00575 !
00576 !     BEFORE PARCOM_BORD ,CANCELLING VALUES OF POINTS TREATED
00577 !     IN ANOTHER SUB-DOMAIN
00578 !
00579       IF(NCSIZE.GT.1) THEN
00580         DO K=1,NPTFR
00581           IF(MESH%ELTCAR%I(NBOR(K)).EQ.0) THEN
00582             UBTIL%R(K)=0.D0
00583             VBTIL%R(K)=0.D0
00584             HBTIL%R(K)=0.D0
00585             ZBTIL%R(K)=0.D0
00586           ENDIF
00587         ENDDO
00588         CALL PARCOM_BORD(UBTIL%R,1,MESH)
00589         CALL PARCOM_BORD(VBTIL%R,1,MESH)
00590         CALL PARCOM_BORD(HBTIL%R,1,MESH)
00591         CALL PARCOM_BORD(ZBTIL%R,1,MESH)
00592       ENDIF
00593 !
00594 !----------------------------------------------------------------------
00595 !     COMPUTES THE RIEMANN INVARIANTS W2 CARRIED BY THIS ADVECTION FIELD
00596 !----------------------------------------------------------------------
00597 !
00598       DO K=1,NPTFR
00599         IF(NUMLIQ(K).NE.0) THEN
00600         IF(FRTYPE(NUMLIQ(K)).EQ.2) THEN
00601           N=NBOR(K)
00602           UNORM=SQRT(U%R(N)**2+V%R(N)**2)
00603           IF(UNORM.GT.1.D-12) THEN
00604             DETADX=U%R(N)/UNORM
00605             DETADY=V%R(N)/UNORM
00606           ELSE
00607             DETADX=0.D0
00608             DETADY=0.D0
00609           ENDIF
00610 !         UETABAR=U%R(N)*DETADX+V%R(N)*DETADY
00611           UETABAR=UNORM
00612           CBAR=C%R(N)
00613           HH  =HBTIL%R(K)
00614           UETA=UBTIL%R(K)*DETADX+VBTIL%R(K)*DETADY
00615           W2R(K)=(HH+ZBTIL%R(K)-ZF%R(N))*CBAR+HH*(UETA-UETABAR)
00616         ENDIF
00617         ENDIF
00618       ENDDO
00619 !
00620 !     COMPUTES THE ADVECTION FIELD U-C
00621 !
00622       DO N=1,NPOIN
00623         UNORM=SQRT(U%R(N)**2+V%R(N)**2)
00624         IF(UNORM.GT.1.D-12) THEN
00625           DETADX=U%R(N)/UNORM
00626           DETADY=V%R(N)/UNORM
00627         ELSE
00628           NORMZS=SQRT(GZSX%R(N)**2+GZSY%R(N)**2)
00629           IF(NORMZS.GT.1.D-12) THEN
00630             DETADX=GZSX%R(N)/NORMZS
00631             DETADY=GZSY%R(N)/NORMZS
00632           ELSE
00633             DETADX=0.D0
00634             DETADY=0.D0
00635           ENDIF
00636         ENDIF
00637         UETABAR=U%R(N)*DETADX+V%R(N)*DETADY-C%R(N)
00638         UCONV%R(N)=UETABAR*DETADX
00639         VCONV%R(N)=UETABAR*DETADY
00640       ENDDO
00641 !
00642 !     CHARACTERISTICS FOR THE GROUP OF POINTS, FIELD U + C
00643 !
00644 !     PREPARING THE CALL TO SCARACT
00645 !
00646       DO IFR=1,NPT
00647         IPT=NBOR(LISPFR(IFR))
00648         XCONV%R(IFR) = X(IPT)
00649         YCONV%R(IFR) = Y(IPT)
00650         ELT_T(IFR)   = MESH%ELTCAR%I(IPT)
00651         SHPP(1,IFR)  = SHP(1,IPT)
00652         SHPP(2,IFR)  = SHP(2,IPT)
00653         SHPP(3,IFR)  = SHP(3,IPT)
00654       ENDDO
00655 !
00656       NOMB=4
00657       CALL SCARACT(FNCAR1,FTILD1,UCONV%R,VCONV%R,VCONV%R,VCONV%R,
00658      &             X,Y,ZSTAR,ZSTAR,
00659      &             XCONV%R,YCONV%R,ZCONV,ZCONV,
00660      &             DX_T,DY_T,DZ_T,DZ_T,Z,SHPP,SHZ,SHZ,
00661      &             SURDET,DT,IKLE,IFABOR,ELT_T,ETA,ETA,
00662      &             IT3%I,IT4%I,IELM,
00663      &             IELMU,NELEM,NELMAX,NOMB,NPOIN,NPOIN,NDP,NPLAN,1,
00664      &             MESH,NPT,U%DIM1,-1,
00665      &             SHPBUF%R,SHPBUF%R,SHPBUF%R,FREBUF,SIZEBUF)
00666 !
00667 !----------------------------------------------------------------------
00668 !     UBTIL, VBTIL, HBTIL AT BOUNDARY NODES NUMBERING
00669 !----------------------------------------------------------------------
00670 !
00671 !     BACKWARD LOOP TO AVOID ERASING DATA, K ALWAYS GREATER THAN J
00672       DO J=NPT,1,-1
00673         K=LISPFR(J)
00674         UBTIL%R(K)=UBTIL%R(J)
00675         VBTIL%R(K)=VBTIL%R(J)
00676         HBTIL%R(K)=HBTIL%R(J)
00677         ZBTIL%R(K)=ZBTIL%R(J)
00678       ENDDO
00679 !
00680 !     BEFORE PARCOM_BORD ,CANCELLING VALUES OF POINTS TREATED
00681 !     IN ANOTHER SUB-DOMAIN
00682 !
00683       IF(NCSIZE.GT.1) THEN
00684         DO K=1,NPTFR
00685           IF(MESH%ELTCAR%I(NBOR(K)).EQ.0) THEN
00686             UBTIL%R(K)=0.D0
00687             VBTIL%R(K)=0.D0
00688             HBTIL%R(K)=0.D0
00689             ZBTIL%R(K)=0.D0
00690           ENDIF
00691         ENDDO
00692         CALL PARCOM_BORD(UBTIL%R,1,MESH)
00693         CALL PARCOM_BORD(VBTIL%R,1,MESH)
00694         CALL PARCOM_BORD(HBTIL%R,1,MESH)
00695         CALL PARCOM_BORD(ZBTIL%R,1,MESH)
00696       ENDIF
00697 !
00698 ! COMPUTES THE RIEMANN INVARIANTS W3 CARRIED BY THIS ADVECTION FIELD
00699 !
00700       DO K=1,NPTFR
00701         IF(NUMLIQ(K).NE.0) THEN
00702         IF(FRTYPE(NUMLIQ(K)).EQ.2) THEN
00703           N=NBOR(K)
00704           UNORM=SQRT(U%R(N)**2+V%R(N)**2)
00705           IF(UNORM.GT.1.D-12) THEN
00706             DETADX=U%R(N)/UNORM
00707             DETADY=V%R(N)/UNORM
00708           ELSE
00709             DETADX=0.D0
00710             DETADY=0.D0
00711           ENDIF
00712 !         UETABAR=U%R(N)*DETADX+V%R(N)*DETADY
00713 !         MARCHE AUSSI
00714           UETABAR=UNORM
00715           CBAR=C%R(N)
00716           HH  =HBTIL%R(K)
00717           UETA=UBTIL%R(K)*DETADX+VBTIL%R(K)*DETADY
00718           W3R(K)=(HH+ZBTIL%R(K)-ZF%R(N))*CBAR+HH*(-UETA+UETABAR)
00719         ENDIF
00720         ENDIF
00721       ENDDO
00722 !
00723 !----------------------------------------------------------------------
00724 !       RE-BUILDS THE TELEMAC-2D VARIABLES
00725 !----------------------------------------------------------------------
00726 !
00727       DO K=1,NPTFR
00728         IF(NUMLIQ(K).NE.0) THEN
00729         IF(FRTYPE(NUMLIQ(K)).EQ.2) THEN
00730           N=NBOR(K)
00731           UNORM=SQRT(U%R(N)**2+V%R(N)**2)
00732           IF(UNORM.GT.1.D-12) THEN
00733             DETADX=U%R(N)/UNORM
00734             DETADY=V%R(N)/UNORM
00735           ELSE
00736             NORMZS=SQRT(GZSX%R(N)**2+GZSY%R(N)**2)
00737             IF(NORMZS.GT.1.D-12) THEN
00738               DETADX=GZSX%R(N)/NORMZS
00739               DETADY=GZSY%R(N)/NORMZS
00740             ELSE
00741               DETADX=0.D0
00742               DETADY=0.D0
00743             ENDIF
00744           ENDIF
00745           IF(C%R(N)**2.GT.GRAV*HMIN) THEN
00746             HBOR(K)=(W2R(K)+W3R(K))/(2*C%R(N))
00747             IF(HBOR(K).GT.HMIN) THEN
00748 !             BEWARE TIDAL FLATS, AND HIDDEN PARAMETER 0.1
00749               HHBOR=MAX(0.1D0,HBOR(K))
00750 !             formulas 37 and 38 of proceedings but coefficient 0.5
00751 !             is missing in this document.
00752               UBOR(K)=(       DETADY* W1R(K)
00753      &                 +0.5D0*DETADX*(W2R(K)-W3R(K)))/HHBOR+U%R(N)
00754               VBOR(K)=(      -DETADX* W1R(K)
00755      &                 +0.5D0*DETADY*(W2R(K)-W3R(K)))/HHBOR+V%R(N)
00756               IF(NTRAC.GT.0) THEN
00757                 DO ITRAC=1,NTRAC
00758 !                 REMEMBER THAT W4 IS STORED INTO SECOND DIMENSION OF TBOR
00759                   TBOR%ADR(ITRAC)%P%R(K)=
00760      &            TBOR%ADR(ITRAC)%P%R(K+NPTFR)/HHBOR+T%ADR(ITRAC)%P%R(N)
00761                 ENDDO
00762               ENDIF
00763             ELSE
00764 !             BECOMES DRY
00765               HBOR(K)=MAX(0.D0,HBOR(K))
00766               UBOR(K)=0.D0
00767               VBOR(K)=0.D0
00768             ENDIF
00769           ELSE
00770 !           WAS DRY, H IS GIVEN BY BORD
00771             UBOR(K)=0.D0
00772             VBOR(K)=0.D0
00773           ENDIF
00774         ENDIF
00775         ENDIF
00776       ENDDO
00777 !
00778 !-----------------------------------------------------------------------
00779 !
00780       RETURN
00781       END

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