The TELEMAC-MASCARET system  trunk
bedload_seccurrent_gaia.f
Go to the documentation of this file.
1 ! **********************************
2  SUBROUTINE bedload_seccurrent_gaia
3 ! **********************************
4 !
5  &(ielmu,calfa,salfa)
6 !
7 !***********************************************************************
8 ! GAIA
9 !***********************************************************************
10 !
12 !
13 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
17 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
18 !
20  USE bief
22  IMPLICIT NONE
23 !
24 !
25 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
26 !
27  INTEGER, INTENT(IN) :: IELMU
28  TYPE(bief_obj), INTENT(INOUT) :: CALFA, SALFA
29 !
30 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
31 !
32  INTEGER I
33  DOUBLE PRECISION C, ALPHAL
34 !
35 
36 !
37 !FORMULATION FROM ENGELUND WITH GEOMETRIC RADIUS: TAN(TETA) = 7*(H/R)
38 !
39 
40 ! RADIUS OF CURVATURE MUST BE CALCULATED FROM FREE SURFACE SLOPE
41 ! USE OF ALPHA (0.75 FOR VERY ROUGH BOTTOMS, 1 FOR SMOOTH ONES)
42 
43 ! RK MODIFICATION FOR SECONDARY CURRENTS
44 !
45 ! COMPUTES THE GRADIENT OF THE FREE SURFACE IN X AND Y DIRECTION
46 !
47  CALL vector(t5,'=','GRADF X',ielmu,
48  & 1.d0,z,s,s,s,s,s,mesh,msk,maskel)
49  CALL vector(t6,'=','GRADF Y',ielmu,
50  & 1.d0,z,s,s,s,s,s,mesh,msk,maskel)
51 !
52 ! FOR PARALLEL COMPUTING
53 !
54  IF(ncsize.GT.1) THEN
55  CALL parcom (t5, 2, mesh)
56  CALL parcom (t6, 2, mesh)
57  ENDIF
58 !
59  CALL os('X=XY ',x=t5,y=unsv2d)
60  CALL os('X=XY ',x=t6,y=unsv2d)
61 !
62 ! COMPUTES THE X- AND Y-COMPONENTS OF THE SECONDARY CURRENT
63 ! ACCORDING TO ENGELUNG. TAU_X_SEC = C*QV, TAU_Y_SEC = C*QU
64 !
65 
66 ! BEWARE: THE VARIABLE ALPHA IS MORE THAN THE ALPHA FROM THE THEORY
67  alphal = alpha
68 !RK
69  IF(alpha.GT.0.d0) THEN
70  alphal = 7.d0 / alphal * xmve *grav
71  ELSE
72  alphal = 0.d0
73  ENDIF
74 ! CASE OF RADII FROM FILE
75  IF(havesecfile) THEN
76 ! ALPHA IS CALIBRATION FACTOR
77  alphal = 7.d0/alpha
78 ! CALCULATING
79  CALL os('X=-Y ' , t1 , radsec , radsec , c ) ! -RADIUS
80  CALL os('X=Y/Z ' , t8 , hn , t1 , c ) ! HN/-RADIUS
81 !RK CALL OS('X=CY ' , T8 , T8 , T8 , ALPHAL ) ! 7*HN/-RADIUS
82  CALL os('X=CY ' , t8 , t8 , t8 , alphal ) ! ALPHAL*HN/-RADIUS
83  CALL os('X=ATN(Y) ' , t8 , t8 , t8 , c ) ! ARCTAN(7*HN/-RADIUS)
84  ELSE
85 ! CALCULATING THE RADII FROM WATER LEVEL SLOPE
86 !
87  CALL os( 'X=YZ ' , t1 , t6 , u2d , c ) ! DZSDY*U2D
88  CALL os( 'X=YZ ' , t2 , t5 , v2d , c ) ! DZSDX*V2D
89  CALL os( 'X=X-Y ',x=t1,y=t2) ! U2D*DZSDY - V2D*DZSDX
90 !
91  CALL os( 'X=YZ ' , t2 , u2d , u2d , c ) ! U2D**2
92  CALL os( 'X=YZ ' , t3 , v2d , v2d , c ) ! V2D**2
93  CALL os( 'X=X+Y ' , t2 , t3 , t3 , c ) ! U2D**2+V2D**2
94 !
95  CALL os('X=Y/Z ' , t1 , t1 , t2, c ,2 , 0.d0,1.d-12) !(U2D*DZSDY - V2D*DZSDX)/(U2D**2+V2D**2)
96 !
97  CALL os( 'X=CX ' , t1 , t2 , t3 , alphal ) ! T1 * 7/ALPHAL*XMVE*GRAV
98  CALL os( 'X=XY ' , t1 , hn , t3 , c ) ! T1*HN
99 !
100 ! FOR ALL ROUGHNESS LAWS
101 !
102  CALL os( 'X=CXY ' , x=t1 , y=cf , c=0.5d0 )!T1*CF/2
103 !
104 ! TAU_X_SEK = -C*V2D : T5
105 ! TAU_Y_SEK = C*U2D : T6
106 !
107  CALL os('X=YZ ' , t5 , t1 , v2d, c ) ! C*V2D
108  CALL os('X=YZ ' , t6 , t1 , u2d, c ) ! C*U2D
109  CALL os('X=-Y ' , t6 , t6 , qv, c ) ! -C*U2D
110 !
111 ! SQRT(TAU_X_SEK**2+TAU_Y_SEK**2) : T3
112 !
113  CALL os('X=YZ ' , t2 , t5 , t5, c ) ! T2 = (C*V2D)**2
114  CALL os('X=YZ ' , t3 , t6 , t6, c ) ! T3 = (C*U2D)**2
115  CALL os('X=X+Y ' , t2 , t3 , t3, c ) ! T2 = (C*V2D)**2+(C*U2D)**2
116  CALL os('X=SQR(Y)' , x=t3 , y=t2 ) ! T3 = SQRT((C*U2D)**2+(C*V2D)**2
117 !
118 ! TAU_X_GES = TOB*EFFPNT*CALFA + TAU_X_SEK : T1
119 ! TAU_Y_GES = TOB*EFFPNT*SALFA + TAU_Y_SEK : T2
120 !
121  ENDIF
122 
123  CALL os( 'X=YZ ' , t1 , tob , coefpn , c ) ! TOB*EFFPNT
124  CALL os( 'X=YZ ' , t2 , t1 , salfa , c ) ! TOB*EFFPNT*SALFA
125  CALL os( 'X=YZ ' , t1 , t1 , calfa , c ) ! TOB*EFFPNT*CALFA
126 !
127  IF(havesecfile) THEN
128 !CALCULATION TAU_X_COR: TAU_X_COR=TAU_X*COS(TETA)-Y*SIN(TETA)
129  CALL os( 'X=COS(Y)' , t3 , t8 , t8 , c ) !COS(TETA)
130  CALL os( 'X=SIN(Y)' , t4 , t8 , t8 , c ) !SIN(TETA)
131  CALL os( 'X=YZ ' , t5 , t3 , t1 , c ) !TAU_X*COS(TETA)
132  CALL os( 'X=YZ ' , t6 , t4 , t2 , c ) !TAU_Y*SIN(TETA)
133  CALL os( 'X=Y-Z ' , t5 , t5 , t6 , c ) ! TAU_X_COR=TAU_X*COS(TETA)-TAU_Y*SIN(TETA)
134 !CALCULATION TAU_Y_COR: TAU_Y_COR=TAU_Y*SIN(TETA)+TAU_Y*COS(TETA)
135  CALL os( 'X=YZ ' , t6 , t4 , t1 , c ) !TAU_X*SIN(TETA)
136  CALL os( 'X=YZ ' , t7 , t3 , t2 , c ) !TAU_Y*COS(TETA)
137  CALL os( 'X=X+Y ' , t6 , t7 , t7 , c ) ! TAU_Y_COR=TAU_Y*COS(TETA)+TAU_Y*SIN(TETA)
138 !RK
139  CALL os( 'X=Y ' , t1 , t5 , t1 , c )
140  CALL os( 'X=Y ' , t2 , t6 , t1 , c )
141  ELSE
142  CALL os('X=X+Y ' , t1 , t5 , t3, c ) ! TAU_X_GES = TOB*CALFA+TAU_X_SEK
143  CALL os('X=X+Y ' , t2 , t6 , t3, c ) ! TAU_Y_GES = TOB*SALFA+TAU_Y_SEK
144 !
145  ENDIF
146 !
147 ! TAU_GES=SQRT(TAU_X_GES**2+TAU_Y_GES**2)
148 !
149  CALL os( 'X=YZ ' , t3 , t1 , t1 , c ) ! TAU_X_GES**2
150  CALL os( 'X=YZ ' , t4 , t2 , t2 , c ) ! TAU_Y_GES**2
151  CALL os('X=X+Y ' , t4 , t3 , t3, c ) ! TAU_X_GES**2+TAU_Y_GES**2
152  CALL os('X=SQR(Y)' , t4 , t4 , t3, c ) ! SQRT(TAU_X_GES**2+TAU_Y_GES**2)
153 !
154 ! NEW ANGLE
155 ! CALFA_NEW = COS(TAU_X_GES/TAU_GES)
156 ! SALFA_NEW = SIN(TAU_Y_GES/TAU_GES)
157 !
158  CALL os('X=Y/Z ' , t1 , t1 , t4, c ,2 , 0.d0,1.d-12) !TAU_X_GES/TAU_GES
159  CALL os('X=Y/Z ' , t2 , t2 , t4, c ,2 , 0.d0,1.d-12) !TAU_Y_GES/TAU_GES
160 !
161 ! TAKEN FROM EFFPNT UEBER
162 ! TO MAKE SURE THAT TAU_X_GES/TAU_GES IS IN RANGE [-1,1]
163 !
164  DO i=1,npoin
165  IF(t1%R(i).LT.-1.d0.OR.t1%R(i).GT.1.d0.OR.
166  & t2%R(i).LT.-1.d0.OR.t2%R(i).GT.1.d0) THEN
167  WRITE(lu,*) 'NOT ACCEPTABLE BORDER CROSSING',i
168  ENDIF
169  t1%R(i) = min(t1%R(i),1.d0)
170  t1%R(i) = max(t1%R(i),-1.d0)
171  t2%R(i) = min(t2%R(i),1.d0)
172  t2%R(i) = max(t2%R(i),-1.d0)
173  ENDDO
174 !
175 ! COEFPN_NEW = TAU_GES / TOB
176 !
177  CALL os('X=Y/Z ' , coefpn , t4 , tob, c ,2 , 0.d0,1.d-12) !COEFPN=TAU_GES/TOB
178  CALL os( 'X=Y ' ,x=calfa ,y=t1 ) ! (TAU_X_GES/TAU_GES)
179  CALL os( 'X=Y ' ,x=salfa ,y=t2 ) ! (TAU_Y_GES/TAU_GES)
180 !
181 ! FROM EFFPNT
182 ! FOR BOUNDARY NODES WITH IMPOSED FLOW :
183 ! QS IS NOT MODIFIED WHEN USER-DEFINED
184 !
185  DO i = 1 , nptfr
186  IF(liqbor%I(i).EQ.5) THEN
187  coefpn%R(mesh%NBOR%I(i)) = 1.d0
188  ENDIF
189  ENDDO
190 !
191 !-----------------------------------------------------------------------
192 !
193  RETURN
194  END
integer, pointer nptfr
Number of boundary points.
type(bief_obj), target s
Void structure.
type(bief_obj), target liqbor
Type of boundary conditions on sand transport rate.
type(bief_obj), target maskel
Mask.
type(bief_obj), target u2d
Components of depth-averaged velocity.
type(bief_obj), target tob
Bed shear stress [n/m2].
type(bief_obj), pointer t7
type(bief_obj), pointer t5
type(bief_obj), target v2d
double precision, dimension(:), pointer x
2d coordinates of the mesh
double precision xmve
Water density (from steering file of T2D or T3D)
double precision, dimension(:), pointer y
double precision, target alpha
Secondary Current Alpha Coefficient.
logical havesecfile
Secondary currents radii file.
type(bief_obj), target coefpn
Correction of transport for sloping bed effect.
type(bief_obj), pointer t8
type(bief_obj), target radsec
Curve radius for secondary currents.
type(bief_obj), target hn
Water depth.
type(bief_obj), target unsv2d
Inverse of integral of bases.
type(bief_obj), pointer t6
double precision grav
Gravity acceleration.
type(bief_obj), pointer t2
type(bief_obj), target z
Free surface elevation.
type(bief_obj), target qv
Y component of the flow rate.
logical msk
Include masking.
type(bief_obj), pointer t1
Aliases for work vectors in tb.
subroutine vector(VEC, OP, FORMUL, IELM1, XMUL, F, G, H, U, V, W, MESH, MSK, MASKEL, LEGO, ASSPAR)
Definition: vector.f:7
type(bief_obj), pointer t4
type(bief_obj), target cf
Quadratic friction coefficient.
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
type(bief_mesh), target mesh
Mesh structure.
subroutine parcom(X, ICOM, MESH)
Definition: parcom.f:7
subroutine bedload_seccurrent_gaia(IELMU, CALFA, SALFA)
integer, pointer npoin
Number of 2d points in the mesh.
type(bief_obj), pointer t3
Definition: bief.f:3