1 | MODULE seddsr |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE seddsr *** |
---|
4 | !! Sediment : dissolution and reaction in pore water related |
---|
5 | !! related to organic matter |
---|
6 | !!===================================================================== |
---|
7 | !! * Modules used |
---|
8 | USE sed ! sediment global variable |
---|
9 | USE sed_oce |
---|
10 | USE sedini |
---|
11 | USE lib_mpp ! distribued memory computing library |
---|
12 | USE lib_fortran |
---|
13 | |
---|
14 | IMPLICIT NONE |
---|
15 | PRIVATE |
---|
16 | |
---|
17 | PUBLIC sed_dsr |
---|
18 | |
---|
19 | !! * Module variables |
---|
20 | |
---|
21 | REAL(wp) :: zadsnh4 |
---|
22 | REAL(wp), DIMENSION(jpsol), PUBLIC :: dens_mol_wgt ! molecular density |
---|
23 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc ! temp. variables |
---|
24 | |
---|
25 | |
---|
26 | !! $Id$ |
---|
27 | CONTAINS |
---|
28 | |
---|
29 | SUBROUTINE sed_dsr( kt, knt ) |
---|
30 | !!---------------------------------------------------------------------- |
---|
31 | !! *** ROUTINE sed_dsr *** |
---|
32 | !! |
---|
33 | !! ** Purpose : computes pore water dissolution and reaction |
---|
34 | !! |
---|
35 | !! ** Methode : Computation of the redox reactions in sediment. |
---|
36 | !! The main redox reactions are solved in sed_dsr whereas |
---|
37 | !! the secondary reactions are solved in sed_dsr_redoxb. |
---|
38 | !! A strand spliting approach is being used here (see |
---|
39 | !! sed_dsr_redoxb for more information). |
---|
40 | !! |
---|
41 | !! History : |
---|
42 | !! ! 98-08 (E. Maier-Reimer, Christoph Heinze ) Original code |
---|
43 | !! ! 04-10 (N. Emprin, M. Gehlen ) f90 |
---|
44 | !! ! 06-04 (C. Ethe) Re-organization |
---|
45 | !! ! 19-08 (O. Aumont) Debugging and improvement of the model. |
---|
46 | !! The original method is replaced by a |
---|
47 | !! Strand splitting method which deals |
---|
48 | !! well with stiff reactions. |
---|
49 | !!---------------------------------------------------------------------- |
---|
50 | !! Arguments |
---|
51 | INTEGER, INTENT(in) :: kt, knt ! number of iteration |
---|
52 | ! --- local variables |
---|
53 | INTEGER :: ji, jk, js, jw, jn ! dummy looop indices |
---|
54 | |
---|
55 | REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2, zrearat3 ! reaction rate in pore water |
---|
56 | REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat ! undersaturation ; indice jpwatp1 is for calcite |
---|
57 | REAL(wp), DIMENSION(jpoce,jpksed) :: zkpoc, zkpos, zkpor, zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite |
---|
58 | REAL(wp), DIMENSION(jpoce) :: zsumtot |
---|
59 | REAL(wp) :: zsolid1, zsolid2, zsolid3, zvolw, zreasat |
---|
60 | REAL(wp) :: zsatur, zsatur2, znusil, zkpoca, zkpocb, zkpocc |
---|
61 | REAL(wp) :: zratio, zgamma, zbeta, zlimtmp, zundsat2 |
---|
62 | !! |
---|
63 | !!---------------------------------------------------------------------- |
---|
64 | |
---|
65 | IF( ln_timing ) CALL timing_start('sed_dsr') |
---|
66 | ! |
---|
67 | IF( kt == nitsed000 .AND. knt == 1 ) THEN |
---|
68 | IF (lwp) THEN |
---|
69 | WRITE(numsed,*) ' sed_dsr : Dissolution reaction ' |
---|
70 | WRITE(numsed,*) ' ' |
---|
71 | ENDIF |
---|
72 | ENDIF |
---|
73 | |
---|
74 | ! Initializations |
---|
75 | !---------------------- |
---|
76 | |
---|
77 | zrearat1(:,:) = 0. ; zundsat(:,:) = 0. ; zkpoc(:,:) = 0. |
---|
78 | zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. ; zrearat2(:,:) = 0. |
---|
79 | zlimso4(:,:) = 0. ; zkpor(:,:) = 0. ; zrearat3(:,:) = 0. |
---|
80 | zkpos (:,:) = 0. |
---|
81 | zsumtot(:) = rtrn |
---|
82 | |
---|
83 | ALLOCATE( zvolc(jpoce, jpksed, jpsol) ) |
---|
84 | zvolc(:,:,:) = 0. |
---|
85 | zadsnh4 = 1.0 / ( 1.0 + adsnh4 ) |
---|
86 | |
---|
87 | ! Inhibition terms for the different redox equations |
---|
88 | ! -------------------------------------------------- |
---|
89 | DO jk = 1, jpksed |
---|
90 | DO ji = 1, jpoce |
---|
91 | zkpoc(ji,jk) = reac_pocl |
---|
92 | zkpos(ji,jk) = reac_pocs |
---|
93 | zkpor(ji,jk) = reac_pocr |
---|
94 | END DO |
---|
95 | END DO |
---|
96 | |
---|
97 | ! Conversion of volume units |
---|
98 | !---------------------------- |
---|
99 | DO js = 1, jpsol |
---|
100 | DO jk = 1, jpksed |
---|
101 | DO ji = 1, jpoce |
---|
102 | zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) / & |
---|
103 | & ( volw3d(ji,jk) * 1.e-3 ) |
---|
104 | ENDDO |
---|
105 | ENDDO |
---|
106 | ENDDO |
---|
107 | |
---|
108 | !---------------------------------------------------------- |
---|
109 | ! 5. Beginning of solid reaction |
---|
110 | !--------------------------------------------------------- |
---|
111 | |
---|
112 | ! Definition of reaction rates [rearat]=sans dim |
---|
113 | ! For jk=1 no reaction (pure water without solid) for each solid compo |
---|
114 | zrearat1(:,:) = 0. |
---|
115 | zrearat2(:,:) = 0. |
---|
116 | zrearat3(:,:) = 0. |
---|
117 | |
---|
118 | zundsat(:,:) = pwcp(:,:,jwoxy) |
---|
119 | |
---|
120 | DO jk = 2, jpksed |
---|
121 | DO ji = 1, jpoce |
---|
122 | zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 ) |
---|
123 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
124 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
125 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
126 | zkpoca = zkpoc(ji,jk) * zlimo2(ji,jk) |
---|
127 | zkpocb = zkpos(ji,jk) * zlimo2(ji,jk) |
---|
128 | zkpocc = zkpor(ji,jk) * zlimo2(ji,jk) |
---|
129 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
130 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
131 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
132 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
133 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
134 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
135 | ENDDO |
---|
136 | ENDDO |
---|
137 | |
---|
138 | ! left hand side of coefficient matrix |
---|
139 | ! DO jn = 1, 5 |
---|
140 | DO jk = 2, jpksed |
---|
141 | DO ji = 1, jpoce |
---|
142 | jflag1: DO jn = 1, 10 |
---|
143 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
144 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
145 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
146 | zbeta = xksedo2 - pwcp(ji,jk,jwoxy) + so2ut * ( zrearat1(ji,jk) & |
---|
147 | & + zrearat2(ji,jk) + zrearat3(ji,jk) ) |
---|
148 | zgamma = - xksedo2 * pwcp(ji,jk,jwoxy) |
---|
149 | zundsat2 = zundsat(ji,jk) |
---|
150 | zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 |
---|
151 | zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 ) |
---|
152 | zkpoca = zkpoc(ji,jk) * zlimo2(ji,jk) |
---|
153 | zkpocb = zkpos(ji,jk) * zlimo2(ji,jk) |
---|
154 | zkpocc = zkpor(ji,jk) * zlimo2(ji,jk) |
---|
155 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
156 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
157 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
158 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
159 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
160 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
161 | IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN |
---|
162 | EXIT jflag1 |
---|
163 | ENDIF |
---|
164 | END DO jflag1 |
---|
165 | END DO |
---|
166 | END DO |
---|
167 | |
---|
168 | ! New solid concentration values (jk=2 to jksed) for each couple |
---|
169 | DO jk = 2, jpksed |
---|
170 | DO ji = 1, jpoce |
---|
171 | zreasat = zrearat1(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) |
---|
172 | solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat |
---|
173 | zreasat = zrearat2(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) |
---|
174 | solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat |
---|
175 | zreasat = zrearat3(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) |
---|
176 | solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat |
---|
177 | ENDDO |
---|
178 | ENDDO |
---|
179 | |
---|
180 | ! New pore water concentrations |
---|
181 | DO jk = 2, jpksed |
---|
182 | DO ji = 1, jpoce |
---|
183 | ! Acid Silicic |
---|
184 | pwcp(ji,jk,jwoxy) = zundsat(ji,jk) |
---|
185 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimo2(ji,jk) * zundsat(ji,jk) ! oxygen |
---|
186 | ! For DIC |
---|
187 | pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat |
---|
188 | zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 |
---|
189 | ! For Phosphate (in mol/l) |
---|
190 | pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r |
---|
191 | ! For iron (in mol/l) |
---|
192 | pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat |
---|
193 | ! For alkalinity |
---|
194 | pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) |
---|
195 | ! Ammonium |
---|
196 | pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 |
---|
197 | ! Ligands |
---|
198 | sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat - reac_ligc * sedligand(ji,jk) |
---|
199 | ENDDO |
---|
200 | ENDDO |
---|
201 | |
---|
202 | !-------------------------------------------------------------------- |
---|
203 | ! Begining POC denitrification and NO3- diffusion |
---|
204 | ! (indice n°5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) |
---|
205 | !-------------------------------------------------------------------- |
---|
206 | |
---|
207 | zrearat1(:,:) = 0. |
---|
208 | zrearat2(:,:) = 0. |
---|
209 | zrearat3(:,:) = 0. |
---|
210 | |
---|
211 | zundsat(:,:) = pwcp(:,:,jwno3) |
---|
212 | |
---|
213 | DO jk = 2, jpksed |
---|
214 | DO ji = 1, jpoce |
---|
215 | zlimno3(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) / ( zundsat(ji,jk) + xksedno3 ) |
---|
216 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
217 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
218 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
219 | zkpoca = zkpoc(ji,jk) * zlimno3(ji,jk) |
---|
220 | zkpocb = zkpos(ji,jk) * zlimno3(ji,jk) |
---|
221 | zkpocc = zkpor(ji,jk) * zlimno3(ji,jk) |
---|
222 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
223 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
224 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
225 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
226 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
227 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
228 | END DO |
---|
229 | END DO |
---|
230 | |
---|
231 | ! DO jn = 1, 5 |
---|
232 | DO jk = 2, jpksed |
---|
233 | DO ji = 1, jpoce |
---|
234 | jflag2: DO jn = 1, 10 |
---|
235 | zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) |
---|
236 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
237 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
238 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
239 | zbeta = xksedno3 - pwcp(ji,jk,jwno3) + srDnit * ( zrearat1(ji,jk) & |
---|
240 | & + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimtmp |
---|
241 | zgamma = - xksedno3 * pwcp(ji,jk,jwno3) |
---|
242 | zundsat2 = zundsat(ji,jk) |
---|
243 | zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 |
---|
244 | zlimno3(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) / ( zundsat(ji,jk) + xksedno3 ) |
---|
245 | zkpoca = zkpoc(ji,jk) * zlimno3(ji,jk) |
---|
246 | zkpocb = zkpos(ji,jk) * zlimno3(ji,jk) |
---|
247 | zkpocc = zkpor(ji,jk) * zlimno3(ji,jk) |
---|
248 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
249 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
250 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
251 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
252 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
253 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
254 | IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN |
---|
255 | EXIT jflag2 |
---|
256 | ENDIF |
---|
257 | END DO jflag2 |
---|
258 | END DO |
---|
259 | END DO |
---|
260 | |
---|
261 | |
---|
262 | ! New solid concentration values (jk=2 to jksed) for each couple |
---|
263 | DO jk = 2, jpksed |
---|
264 | DO ji = 1, jpoce |
---|
265 | zreasat = zrearat1(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) |
---|
266 | solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat |
---|
267 | zreasat = zrearat2(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) |
---|
268 | solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat |
---|
269 | zreasat = zrearat3(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) |
---|
270 | solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat |
---|
271 | ENDDO |
---|
272 | ENDDO |
---|
273 | |
---|
274 | ! New dissolved concentrations |
---|
275 | DO jk = 2, jpksed |
---|
276 | DO ji = 1, jpoce |
---|
277 | ! For nitrates |
---|
278 | pwcp(ji,jk,jwno3) = zundsat(ji,jk) |
---|
279 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimno3(ji,jk) * zundsat(ji,jk) |
---|
280 | ! For DIC |
---|
281 | pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat |
---|
282 | zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 |
---|
283 | ! For Phosphate (in mol/l) |
---|
284 | pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r |
---|
285 | ! Ligands |
---|
286 | sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat |
---|
287 | ! For iron (in mol/l) |
---|
288 | pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat |
---|
289 | ! For alkalinity |
---|
290 | pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srDnit + srno3 * zadsnh4 - 2.* spo4r ) |
---|
291 | ! Ammonium |
---|
292 | pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 |
---|
293 | ENDDO |
---|
294 | ENDDO |
---|
295 | |
---|
296 | !-------------------------------------------------------------------- |
---|
297 | ! Begining POC iron reduction |
---|
298 | ! (indice n�5 for couple POFe(OH)3 ie solcp(:,:,jspoc)/pwcp(:,:,jsfeo)) |
---|
299 | !-------------------------------------------------------------------- |
---|
300 | |
---|
301 | zrearat1(:,:) = 0. |
---|
302 | zrearat2(:,:) = 0. |
---|
303 | zrearat3(:,:) = 0. |
---|
304 | |
---|
305 | zundsat(:,:) = solcp(:,:,jsfeo) |
---|
306 | |
---|
307 | DO jk = 2, jpksed |
---|
308 | DO ji = 1, jpoce |
---|
309 | zlimfeo(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
310 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) / ( zundsat(ji,jk) + xksedfeo ) |
---|
311 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
312 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
313 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
314 | zkpoca = zkpoc(ji,jk) * zlimfeo(ji,jk) |
---|
315 | zkpocb = zkpos(ji,jk) * zlimfeo(ji,jk) |
---|
316 | zkpocc = zkpor(ji,jk) * zlimfeo(ji,jk) |
---|
317 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
318 | & ( 1. + zkpoca * zundsat(ji,jk) * dtsed2 ) |
---|
319 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
320 | & ( 1. + zkpocb * zundsat(ji,jk) * dtsed2 ) |
---|
321 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
322 | & ( 1. + zkpocc * zundsat(ji,jk) * dtsed2 ) |
---|
323 | END DO |
---|
324 | END DO |
---|
325 | |
---|
326 | ! DO jn = 1, 5 |
---|
327 | DO jk = 2, jpksed |
---|
328 | DO ji = 1, jpoce |
---|
329 | jflag3: DO jn = 1, 10 |
---|
330 | zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
331 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) |
---|
332 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
333 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
334 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
335 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) / zvolc(ji,jk,jsfeo) |
---|
336 | zbeta = xksedfeo - solcp(ji,jk,jsfeo) + 4.0 * zreasat * zlimtmp |
---|
337 | zgamma = -xksedfeo * solcp(ji,jk,jsfeo) |
---|
338 | zundsat2 = zundsat(ji,jk) |
---|
339 | zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 |
---|
340 | zlimfeo(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
341 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) / ( zundsat(ji,jk) + xksedfeo ) |
---|
342 | zkpoca = zkpoc(ji,jk) * zlimfeo(ji,jk) |
---|
343 | zkpocb = zkpos(ji,jk) * zlimfeo(ji,jk) |
---|
344 | zkpocc = zkpor(ji,jk) * zlimfeo(ji,jk) |
---|
345 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
346 | & ( 1. + zkpoca * zundsat(ji,jk) * dtsed2 ) |
---|
347 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
348 | & ( 1. + zkpocb * zundsat(ji,jk) * dtsed2 ) |
---|
349 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
350 | & ( 1. + zkpocc * zundsat(ji,jk) * dtsed2 ) |
---|
351 | IF ( ABS( (zundsat(ji,jk)-zundsat2)/( MAX(0.,zundsat2)+rtrn)) < 1E-8 ) THEN |
---|
352 | EXIT jflag3 |
---|
353 | ENDIF |
---|
354 | END DO jflag3 |
---|
355 | END DO |
---|
356 | END DO |
---|
357 | |
---|
358 | |
---|
359 | ! New solid concentration values (jk=2 to jksed) for each couple |
---|
360 | DO jk = 2, jpksed |
---|
361 | DO ji = 1, jpoce |
---|
362 | zreasat = zrearat1(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) |
---|
363 | solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat |
---|
364 | zreasat = zrearat2(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) |
---|
365 | solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat |
---|
366 | zreasat = zrearat3(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) |
---|
367 | solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat |
---|
368 | END DO |
---|
369 | END DO |
---|
370 | |
---|
371 | ! New dissolved concentrations |
---|
372 | DO jk = 2, jpksed |
---|
373 | DO ji = 1, jpoce |
---|
374 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimfeo(ji,jk) * zundsat(ji,jk) |
---|
375 | ! For FEOH |
---|
376 | solcp(ji,jk,jsfeo) = zundsat(ji,jk) |
---|
377 | ! For DIC |
---|
378 | pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat |
---|
379 | zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 |
---|
380 | ! For Phosphate (in mol/l) |
---|
381 | pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * ( spo4r + 4.0 * redfep ) |
---|
382 | ! Ligands |
---|
383 | sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat |
---|
384 | ! For iron (in mol/l) |
---|
385 | pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat |
---|
386 | ! For alkalinity |
---|
387 | pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) + 8.0 * zreasat |
---|
388 | ! Ammonium |
---|
389 | pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 |
---|
390 | pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + zreasat * 4.0 |
---|
391 | ENDDO |
---|
392 | ENDDO |
---|
393 | |
---|
394 | !-------------------------------------------------------------------- |
---|
395 | ! Begining POC denitrification and NO3- diffusion |
---|
396 | ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) |
---|
397 | !-------------------------------------------------------------------- |
---|
398 | |
---|
399 | zrearat1(:,:) = 0. |
---|
400 | zrearat2(:,:) = 0. |
---|
401 | zrearat3(:,:) = 0. |
---|
402 | |
---|
403 | zundsat(:,:) = pwcp(:,:,jwso4) |
---|
404 | |
---|
405 | DO jk = 2, jpksed |
---|
406 | DO ji = 1, jpoce |
---|
407 | zlimso4(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
408 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo) & |
---|
409 | & / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) / ( zundsat(ji,jk) + xksedso4 ) |
---|
410 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
411 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
412 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
413 | zkpoca = zkpoc(ji,jk) * zlimso4(ji,jk) |
---|
414 | zkpocb = zkpos(ji,jk) * zlimso4(ji,jk) |
---|
415 | zkpocc = zkpor(ji,jk) * zlimso4(ji,jk) |
---|
416 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
417 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
418 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
419 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
420 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
421 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
422 | END DO |
---|
423 | END DO |
---|
424 | ! |
---|
425 | ! DO jn = 1, 5 |
---|
426 | DO jk = 2, jpksed |
---|
427 | DO ji = 1, jpoce |
---|
428 | jflag4: DO jn = 1, 10 |
---|
429 | zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
430 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo) & |
---|
431 | & / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) |
---|
432 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
433 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
434 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
435 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) |
---|
436 | zbeta = xksedso4 - pwcp(ji,jk,jwso4) + 0.5 * zreasat * zlimtmp |
---|
437 | zgamma = - xksedso4 * pwcp(ji,jk,jwso4) |
---|
438 | zundsat2 = zundsat(ji,jk) |
---|
439 | zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 |
---|
440 | zlimso4(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
441 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo) & |
---|
442 | & / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) / ( zundsat(ji,jk) + xksedso4 ) |
---|
443 | zkpoca = zkpoc(ji,jk) * zlimso4(ji,jk) |
---|
444 | zkpocb = zkpos(ji,jk) * zlimso4(ji,jk) |
---|
445 | zkpocc = zkpor(ji,jk) * zlimso4(ji,jk) |
---|
446 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
447 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
448 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
449 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
450 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
451 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
452 | IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN |
---|
453 | EXIT jflag4 |
---|
454 | ENDIF |
---|
455 | END DO jflag4 |
---|
456 | END DO |
---|
457 | END DO |
---|
458 | |
---|
459 | ! New solid concentration values (jk=2 to jksed) for each couple |
---|
460 | DO jk = 2, jpksed |
---|
461 | DO ji = 1, jpoce |
---|
462 | zreasat = zrearat1(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) |
---|
463 | solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat |
---|
464 | zreasat = zrearat2(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) |
---|
465 | solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat |
---|
466 | zreasat = zrearat3(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) |
---|
467 | solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat |
---|
468 | ENDDO |
---|
469 | ENDDO |
---|
470 | ! |
---|
471 | ! New dissolved concentrations |
---|
472 | DO jk = 2, jpksed |
---|
473 | DO ji = 1, jpoce |
---|
474 | ! For sulfur |
---|
475 | pwcp(ji,jk,jwh2s) = pwcp(ji,jk,jwh2s) - ( zundsat(ji,jk) - pwcp(ji,jk,jwso4) ) |
---|
476 | pwcp(ji,jk,jwso4) = zundsat(ji,jk) |
---|
477 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimso4(ji,jk) * zundsat(ji,jk) |
---|
478 | ! For DIC |
---|
479 | pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat |
---|
480 | zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 |
---|
481 | ! For Phosphate (in mol/l) |
---|
482 | pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r |
---|
483 | ! Ligands |
---|
484 | sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat |
---|
485 | ! For iron (in mol/l) |
---|
486 | pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat |
---|
487 | ! For alkalinity |
---|
488 | pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) + zreasat |
---|
489 | ! Ammonium |
---|
490 | pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 |
---|
491 | ENDDO |
---|
492 | ENDDO |
---|
493 | |
---|
494 | ! Oxydation of the reduced products. Here only ammonium and ODU is accounted for |
---|
495 | ! There are two options here: A simple time splitting scheme and a modified |
---|
496 | ! Patankar scheme |
---|
497 | ! ------------------------------------------------------------------------------ |
---|
498 | |
---|
499 | call sed_dsr_redoxb |
---|
500 | |
---|
501 | ! -------------------------------------------------------------- |
---|
502 | ! 4/ Computation of the bioturbation coefficient |
---|
503 | ! This parameterization is taken from Archer et al. (2002) |
---|
504 | ! -------------------------------------------------------------- |
---|
505 | |
---|
506 | DO ji = 1, jpoce |
---|
507 | db(ji,:) = dbiot * zsumtot(ji) * pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6) |
---|
508 | END DO |
---|
509 | |
---|
510 | ! ------------------------------------------------------ |
---|
511 | ! Vertical variations of the bioturbation coefficient |
---|
512 | ! ------------------------------------------------------ |
---|
513 | IF (ln_btbz) THEN |
---|
514 | DO ji = 1, jpoce |
---|
515 | db(ji,:) = db(ji,:) * exp( -(profsedw(:) / dbtbzsc)**2 ) / (365.0 * 86400.0) |
---|
516 | END DO |
---|
517 | ELSE |
---|
518 | DO jk = 1, jpksed |
---|
519 | IF (profsedw(jk) > dbtbzsc) THEN |
---|
520 | db(:,jk) = 0.0 |
---|
521 | ENDIF |
---|
522 | END DO |
---|
523 | ENDIF |
---|
524 | |
---|
525 | IF (ln_irrig) THEN |
---|
526 | DO jk = 1, jpksed |
---|
527 | DO ji = 1, jpoce |
---|
528 | irrig(ji,jk) = ( 7.63752 - 7.4465 * exp( -0.89603 * zsumtot(ji) ) ) * pwcp(ji,1,jwoxy) & |
---|
529 | & / (pwcp(ji,1,jwoxy) + 20.E-6) |
---|
530 | irrig(ji,jk) = irrig(ji,jk) * exp( -(profsedw(jk) / xirrzsc) ) |
---|
531 | END DO |
---|
532 | END DO |
---|
533 | ELSE |
---|
534 | irrig(:,:) = 0.0 |
---|
535 | ENDIF |
---|
536 | |
---|
537 | DEALLOCATE( zvolc ) |
---|
538 | |
---|
539 | IF( ln_timing ) CALL timing_stop('sed_dsr') |
---|
540 | ! |
---|
541 | END SUBROUTINE sed_dsr |
---|
542 | |
---|
543 | SUBROUTINE sed_dsr_redoxb |
---|
544 | !!---------------------------------------------------------------------- |
---|
545 | !! *** ROUTINE sed_dsr_redox *** |
---|
546 | !! |
---|
547 | !! ** Purpose : computes secondary redox reactions |
---|
548 | !! |
---|
549 | !! ** Methode : It uses Strand splitter algorithm proposed by |
---|
550 | !! Nguyen et al. (2009) and modified by Wang et al. (2018) |
---|
551 | !! Basically, each equation is solved analytically when |
---|
552 | !! feasible, otherwise numerically at t+1/2. Then |
---|
553 | !! the last equation is solved at t+1. The other equations |
---|
554 | !! are then solved at t+1 starting in the reverse order. |
---|
555 | !! Ideally, it's better to start from the fastest reaction |
---|
556 | !! to the slowest and then reverse the order to finish up |
---|
557 | !! with the fastest one. But random order works well also. |
---|
558 | !! The scheme is second order, positive and mass |
---|
559 | !! conserving. It works well for stiff systems. |
---|
560 | !! |
---|
561 | !! History : |
---|
562 | !! ! 18-08 (O. Aumont) Original code |
---|
563 | !!---------------------------------------------------------------------- |
---|
564 | !! Arguments |
---|
565 | ! --- local variables |
---|
566 | INTEGER :: ji, jk, jn ! dummy looop indices |
---|
567 | |
---|
568 | REAL, DIMENSION(6) :: zsedtrn, zsedtra |
---|
569 | REAL(wp) :: zalpha, zbeta, zgamma, zdelta, zepsi, zsedfer |
---|
570 | !! |
---|
571 | !!---------------------------------------------------------------------- |
---|
572 | |
---|
573 | IF( ln_timing ) CALL timing_start('sed_dsr_redoxb') |
---|
574 | |
---|
575 | DO ji = 1, jpoce |
---|
576 | DO jk = 2, jpksed |
---|
577 | zsedtrn(1) = pwcp(ji,jk,jwoxy) |
---|
578 | zsedtrn(2) = MAX(0., pwcp(ji,jk,jwh2s) ) |
---|
579 | zsedtrn(3) = pwcp(ji,jk,jwnh4) |
---|
580 | zsedtrn(4) = MAX(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) |
---|
581 | zsedfer = MIN(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) |
---|
582 | zsedtrn(5) = solcp(ji,jk,jsfeo) * zvolc(ji,jk,jsfeo) |
---|
583 | zsedtrn(6) = solcp(ji,jk,jsfes) * zvolc(ji,jk,jsfes) |
---|
584 | zsedtra(:) = zsedtrn(:) |
---|
585 | |
---|
586 | ! First pass of the scheme. At the end, it is 1st order |
---|
587 | ! ----------------------------------------------------- |
---|
588 | ! Fe + O2 |
---|
589 | zalpha = zsedtra(1) - 0.25 * zsedtra(4) |
---|
590 | zbeta = zsedtra(4) + zsedtra(5) |
---|
591 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) |
---|
592 | zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) |
---|
593 | IF ( zalpha == 0. ) THEN |
---|
594 | zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) |
---|
595 | ELSE |
---|
596 | zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
597 | & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) |
---|
598 | ENDIF |
---|
599 | zsedtra(1) = zalpha + 0.25 * zsedtra(4) |
---|
600 | zsedtra(5) = zbeta - zsedtra(4) |
---|
601 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) |
---|
602 | pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) |
---|
603 | ! H2S + O2 |
---|
604 | zalpha = zsedtra(1) - 2.0 * zsedtra(2) |
---|
605 | zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) |
---|
606 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) |
---|
607 | IF ( zalpha == 0. ) THEN |
---|
608 | zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) |
---|
609 | ELSE |
---|
610 | zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
611 | & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) |
---|
612 | ENDIF |
---|
613 | zsedtra(1) = zalpha + 2.0 * zsedtra(2) |
---|
614 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) |
---|
615 | pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) |
---|
616 | ! NH4 + O2 |
---|
617 | zalpha = zsedtra(1) - 2.0 * zsedtra(3) |
---|
618 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) |
---|
619 | IF ( zalpha == 0. ) THEN |
---|
620 | zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0 ) |
---|
621 | ELSE |
---|
622 | zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
623 | & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) |
---|
624 | ENDIF |
---|
625 | zsedtra(1) = zalpha + 2.0 * zsedtra(3) |
---|
626 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) |
---|
627 | ! FeS - O2 |
---|
628 | zalpha = zsedtra(1) - 2.0 * zsedtra(6) |
---|
629 | zbeta = zsedtra(4) + zsedtra(6) |
---|
630 | zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) |
---|
631 | IF ( zalpha == 0. ) THEN |
---|
632 | zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2.0 ) |
---|
633 | ELSE |
---|
634 | zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
635 | & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) |
---|
636 | ENDIF |
---|
637 | zsedtra(1) = zalpha + 2.0 * zsedtra(6) |
---|
638 | zsedtra(4) = zbeta - zsedtra(6) |
---|
639 | pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) |
---|
640 | ! ! Fe - H2S |
---|
641 | zalpha = zsedtra(2) - zsedtra(4) |
---|
642 | zbeta = zsedtra(4) + zsedtra(6) |
---|
643 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) |
---|
644 | IF ( zalpha == 0. ) THEN |
---|
645 | zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) |
---|
646 | ELSE |
---|
647 | zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
648 | & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) |
---|
649 | ENDIF |
---|
650 | zsedtra(2) = zalpha + zsedtra(4) |
---|
651 | zsedtra(6) = zbeta - zsedtra(4) |
---|
652 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) |
---|
653 | ! FEOH + H2S |
---|
654 | zalpha = zsedtra(5) - 2.0 * zsedtra(2) |
---|
655 | zbeta = zsedtra(5) + zsedtra(4) |
---|
656 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) |
---|
657 | zdelta = pwcp(ji,jk,jwso4) + zsedtra(2) |
---|
658 | zepsi = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) |
---|
659 | IF ( zalpha == 0. ) THEN |
---|
660 | zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_feh2s * dtsed2 ) |
---|
661 | ELSE |
---|
662 | zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 ) & |
---|
663 | & + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) ) |
---|
664 | ENDIF |
---|
665 | zsedtra(5) = zalpha + 2.0 * zsedtra(2) |
---|
666 | zsedtra(4) = zbeta - zsedtra(5) |
---|
667 | pwcp(ji,jk,jwso4) = zdelta - zsedtra(2) |
---|
668 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) |
---|
669 | pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5) |
---|
670 | ! Fe - H2S |
---|
671 | zalpha = zsedtra(2) - zsedtra(4) |
---|
672 | zbeta = zsedtra(4) + zsedtra(6) |
---|
673 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) |
---|
674 | IF ( zalpha == 0. ) THEN |
---|
675 | zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) |
---|
676 | ELSE |
---|
677 | zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
678 | & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) |
---|
679 | ENDIF |
---|
680 | zsedtra(2) = zalpha + zsedtra(4) |
---|
681 | zsedtra(6) = zbeta - zsedtra(4) |
---|
682 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) |
---|
683 | ! FeS - O2 |
---|
684 | zalpha = zsedtra(1) - 2.0 * zsedtra(6) |
---|
685 | zbeta = zsedtra(4) + zsedtra(6) |
---|
686 | zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) |
---|
687 | IF (zalpha == 0.) THEN |
---|
688 | zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2. ) |
---|
689 | ELSE |
---|
690 | zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
691 | & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) |
---|
692 | ENDIF |
---|
693 | zsedtra(1) = zalpha + 2.0 * zsedtra(6) |
---|
694 | zsedtra(4) = zbeta - zsedtra(6) |
---|
695 | pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) |
---|
696 | ! NH4 + O2 |
---|
697 | zalpha = zsedtra(1) - 2.0 * zsedtra(3) |
---|
698 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) |
---|
699 | IF (zalpha == 0.) THEN |
---|
700 | zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0) |
---|
701 | ELSE |
---|
702 | zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
703 | & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) |
---|
704 | ENDIF |
---|
705 | zsedtra(1) = zalpha + 2.0 * zsedtra(3) |
---|
706 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) |
---|
707 | ! H2S + O2 |
---|
708 | zalpha = zsedtra(1) - 2.0 * zsedtra(2) |
---|
709 | zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) |
---|
710 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) |
---|
711 | IF ( zalpha == 0. ) THEN |
---|
712 | zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) |
---|
713 | ELSE |
---|
714 | zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
715 | & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) |
---|
716 | ENDIF |
---|
717 | zsedtra(1) = zalpha + 2.0 * zsedtra(2) |
---|
718 | pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) |
---|
719 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) |
---|
720 | ! Fe + O2 |
---|
721 | zalpha = zsedtra(1) - 0.25 * zsedtra(4) |
---|
722 | zbeta = zsedtra(4) + zsedtra(5) |
---|
723 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) |
---|
724 | zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) |
---|
725 | IF ( zalpha == 0. ) THEN |
---|
726 | zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) |
---|
727 | ELSE |
---|
728 | zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
729 | & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) |
---|
730 | ENDIF |
---|
731 | zsedtra(1) = zalpha + 0.25 * zsedtra(4) |
---|
732 | zsedtra(5) = zbeta - zsedtra(4) |
---|
733 | pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) |
---|
734 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) |
---|
735 | pwcp(ji,jk,jwoxy) = zsedtra(1) |
---|
736 | pwcp(ji,jk,jwh2s) = zsedtra(2) |
---|
737 | pwcp(ji,jk,jwnh4) = zsedtra(3) |
---|
738 | pwcp(ji,jk,jwfe2) = zsedtra(4) + sedligand(ji,jk) + zsedfer |
---|
739 | pwcp(ji,jk,jwno3) = pwcp(ji,jk,jwno3) + ( zsedtrn(3) - pwcp(ji,jk,jwnh4) ) |
---|
740 | solcp(ji,jk,jsfeo) = zsedtra(5) / zvolc(ji,jk,jsfeo) |
---|
741 | solcp(ji,jk,jsfes) = zsedtra(6) / zvolc(ji,jk,jsfes) |
---|
742 | END DO |
---|
743 | END DO |
---|
744 | |
---|
745 | IF( ln_timing ) CALL timing_stop('sed_dsr_redoxb') |
---|
746 | |
---|
747 | END SUBROUTINE sed_dsr_redoxb |
---|
748 | |
---|
749 | END MODULE seddsr |
---|