source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/seddsr.F90 @ 10345

Last change on this file since 10345 was 10345, checked in by smasson, 2 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10344, see #2133

  • Property svn:keywords set to Id
File size: 33.7 KB
Line 
1MODULE 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$
27CONTAINS
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
142jflag1:     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
234jflag2:    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
329jflag3:     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
428jflag4:     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            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  &
594            &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn )
595            zsedtra(1) = zalpha + 0.25 * zsedtra(4) 
596            zsedtra(5) = zbeta  - zsedtra(4)
597            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4)
598            pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4)
599            ! H2S + O2
600            zalpha = zsedtra(1) - 2.0 * zsedtra(2)
601            zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2)
602            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2)
603            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  &
604            &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn )
605            zsedtra(1) = zalpha + 2.0 * zsedtra(2)
606            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2)
607            pwcp(ji,jk,jwso4) = zbeta - zsedtra(2)
608            ! NH4 + O2
609            zalpha = zsedtra(1) - 2.0 * zsedtra(3)
610            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3)
611            zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  &
612            &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn )
613            zsedtra(1) = zalpha + 2.0 * zsedtra(3)
614            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3)
615            ! FeS - O2
616            zalpha = zsedtra(1) - 2.0 * zsedtra(6)
617            zbeta  = zsedtra(4) + zsedtra(6)
618            zgamma = pwcp(ji,jk,jwso4) + zsedtra(6)
619            zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  &
620            &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn )
621            zsedtra(1) = zalpha + 2.0 * zsedtra(6)
622            zsedtra(4) = zbeta  - zsedtra(6)
623            pwcp(ji,jk,jwso4) = zgamma - zsedtra(6)
624!            ! Fe - H2S
625            zalpha = zsedtra(2) - zsedtra(4)
626            zbeta  = zsedtra(4) + zsedtra(6)
627            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
628            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  &
629            &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn )
630            zsedtra(2) = zalpha + zsedtra(4)
631            zsedtra(6) = zbeta  - zsedtra(4)
632            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4)
633            ! FEOH + H2S
634            zalpha = zsedtra(5) - 2.0 * zsedtra(2)
635            zbeta  = zsedtra(5) + zsedtra(4)
636            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
637            zdelta = pwcp(ji,jk,jwso4) + zsedtra(2)
638            zepsi  = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5)
639            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 )  &
640            &            + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) + rtrn )
641            zsedtra(5) = zalpha + 2.0 * zsedtra(2)
642            zsedtra(4) = zbeta  - zsedtra(5)
643            pwcp(ji,jk,jwso4) = zdelta - zsedtra(2)
644            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4)
645            pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5)
646            ! Fe - H2S
647            zalpha = zsedtra(2) - zsedtra(4)
648            zbeta  = zsedtra(4) + zsedtra(6)
649            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
650            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  &
651            &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn )
652            zsedtra(2) = zalpha + zsedtra(4)
653            zsedtra(6) = zbeta  - zsedtra(4)
654            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4)
655            ! FeS - O2
656            zalpha = zsedtra(1) - 2.0 * zsedtra(6)
657            zbeta  = zsedtra(4) + zsedtra(6)
658            zgamma = pwcp(ji,jk,jwso4) + zsedtra(6)
659            zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  &
660            &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn )
661            zsedtra(1) = zalpha + 2.0 * zsedtra(6)
662            zsedtra(4) = zbeta  - zsedtra(6)
663            pwcp(ji,jk,jwso4) = zgamma - zsedtra(6)
664            ! NH4 + O2
665            zalpha = zsedtra(1) - 2.0 * zsedtra(3)
666            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3)
667            zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  &
668            &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn )
669            zsedtra(1) = zalpha + 2.0 * zsedtra(3)
670            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3)
671            ! H2S + O2
672            zalpha = zsedtra(1) - 2.0 * zsedtra(2)
673            zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2)
674            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2)
675            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  &
676            &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn )
677            zsedtra(1) = zalpha + 2.0 * zsedtra(2)
678            pwcp(ji,jk,jwso4) = zbeta - zsedtra(2)
679            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2)
680            ! Fe + O2
681            zalpha = zsedtra(1) - 0.25 * zsedtra(4)
682            zbeta  = zsedtra(4) + zsedtra(5)
683            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
684            zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4)
685            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  &
686            &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn )
687            zsedtra(1) = zalpha + 0.25 * zsedtra(4)
688            zsedtra(5) = zbeta  - zsedtra(4)
689            pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4)
690            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4)
691            pwcp(ji,jk,jwoxy)  = zsedtra(1)
692            pwcp(ji,jk,jwh2s)  = zsedtra(2)
693            pwcp(ji,jk,jwnh4)  = zsedtra(3)
694            pwcp(ji,jk,jwfe2)  = zsedtra(4) + sedligand(ji,jk) + zsedfer
695            pwcp(ji,jk,jwno3)  = pwcp(ji,jk,jwno3)  + ( zsedtrn(3) - pwcp(ji,jk,jwnh4) )
696            solcp(ji,jk,jsfeo) = zsedtra(5) / zvolc(ji,jk,jsfeo)
697            solcp(ji,jk,jsfes) = zsedtra(6) / zvolc(ji,jk,jsfes)
698         END DO
699      END DO
700
701      IF( ln_timing )  CALL timing_stop('sed_dsr_redoxb')
702
703  END SUBROUTINE sed_dsr_redoxb
704
705END MODULE seddsr
Note: See TracBrowser for help on using the repository browser.