New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
seddsr.F90 in NEMO/branches/UKMO/NEMO_4.0.1_mirror/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_mirror/src/TOP/PISCES/SED/seddsr.F90 @ 11715

Last change on this file since 11715 was 11715, checked in by davestorkey, 5 years ago

UKMO/NEMO_4.0.1_mirror : Remove SVN keywords.
File size: 35.4 KB
RevLine 
[3443]1MODULE seddsr
2   !!======================================================================
3   !!              ***  MODULE  seddsr  ***
[10222]4   !!    Sediment : dissolution and reaction in pore water related
5   !!    related to organic matter
[3443]6   !!=====================================================================
7   !! * Modules used
8   USE sed     ! sediment global variable
[10222]9   USE sed_oce
10   USE sedini
11   USE lib_mpp         ! distribued memory computing library
12   USE lib_fortran
[3443]13
[10222]14   IMPLICIT NONE
15   PRIVATE
16
[3443]17   PUBLIC sed_dsr
18
19   !! * Module variables
20
[10222]21   REAL(wp) :: zadsnh4
22   REAL(wp), DIMENSION(jpsol), PUBLIC      :: dens_mol_wgt  ! molecular density
23   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc    ! temp. variables
[3443]24
[10222]25
[5215]26   !! $Id$
[3443]27CONTAINS
28   
[10222]29   SUBROUTINE sed_dsr( kt, knt ) 
[3443]30      !!----------------------------------------------------------------------
31      !!                   ***  ROUTINE sed_dsr  ***
32      !!
33      !!  ** Purpose :  computes pore water dissolution and reaction
34      !!
[10222]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).
[3443]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
[10222]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.
[3443]49      !!----------------------------------------------------------------------
50      !! Arguments
[10222]51      INTEGER, INTENT(in) ::   kt, knt       ! number of iteration
[3443]52      ! --- local variables
[10222]53      INTEGER :: ji, jk, js, jw, jn   ! dummy looop indices
[3443]54
[10222]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
[3443]59      REAL(wp)  ::  zsolid1, zsolid2, zsolid3, zvolw, zreasat
[10222]60      REAL(wp)  ::  zsatur, zsatur2, znusil, zkpoca, zkpocb, zkpocc
61      REAL(wp)  ::  zratio, zgamma, zbeta, zlimtmp, zundsat2
[3443]62      !!
63      !!----------------------------------------------------------------------
64
[10222]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
[3443]72      ENDIF
73
[10222]74     ! Initializations
75     !----------------------
[3443]76     
[10222]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 )
[3443]86
[10222]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
[3443]96
97      ! Conversion of volume units
98      !----------------------------
99      DO js = 1, jpsol
100         DO jk = 1, jpksed
[10222]101            DO ji = 1, jpoce
[3443]102               zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) /  &
[10222]103                  &              ( volw3d(ji,jk) * 1.e-3 )
[3443]104            ENDDO
105         ENDDO
106      ENDDO
107
108      !----------------------------------------------------------
[10222]109      ! 5.  Beginning of solid reaction
[3443]110      !---------------------------------------------------------
111     
112      ! Definition of reaction rates [rearat]=sans dim
113      ! For jk=1 no reaction (pure water without solid) for each solid compo
[10222]114      zrearat1(:,:) = 0.
115      zrearat2(:,:) = 0.
116      zrearat3(:,:) = 0.
[3443]117
[10222]118      zundsat(:,:) = pwcp(:,:,jwoxy)
[3443]119
120      DO jk = 2, jpksed
121         DO ji = 1, jpoce
[10222]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 )
[3443]135         ENDDO
136      ENDDO
137
[10222]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
[3443]167
168      ! New solid concentration values (jk=2 to jksed) for each couple
[10222]169      DO jk = 2, jpksed
[3443]170         DO ji = 1, jpoce
[10222]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
[3443]177         ENDDO
178      ENDDO
179
180      ! New pore water concentrations   
[10222]181      DO jk = 2, jpksed
[3443]182         DO ji = 1, jpoce
183            ! Acid Silicic
[10222]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         
[3443]186            ! For DIC
187            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
[10222]188            zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3
[3443]189            ! For Phosphate (in mol/l)
[10222]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
[3443]193            ! For alkalinity
[10222]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)
[3443]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
[10222]207      zrearat1(:,:) = 0.
208      zrearat2(:,:) = 0.
209      zrearat3(:,:) = 0.
210
211      zundsat(:,:) = pwcp(:,:,jwno3)
212
213      DO jk = 2, jpksed
[3443]214         DO ji = 1, jpoce
[10222]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
[3443]232      DO jk = 2, jpksed
233         DO ji = 1, jpoce
[10222]234jflag2:    DO jn = 1, 10
235               zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) )
[3443]236               zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc)
[10222]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
[3443]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
[10222]265            zreasat = zrearat1(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc)
[3443]266            solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat
[10222]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
[3443]271         ENDDO
272      ENDDO
273
274      ! New dissolved concentrations
[10222]275      DO jk = 2, jpksed
[3443]276         DO ji = 1, jpoce
277            ! For nitrates
[10222]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)
[3443]280            ! For DIC
281            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
[10222]282            zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3
[3443]283            ! For Phosphate (in mol/l)
284            pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * spo4r           
[10222]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
[3443]289            ! For alkalinity
[10222]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
[3443]293         ENDDO
294      ENDDO
295
[10222]296      !--------------------------------------------------------------------
297      ! Begining POC iron reduction
298      ! (indice n�5 for couple POFe(OH)3 ie solcp(:,:,jspoc)/pwcp(:,:,jsfeo))
299      !--------------------------------------------------------------------
[3443]300
[10222]301      zrearat1(:,:) = 0.
302      zrearat2(:,:) = 0.
303      zrearat3(:,:) = 0.
[3443]304
[10222]305      zundsat(:,:) = solcp(:,:,jsfeo)
[3443]306
[10222]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
[3443]325
[10222]326!      DO jn = 1, 5
327      DO jk = 2, jpksed
[3443]328         DO ji = 1, jpoce
[10222]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
[3443]357
358
[10222]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
[3443]370
[10222]371      ! New dissolved concentrations
372      DO jk = 2, jpksed
[3443]373         DO ji = 1, jpoce
[10222]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
[3443]391         ENDDO
392      ENDDO
393
[10222]394      !--------------------------------------------------------------------
395      ! Begining POC denitrification and NO3- diffusion
396      ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3))
397      !--------------------------------------------------------------------
[3443]398
[10222]399      zrearat1(:,:) = 0.
400      zrearat2(:,:) = 0.
401      zrearat3(:,:) = 0.
[3443]402
[10222]403      zundsat(:,:) = pwcp(:,:,jwso4)
[3443]404
[10222]405      DO jk = 2, jpksed
[3443]406         DO ji = 1, jpoce
[10222]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
[3443]426      DO jk = 2, jpksed
427         DO ji = 1, jpoce
[10222]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
[3443]456         END DO
457      END DO
458
[10222]459     ! New solid concentration values (jk=2 to jksed) for each couple
[3443]460      DO jk = 2, jpksed
461         DO ji = 1, jpoce
[10222]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
[3443]468         ENDDO
469      ENDDO
[10222]470!
[3443]471      ! New dissolved concentrations
[10222]472      DO jk = 2, jpksed
[3443]473         DO ji = 1, jpoce
[10222]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)
[3443]478            ! For DIC
479            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
[10222]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
[3443]487            ! For alkalinity
[10222]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
[3443]491         ENDDO
492      ENDDO
493
[10222]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      ! ------------------------------------------------------------------------------
[3443]498
[10222]499      call sed_dsr_redoxb
[3443]500
[10222]501      ! --------------------------------------------------------------
502      !    4/ Computation of the bioturbation coefficient
503      !       This parameterization is taken from Archer et al. (2002)
504      ! --------------------------------------------------------------
[3443]505
[10222]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
[3443]509
[10222]510      ! ------------------------------------------------------
511      !    Vertical variations of the bioturbation coefficient
512      ! ------------------------------------------------------
513      IF (ln_btbz) THEN
[3443]514         DO ji = 1, jpoce
[10222]515            db(ji,:) = db(ji,:) * exp( -(profsedw(:) / dbtbzsc)**2 ) / (365.0 * 86400.0)
[3443]516         END DO
[10222]517      ELSE
518         DO jk = 1, jpksed
519            IF (profsedw(jk) > dbtbzsc) THEN
520                db(:,jk) = 0.0
521            ENDIF
[3443]522         END DO
[10222]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
[3443]532         END DO
[10222]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
[3443]575      DO ji = 1, jpoce
[10222]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)
[10362]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
[10222]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)
[10362]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
[10222]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)
[10362]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
[10222]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)
[10362]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
[10222]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)
[10362]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
[10222]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)
[10362]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
[10222]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)
[10362]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
[10222]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)
[10362]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
[10222]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)
[10362]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
[10222]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)
[10362]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
[10222]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)
[10362]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
[10222]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)
[3443]742         END DO
[10222]743      END DO
[3443]744
[10222]745      IF( ln_timing )  CALL timing_stop('sed_dsr_redoxb')
746
747  END SUBROUTINE sed_dsr_redoxb
748
[3443]749END MODULE seddsr
Note: See TracBrowser for help on using the repository browser.