bedload_seccurrent.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\sisyphe\bedload_seccurrent.f
00002 !
00089                      SUBROUTINE BEDLOAD_SECCURRENT
00090 !                    *****************************
00091 !
00092      &(IELMU)
00093 !
00094 !***********************************************************************
00095 ! SISYPHE   V6P2                                   21/07/2011
00096 !***********************************************************************
00097 !
00098 !
00099 !
00100 !
00101 !
00102 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00103 !| IELMU          |-->| TYPE OF ELEMENT FOR VELOCITY
00104 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00105 !
00106       USE DECLARATIONS_SISYPHE
00107       USE BIEF
00108       IMPLICIT NONE
00109 !
00110       INTEGER LNG,LU
00111       COMMON/INFO/LNG,LU
00112 !
00113 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00114 !
00115       INTEGER, INTENT(IN) :: IELMU
00116 !
00117 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00118 !
00119       INTEGER I
00120       DOUBLE PRECISION C, ALPHAL
00121 !
00122 !     RK MODIFICATION FOR SECONDARY CURRENTS
00123 !
00124 !     COMPUTES THE GRADIENT OF THE FREE SURFACE IN X AND Y DIRECTION
00125 !
00126       CALL VECTOR(T5,'=','GRADF          X',IELMU,
00127      &            1.D0,Z,S,S,S,S,S,MESH,MSK,MASKEL)
00128       CALL VECTOR(T6,'=','GRADF          Y',IELMU,
00129      &            1.D0,Z,S,S,S,S,S,MESH,MSK,MASKEL)
00130 !
00131 !     FOR PARALLEL COMPUTING
00132 !
00133       IF(NCSIZE.GT.1) THEN
00134         CALL PARCOM (T5, 2, MESH)
00135         CALL PARCOM (T6, 2, MESH)
00136       ENDIF
00137 !
00138       CALL OS('X=XY    ',X=T5,Y=UNSV2D)
00139       CALL OS('X=XY    ',X=T6,Y=UNSV2D)
00140 !
00141 !     COMPUTES THE X- AND Y-COMPONENTS OF THE SECONDARY CURRENT
00142 !     ACCORDING TO ENGELUNG. TAU_X_SEC = C*QV, TAU_Y_SEC = C*QU
00143 !
00144 ! AT THE MOMENT ALPHA MUST BE SET HERE  (0.75 FOR VERY ROUGH BOTTOMS, 1 FOR SMOOTH ONES)
00145 ! BEWARE: THE VARIABLE ALPHA IS MORE THAN THE ALPHA FROM THE THEORY
00146 ! Now included as Keyword in CAS - FILE
00147       ALPHAL = ALPHA
00148 !RK
00149       IF(ALPHA.GT.0.D0) THEN
00150         ALPHAL = 7.D0 / ALPHAL * XMVE *GRAV
00151       ELSE
00152         ALPHAL = 0.D0
00153       ENDIF
00154 !
00155       CALL OS( 'X=YZ    ' , T1 , T6      , U2D   , C   ) ! DZSDY*U2D
00156       CALL OS( 'X=YZ    ' , T2 , T5      , V2D   , C   ) ! DZSDX*V2D
00157       CALL OS( 'X=-Y    ' , T2 , T2      , T3   , C   )
00158       CALL OS( 'X=X+Y   ' , T1 , T2      , T3   , C   ) ! U2D*DZSDY - V2D*DZSDX
00159 ! NOTE JMH : WHY NOT THE FOLLOWING LINE INSTEAD OF THE 2 PREVIOUS ???
00160 !     CALL OS( 'X=X-Y   ',X=T1,Y=T2) ! U2D*DZSDY - V2D*DZSDX
00161 !
00162       CALL OS( 'X=YZ    ' , T2 , U2D      , U2D   , C   ) ! U2D**2
00163       CALL OS( 'X=YZ    ' , T3 , V2D      , V2D   , C   ) ! V2D**2
00164       CALL OS( 'X=X+Y   ' , T2 , T3      , T3   , C   ) ! U2D**2+V2D**2
00165 !
00166       CALL OS('X=Y/Z   ' , T1 , T1 , T2, C ,2 , 0.D0,1.D-12) !(U2D*DZSDY - V2D*DZSDX)/(U2D**2+V2D**2)
00167 !
00168       CALL OS( 'X=CX    ' , T1 , T2      , T3   , ALPHAL   ) ! T1 * 7/ALPHAL*XMVE*GRAV
00169       CALL OS( 'X=XY    ' , T1 , HN      , T3   , C   ) ! T1*HN
00170 !
00171 !     FOR ALL ROUGHNESS LAWS
00172 !
00173       CALL OS( 'X=CXY   ' , X=T1 , Y=CF , C=0.5D0 )!T1*CF/2
00174 !
00175 !     TAU_X_SEK = -C*V2D : T5
00176 !     TAU_Y_SEK =  C*U2D : T6
00177 !
00178       CALL OS('X=YZ    ' , T5 , T1    , V2D,  C ) ! C*V2D
00179       CALL OS('X=YZ    ' , T6 , T1    , U2D,  C ) ! C*U2D
00180       CALL OS('X=-Y    ' , T6 , T6    , QV,   C ) ! -C*U2D
00181 !
00182 !     SQRT(TAU_X_SEK**2+TAU_Y_SEK**2) : T3
00183 !
00184       CALL OS('X=YZ    ' , T2 , T5    , T5,  C ) ! T2 = (C*V2D)**2
00185       CALL OS('X=YZ    ' , T3 , T6    , T6,  C ) ! T3 = (C*U2D)**2
00186       CALL OS('X=X+Y   ' , T2 , T3    , T3,  C ) ! T2 = (C*V2D)**2+(C*U2D)**2
00187       CALL OS('X=SQR(Y)' , X=T3 , Y=T2 ) ! T3 = SQRT((C*U2D)**2+(C*V2D)**2
00188 !
00189 !     TAU_X_GES = TOB*EFFPNT*CALFA + TAU_X_SEK : T1
00190 !     TAU_Y_GES = TOB*EFFPNT*SALFA + TAU_Y_SEK : T2
00191 !
00192       CALL OS( 'X=YZ    ' , T1 , TOB      , COEFPN   , C   ) ! TOB*EFFPNT
00193       CALL OS( 'X=YZ    ' , T2 , T1      ,  SALFA   , C   ) ! TOB*EFFPNT*SALFA
00194       CALL OS( 'X=YZ    ' , T1 , T1      , CALFA   , C   ) ! TOB*EFFPNT*CALFA
00195       CALL OS('X=X+Y   ' , T1 , T5    , T3,  C ) ! TAU_X_GES = TOB*CALFA+TAU_X_SEK
00196       CALL OS('X=X+Y   ' , T2 , T6    , T3,  C ) ! TAU_Y_GES = TOB*SALFA+TAU_Y_SEK
00197 !
00198 !     TAU_GES=SQRT(TAU_X_GES**2+TAU_Y_GES**2)
00199 !
00200       CALL OS( 'X=YZ    ' , T3 , T1 , T1 , C   ) ! TAU_X_GES**2
00201       CALL OS( 'X=YZ    ' , T4 , T2 , T2 , C   ) ! TAU_Y_GES**2
00202       CALL OS('X=X+Y   ' , T4 , T3    , T3,  C ) ! TAU_X_GES**2+TAU_Y_GES**2
00203       CALL OS('X=SQR(Y)' , T4 , T4    , T3,  C ) ! SQRT(TAU_X_GES**2+TAU_Y_GES**2)
00204 !
00205 !     NEW ANGLE
00206 !     CALFA_NEW = COS(TAU_X_GES/TAU_GES)
00207 !     SALFA_NEW = SIN(TAU_Y_GES/TAU_GES)
00208 !
00209       CALL OS('X=Y/Z   ' , T1 , T1 , T4, C ,2 , 0.D0,1.D-12) !TAU_X_GES/TAU_GES
00210       CALL OS('X=Y/Z   ' , T2 , T2 , T4, C ,2 , 0.D0,1.D-12) !TAU_Y_GES/TAU_GES
00211 !
00212 !     TAKEN FROM EFFPNT UEBER
00213 !     TO MAKE SURE THAT TAU_X_GES/TAU_GES IS IN RANGE [-1,1]
00214 !
00215       DO I=1,NPOIN
00216         IF(T1%R(I).LT.-1.D0.OR.T1%R(I).GT.1.D0.OR.
00217      &     T2%R(I).LT.-1.D0.OR.T2%R(I).GT.1.D0) THEN
00218           WRITE(LU,*) 'NOT ACCEPTABLE BORDER CROSSING',I
00219         ENDIF
00220         T1%R(I) = MIN(T1%R(I),1.D0)
00221         T1%R(I) = MAX(T1%R(I),-1.D0)
00222         T2%R(I) = MIN(T2%R(I),1.D0)
00223         T2%R(I) = MAX(T2%R(I),-1.D0)
00224       ENDDO
00225 !
00226 !     COEFPN_NEW = TAU_GES / TOB
00227 !
00228       CALL OS('X=Y/Z   ' , COEFPN , T4 , TOB, C ,2 , 0.D0,1.D-12) !COEFPN=TAU_GES/TOB
00229       CALL OS( 'X=Y     ' ,X=CALFA ,Y=T1 ) ! (TAU_X_GES/TAU_GES)
00230       CALL OS( 'X=Y     ' ,X=SALFA ,Y=T2 ) ! (TAU_Y_GES/TAU_GES)
00231 !
00232 !     FROM EFFPNT
00233 !     FOR BOUNDARY NODES WITH IMPOSED FLOW :
00234 !     QS IS NOT MODIFIED WHEN USER-DEFINED
00235 !
00236       DO I = 1 , NPTFR
00237         IF(LIQBOR%I(I).EQ.5) THEN
00238           COEFPN%R(MESH%NBOR%I(I)) = 1.D0
00239         ENDIF
00240       ENDDO
00241 !
00242 !-----------------------------------------------------------------------
00243 !
00244       RETURN
00245       END

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