The TELEMAC-MASCARET system  trunk
calcs3d_micropol.f
Go to the documentation of this file.
1 ! ***************************
2  SUBROUTINE calcs3d_micropol
3 ! **************************
4  & (npoin2,tn,texp,timp,zprop,cf,un,vn,
5  & t2_1,t2_2,t2_3,t3_1,t3_2,debug)
6 !
7 !
8 !***********************************************************************
9 ! TELEMAC2D V7P2 21/05/2016
10 !***********************************************************************
11 !
12 !brief COMPUTES SOURCE TERMS FOR MICROPOL WAQ PROCESS IN 3D
13 ! WAQ PROCESS OF CODE_TRACER (MASCARET SYSTEM)
14 !
15 !history R. ATA
16 !+ 21/05/2016
17 !+ V7P0
18 !+ REAL CREATION
19 !
20 !history S.E. BOURBAN (HRW)
21 !+ 07/06/2017
22 !+ V7P3
23 !+ Indexing tracer (IND_*) to avoid conflicting naming convention
24 !+ between user defined tracers, water quality processes and
25 !+ ice processes. Introduction of the array RANK_*.
26 !
27 !history S.E. BOURBAN (HRW)
28 !+ 25/09/2017
29 !+ V7P3
30 !+ TEXP and TIMP are now additive to account for a variety of
31 !+ of sources / sinks on a given TRACER
32 !
33 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
34 !| CCSEDIM |-->| CONSTANT OF EXPONENTIAL DESINTEGRATION
35 !| CDISTRIB |-->| COEFFICIENT OF DISTRIBUTION (KD)
36 !| DEBUG |-->| IF NE.0 THEN DEBUG MODE
37 !| DT |-->| TIME STEP
38 !| ERO |-->| EROSION RATE
39 !| KDESORP |-->| KINETIC CONSTANT OF DESORPTION
40 !| NPOIN |-->| TOTAL NUMBER OF MESH NODES
41 !| TAUB |-->| BED SHEAR
42 !| TAUS |-->| CRITICAL STRESS OF RESUSPENSION
43 !| TAUR |-->| SEDIMENTATION CRITICAL STRESS
44 !| TEXP |<--| EXPLICIT SOURCE TERMS OF TRACERS
45 !| TIMP |<--| IMPLICIT SOURCE TERMS OF TRACERS
46 !| TN |-->| TRACERS
47 !| VITCHU |-->| SEDIMENT SETTLING VELOCITY
48 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
49 !
50  USE bief
54  & ro0,kdesorp,ccsedim,
55  & ind_ss,ind_sf,ind_c,ind_css,ind_csf
56  USE interface_waqtel, ex_calcs3d_micropol => calcs3d_micropol
57 !
58  IMPLICIT NONE
59 !
60 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
61 !
62  INTEGER , INTENT(IN ) :: npoin2
63  INTEGER , INTENT(IN ) :: debug
64  TYPE(bief_obj) , INTENT(IN ) :: tn,zprop,cf,un,vn
65  TYPE(bief_obj) , INTENT(INOUT) :: texp,timp
66  TYPE(bief_obj) , INTENT(INOUT) :: t2_1,t2_2,t2_3,t3_1,t3_2
67 !
68 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
69 !
70 ! LOCAL VARIABLES
71  INTEGER :: i
72  DOUBLE PRECISION, PARAMETER :: eps=1.d-3
73  DOUBLE PRECISION :: cc
74 !
75 !-----------------------------------------------------------------------
76 !
77  IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 0'
78 !
79 ! =======================================
80 ! PRELIMINARY COMPUTATIONS
81 ! =======================================
82 !
83  IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 1'
84 !
85 ! BED SHEAR STRESS (TAUB-STOCKED IN T2_1==>2D TABLE)
86 !
87  CALL taub_waqtel(cf,ro0,t2_1,npoin2,un,vn)
88 !
89 ! DEPOTION PROBABILITY (SED): STOCKED IN T2_2==>2D TABLE
90 !
91  CALL depos_fx(t2_2,t2_1,tn%ADR(ind_ss)%P,taus,vitchu,npoin2)
92 !
93 ! EROSION FLUX (RS): STOCKED IN T2_3 ==> 2D TABLE
94 !
95  CALL erosion_fx(t2_3,t2_1,tn%ADR(ind_sf)%P,taur,ero,1.d-10,
96  & npoin2)
97 !
98  IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 4'
99 !
100 !
101 ! =======================================
102 ! LET'S NOW COMPUTE SOURCE TERMS
103 ! =======================================
104 !
105 ! FIRST TRACER: SUSPENDED LOAD [SS] (IND_SS)
106 !
107 ! BED SOURCES
108 ! DO I=1,NPOIN2
109 ! TEXP%ADR(IND_SS)%P%R(I)=T2_3%R(I)-T2_2%R(I)
110 ! ENDDO
111 ! CALL OVD('X=Y/Z ',TEXP%ADR(IND_SS)%P%R,TEXP%ADR(IND_SS)%P%R,
112 ! & ZPROP%R,0.D0,NPOIN2,2,0.D0,EPS )
113  CALL os ('X=Y-Z ', x=t2_1,y=t2_3,z=t2_2 )
114  CALL ovd('X=X+CY/Z',texp%ADR(ind_ss)%P%R,t2_1%R,
115  & zprop%R,1.d0,npoin2,2,0.d0,eps )
116 !
117  IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 5'
118 !
119 ! SECOND TRACER: BED SEDIMENT [SF] (IND_SF)
120 ! warning: no advection neither diffusion for this tracer
121 !
122  DO i=1,npoin2
123 ! TEXP%ADR(IND_SF)%P%R(I)=T2_2%R(I)-T2_3%R(I)
124  texp%ADR(ind_sf)%P%R(i) = texp%ADR(ind_sf)%P%R(i) +
125  & t2_2%R(i) - t2_3%R(i)
126  ENDDO
127 !
128  IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 6'
129 !
130 ! THIRD TRACER: POLLUTANT DENSITY [C] (IND_C)
131 !
132 ! implicit part
133  CALL os( 'X=X+C ' ,x=timp%ADR(ind_c)%P,c=ccsedim )
134 ! explicit part
135 ! CALL OS( 'X=CY ' ,X=TEXP%ADR(IND_C)%P,Y=TN%ADR(IND_CSS)%P,
136 ! & C=KDESORP )
137  CALL os( 'X=CY ' ,x=t3_2,y=tn%ADR(ind_css)%P ,c=kdesorp )
138  CALL os( 'X=X+Y ' ,x=texp%ADR(ind_c)%P,y=t3_2 )
139  cc =-kdesorp*cdistrib
140 ! warning: the following term causes divergence of the code, it should
141 ! be traited implicitly- it is commented: to be investigated
142 ! more in depth
143 ! CALL OS( 'X=X+CYZ ' ,X=TEXP%ADR(IND_C)%P,Y=TN%ADR(IND_C)%P,
144 ! & Z=TN%ADR(IND_SS)%P ,C=CC )
145 !
146 !
147 ! FORTH TRACER: ABSORBED POLLUTANT BY SUSPENDED LOAD [CSS] (IND_CSS)
148 !
149 ! implicit part
150  CALL os( 'X=X+C ' ,x=timp%ADR(ind_css)%P,c=ccsedim )
151 ! explicit part
152 ! CALL OS( 'X=-Y ' ,X=TEXP%ADR(IND_CSS)%P,Y=TEXP%ADR(IND_C)%P )
153  CALL os( 'X=X-Y ' ,x=texp%ADR(ind_css)%P,y=t3_2 )
154  DO i=1,npoin2
155  t3_1%R(i) = t2_3%R(i)*tn%ADR(ind_csf)%P%R(i) -
156  & t2_2%R(i)*tn%ADR(ind_css)%P%R(i)
157  ENDDO
158  CALL ovd('X=Y/Z ' ,t3_1%R,t3_1%R,
159  & zprop%R,0.d0,npoin2,2,0.d0,eps )
160  CALL os( 'X=X+Y ' ,x=texp%ADR(ind_css)%P,y=t3_1 )
161 !
162  IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 8'
163 !
164 ! FIFTH TRACER: ABSORBED POLLUTANT BY BED SEDIMENT [CFF] (IND_CSF)
165 !
166  DO i=1,npoin2
167  texp%ADR(ind_csf)%P%R(i) = texp%ADR(ind_csf)%P%R(i) +
168  & t2_2%R(i)*tn%ADR(ind_css)%P%R(i) -
169  & t2_3%R(i)*tn%ADR(ind_csf)%P%R(i)
170  ENDDO
171  CALL os( 'X=X+CY ' ,x=texp%ADR(ind_csf)%P,y=tn%ADR(ind_csf)%P,
172  & c=-ccsedim )
173 !
174  IF(debug.GT.0)WRITE(lu,*)'IN MICROPOL3D, STEP 9'
175 !
176 ! MASS BALANCE: MASS ADDED BY EXPLICIT TERMS
177 ! (IMPLICIT PART IS ADDED IN CVDFTR)
178 !
179 !-----------------------------------------------------------------------
180 !
181  RETURN
182  END
subroutine taub_waqtel(CF, DENSITY, TAUB, NPOIN, UN, VN)
Definition: taub_waqtel.f:7
double precision vitchu
subroutine ovd(OP, X, Y, Z, C, NPOIN, IOPT, D, EPS)
Definition: ovd.f:7
double precision cdistrib
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
Definition: bief.f:3