The TELEMAC-MASCARET system  trunk
activelayer.f
Go to the documentation of this file.
1 ! **********************************
2  SUBROUTINE activelayer
3 ! **********************************
4 !
5 !***********************************************************************
6 ! GAIA V8P1 05/02/2020
7 !***********************************************************************
8 !
11 !
15 
17 
18  IMPLICIT NONE
19 
20  INTEGER :: ICLA, IPOIN, I, K, DMIN_I(1)
21  DOUBLE PRECISION :: TAUC,DMAX,SUMME,DSTAR,FMAX,DENS,DMIN
22  DOUBLE PRECISION :: ACMIN
23 
24 !CALCULATION OF THE ACTIVE LAYER THICKNESS
25  SELECT CASE (alt_model)
26  CASE (0)
27  DO ipoin=1,npoin
28  elay%R(ipoin)= elay0
29  END DO
30 !-----------------------------------------------------------------------
31 ! HUNZIKER & GUENTHER
32 !-----------------------------------------------------------------------
33 !
34  CASE (1)
35 ! In the first time step the active layer is only virtual and the values in the
36 ! first layer are zero
37  k=1
38  IF(lt.LE.1) k=2
39  DO ipoin=1,npoin
40 !-----------------------------------------------------------------------
41 ! DMAX - NEW APPROXIMATION: DMAX = D99 / BEFORE DMAX = MAX(DCLA)
42 ! ALL CLASSES MUST BE NON-COHESIVE
43 !-----------------------------------------------------------------------
44  fmax = 0.99d0*mass_sand_tot(k,ipoin)
45  summe = mass_sand(1,k,ipoin)
46  dmax = 0.d0
47  IF (summe.GE.fmax) THEN
48  dmax = dcla(1)
49  ELSE
50  DO i=2,nsand
51  IF(mass_sand(i,k,ipoin).GT.0.d0) THEN
52  summe = mass_sand(i,k,ipoin) + summe
53  IF(summe.GE.fmax.AND.dmax.EQ.0.d0) THEN
54  dmax = dcla(i-1) + ( (dcla(i)-dcla(i-1)) /
55  & mass_sand(i,k,ipoin))*
56  & (fmax-(summe-mass_sand(i,k,ipoin)))
57  ENDIF
58  ENDIF
59  ENDDO
60  ENDIF
61 !DMAX = MAX(DCLA(1), DMAX)
62  IF (dmax.LE.0.d0.AND.summe.GT.zero) THEN
63  WRITE(lu,*)'UHM DMAXERROR'
64  DO i=1,nsicla
65  WRITE(lu,*)'DMAXERROR',ipoin,i,dcla(i),mass_sand(i,k,ipoin)
66  & ,summe
67  WRITE(lu,*)zf%R(ipoin)-zr%R(ipoin)
68  ENDDO
69  stop
70  ENDIF
71  IF(dmax.GT.maxval(dcla(1:nsicla)).AND.summe.GT.zero)THEN
72  print*,'DMAX',dmax,dcla(1:nsicla)
73  stop
74  ENDIF
75  elay%R(ipoin)= 5.d0 * dmax
76  END DO
77 !
78 !-----------------------------------------------------------------------
79 ! FREDSOE & DEIGAARD 1992
80 ! DENSITY IS TAKEN FROM THE FIRST CLASS
81 !-----------------------------------------------------------------------
82 !
83  CASE (2)
84  DO ipoin=1,npoin
85  elay%R(ipoin)= 2.d0 * tob%R(ipoin) / (grav*(xmvs0(1)-xmve))
86  & / tan(phised/180.0d0*pi) / (1.d0-xkv0(1))
87  END DO
88 
89 !
90 !-----------------------------------------------------------------------
91 ! VAN RIJN 1993
92 !-----------------------------------------------------------------------
93 !
94  CASE (3)
96 ! CRICITAL SHEAR STRESS AND DENSITY AND DIAMETER FOR THE FINEST CLASS
97  dmin = minval(dcla(1:nsicla))
98  dmin_i = minloc(dcla(1:nsicla))
99  acmin = ac(dmin_i(1))
100  tauc = acmin*((xmvs0(1)-xmve)*grav*dmin)
101 !
102  DO ipoin=1,npoin
103  elay%R(ipoin) = 0.d0
104  IF(tauc.GT.0.d0) THEN
105  IF(tob%R(ipoin).GT.tauc) THEN
106  dens = (xmvs0(1) - xmve) / xmve
107  dstar = acladm%R(ipoin)*(grav*dens/vce**2)**(1.d0/3.d0)
108  elay%R(ipoin) = 0.3d0*(dstar**0.7d0)*
109  & ((tob%R(ipoin)-tauc)/tauc)**0.5*acladm%R(ipoin)
110  ENDIF
111  ENDIF
112  ENDDO
113 !-----------------------------------------------------------------------
114 ! WONG 2006
115 ! DENSITY IS TAKEN FROM THE FIRST CLASS!
116 ! HARD CODED SHIELDS THRESHOLD -- THIS FORMULA IS NOT SUITABLE FOR
117 ! FINE SEDIMENTS
118 !-----------------------------------------------------------------------
119  CASE (4)
121  dmin = minval(dcla(1:nsicla))
122  DO ipoin=1,npoin
123  IF((tob%R(ipoin)/(xmvs0(1)-xmve)/grav/
124  & dmin).LT.0.0549d0) THEN
125  elay%R(ipoin) = 0.d0
126  ELSE
127  elay%R(ipoin)=5.0d0*acladm%R(ipoin)*
128  & ((tob%R(ipoin)/(xmvs0(1)-xmve)/
129  & grav/dmin) -0.0549d0)**0.56d0
130  ENDIF
131  ENDDO
132 !-----------------------------------------------------------------------
133 ! MALCHEREK 2003
134 !-----------------------------------------------------------------------
135  CASE(5)
136 ! CRICITAL SHEAR STRESS FOR THE FINEST SEDIMENT
137 ! DENSITY FROM THE FIRST CLASS
138  acmin = minval(ac(1:nsicla))
139  dmin = minval(dcla(1:nsicla))
140  tauc = acmin*((xmvs0(1)-xmve)*grav*dmin)
141 !-----------------------------------------------------------------------
142 ! D90 - ONLY FIRST APPROXIMATION!
143 ! ALL CLASSES MUST BE NON-COHESIVE
144 !-----------------------------------------------------------------------
145 ! In the first time step the active layer is only virtual and the values in the
146 ! first layer are zero
147  k=1
148  IF(lt.LE.1) k=2
149  DO ipoin=1,npoin
150  fmax = 0.90d0*mass_sand_tot(k,ipoin)
151  summe = mass_sand(1,k,ipoin)
152  d90 = 0.d0
153  IF (summe.GE.fmax) THEN
154  d90 = dcla(1)
155  ELSE
156  DO i=2,nsand
157  IF(mass_sand(i,k,ipoin).GT.0.d0) THEN
158  summe = mass_sand(i,k,ipoin) + summe
159  IF(summe.GE.fmax.AND.d90.EQ.0.d0) THEN
160  d90 = dcla(i-1) + ( (dcla(i)-dcla(i-1)) /
161  & mass_sand(i,k,ipoin)) *
162  & (fmax-(summe-mass_sand(i,k,ipoin)))
163  ENDIF
164  ENDIF
165  ENDDO
166  ENDIF
167 !D90 = MAX(DCLA(1), D90)
168  IF (d90.LE.0.d0.AND.summe.GT.zero) THEN
169  WRITE(lu,*)'UHM D90ERROR', d90, fmax, summe
170  DO i=1,nsicla
171  WRITE(lu,*)'D90ERROR',ipoin,i,dcla(i),
172  & mass_sand(i,k,ipoin),summe
173  ENDDO
174  WRITE(lu,*)'D90ERROR'
175  ENDIF
176  elay%R(ipoin) = 0.d0
177  IF(tauc.GT.0.d0)
178  & elay%R(ipoin) = d90 / (1.d0-xkv0(1)) *
179  & max(1.d0,tob%R(ipoin)/tauc)
180  ENDDO
181 !-----------------------------------------------------------------------
182 ! SISYPHE
183 !-----------------------------------------------------------------------
184  CASE(6)
186  DO ipoin=1,npoin
187  elay%R(ipoin)=3*acladm%R(ipoin)
188  ENDDO
189  END SELECT
190 
191  END SUBROUTINE
type(bief_obj), target acladm
Mean diameter of active-layer.
type(bief_obj), target tob
Bed shear stress [n/m2].
double precision pi
Pi.
double precision elay0
Wanted active layer thickness; ELAYO is a target value for ELAY, but ELAY may be lower than ELAY0 if ...
type(bief_obj), target zr
Non erodable (rigid) bottom elevation.
subroutine activelayer
Definition: activelayer.f:4
double precision, dimension(:,:), allocatable mass_sand_tot
Surface total mass of sand (kg/m2), for ilayer,ipoin.
double precision xmve
Water density (from steering file of T2D or T3D)
double precision, dimension(:,:,:), allocatable, target mass_sand
Surface mass of sand (kg/m2), for isand,ilayer,ipoin.
double precision vce
Water viscosity: it is defined here because the viscosity set in TELEMAC2D or TELEMAC3D may not b the...
subroutine mean_grain_size_gaia
double precision, dimension(nsiclm) xmvs0
Sand density.
type(bief_obj), target elay
Active layer thickness.
integer nsand
Total number of sand.
double precision zero
Parameter used for clipping variables or testing values against zero.
integer alt_model
CHOOSE A MODEL FOR ESTIMATION OF A DYNAMIC ACTIVE LAYER THICKNESS.
double precision, dimension(nsiclm), target dcla
Sediment diameter for each class It is only relevant for non-cohesive sediments. For the bedload...
double precision, target phised
Friction angle of the sediment.
integer, target nsicla
Number of sediment classes of bed material (less than NISCLM)
double precision grav
Gravity acceleration.
double precision, dimension(nlaymax), target xkv0
Initial porosity by layers.
integer, target lt
Numero du pas de temps.
double precision, dimension(nsiclm), target ac
Critical shields parameter.
double precision, target d90
Sediment diameter D90, for sand when only.
type(bief_obj), target zf
Bottom elevation.
integer, pointer npoin
Number of 2d points in the mesh.