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.
seddsrjac.F90 in NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddsrjac.F90 @ 15297

Last change on this file since 15297 was 15297, checked in by aumont, 2 years ago

major update of the sediment module
File size: 22.2 KB
Line 
1MODULE seddsrjac
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 seddsr
12   USE lib_mpp         ! distribued memory computing library
13   USE lib_fortran
14
15   IMPLICIT NONE
16   PRIVATE
17
18   PUBLIC sed_dsrjac
19
20   !! * Module variables
21
22   !! $Id: seddsr.F90 10362 2018-11-30 15:38:17Z aumont $
23CONTAINS
24   
25   SUBROUTINE sed_dsrjac( NEQ, NROWPD, jacvode, accmask )
26      !!----------------------------------------------------------------------
27      !!                   ***  ROUTINE sed_dsr  ***
28      !!
29      !!  ** Purpose :  computes pore water dissolution and reaction
30      !!
31      !!  ** Methode :  Computation of the redox reactions in sediment.
32      !!                The main redox reactions are solved in sed_dsr whereas
33      !!                the secondary reactions are solved in sed_dsr_redoxb.
34      !!                A strand spliting approach is being used here (see
35      !!                sed_dsr_redoxb for more information).
36      !!
37      !!   History :
38      !!        !  98-08 (E. Maier-Reimer, Christoph Heinze )  Original code
39      !!        !  04-10 (N. Emprin, M. Gehlen ) f90
40      !!        !  06-04 (C. Ethe)  Re-organization
41      !!        !  19-08 (O. Aumont) Debugging and improvement of the model.
42      !!                             The original method is replaced by a
43      !!                              Strand splitting method which deals
44      !!                              well with stiff reactions.
45      !!----------------------------------------------------------------------
46      !! Arguments
47      ! --- local variables
48      INTEGER, INTENT(in) :: NEQ, NROWPD
49      INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask
50      REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode
51      INTEGER :: ji, jni, jnj, jnij, jk, js, jw, jn   ! dummy looop indices
52
53      REAL(wp), DIMENSION(jpoce, jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo    ! undersaturation ; indice jpwatp1 is for calcite   
54      REAL(wp), DIMENSION(jpoce, jpksed) ::zlimdo2, zlimdno3, zlimdso4, zlimdfeo    ! undersaturation ; indice jpwatp1 is for calcite   
55      REAL(wp)  ::  zvolw, zreasat, zlimtmpo2, zlimtmpno3, zlimtmpfeo, zlimtmpso4
56      !!
57      !!----------------------------------------------------------------------
58
59      IF( ln_timing )  CALL timing_start('sed_dsrjac')
60!
61     ! Initializations
62     !----------------------
63     
64      zlimo2 (:,:)    = 0.    ;   zlimno3(:,:)  = 0.
65      zlimso4(:,:)    = 0.
66 
67      !----------------------------------------------------------
68      ! 5.  Beginning of solid reaction
69      !---------------------------------------------------------
70     
71      ! Computes SMS of oxygen
72      DO jk = 2, jpksed
73         DO ji = 1, jpoce
74            IF (accmask(ji) == 0) THEN
75               zlimo2(ji,jk) = pwcp(ji,jk,jwoxy) / ( pwcp(ji,jk,jwoxy) + xksedo2 )
76               zlimdo2(ji,jk) = xksedo2 / ( pwcp(ji,jk,jwoxy) + xksedo2 )**2
77               ! Acid Silicic
78               jni = (jk-1)*jpvode+isvode(jwoxy)
79               jnij = jpvode + 1
80               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - so2ut * rearatpom(ji,jk) * zlimdo2(ji,jk)
81            ENDIF
82         END DO
83      ENDDO
84
85      !--------------------------------------------------------------------
86      ! Denitrification
87      !--------------------------------------------------------------------
88      DO jk = 2, jpksed
89         DO ji = 1, jpoce
90            IF (accmask(ji) == 0) THEN
91
92               zlimno3(ji,jk) = pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 ) 
93               zlimdno3(ji,jk) = xksedno3 / ( pwcp(ji,jk,jwno3) + xksedno3 )**2
94               ! For nitrates
95               jni = (jk-1)*jpvode+isvode(jwno3)
96               jnj = (jk-1)*jpvode+isvode(jwoxy)
97               jnij = jpvode + 1
98               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - srDnit * rearatpom(ji,jk) * (1.0 - zlimo2(ji,jk) ) &
99               &         * zlimdno3(ji,jk)
100               jnij = jni - jnj + jpvode + 1
101               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + srDnit * rearatpom(ji,jk) * zlimno3(ji,jk) * zlimdo2(ji,jk)
102            ENDIF
103         END DO
104      ENDDO
105
106      !--------------------------------------------------------------------
107      ! Begining POC iron reduction
108      !--------------------------------------------------------------------
109
110      DO jk = 2, jpksed
111         DO ji = 1, jpoce
112            IF (accmask(ji) == 0) THEN
113               zlimfeo(ji,jk)  = solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo )
114               zlimdfeo(ji,jk) = xksedfeo / ( solcp(ji,jk,jsfeo) + xksedfeo )**2
115               ! For FEOH
116               jni = (jk-1)*jpvode+isvode(jpwat+jsfeo)
117               jnij = jpvode + 1
118               zlimtmpfeo = ( 1.0 - zlimno3(ji,jk) ) * ( 1.0 - zlimo2(ji,jk) ) * zlimdfeo(ji,jk)
119               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 4.0 * rearatpom(ji,jk) / volc(ji,jk,jsfeo) * zlimtmpfeo
120               jnj = (jk-1)*jpvode+isvode(jwoxy)
121               jnij = jni - jnj + jpvode + 1
122               zlimtmpo2 = zlimfeo(ji,jk) * zlimdo2(ji,jk) * ( 1.0 - zlimno3(ji,jk) )
123               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 4.0 * rearatpom(ji,jk) / volc(ji,jk,jsfeo) * zlimtmpo2 
124               jnj = (jk-1)*jpvode+isvode(jwno3)
125               jnij = jni - jnj + jpvode + 1
126               zlimtmpno3 = zlimfeo(ji,jk) * zlimdno3(ji,jk) * ( 1.0 - zlimo2(ji,jk) )
127               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 4.0 * rearatpom(ji,jk) / volc(ji,jk,jsfeo) * zlimtmpno3
128               ! Iron
129               zreasat = rearatpom(ji,jk) * 4.0 * radssol(jk,jwfe2)
130               jni = (jk-1)*jpvode+isvode(jwfe2)
131               jnj = (jk-1)*jpvode+isvode(jpwat+jsfeo)
132               jnij = jni - jnj + jpvode + 1
133               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + zreasat * zlimtmpfeo
134               jnj = (jk-1)*jpvode+isvode(jwoxy)
135               jnij = jni - jnj + jpvode + 1
136               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - zreasat * zlimtmpo2
137               jnj = (jk-1)*jpvode+isvode(jwno3)
138               jnij = jni - jnj + jpvode + 1
139               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - zreasat * zlimtmpno3
140            ENDIF
141         END DO
142      ENDDO
143
144      !--------------------------------------------------------------------
145      ! Begining POC denitrification and NO3- diffusion
146      !--------------------------------------------------------------------
147
148      DO jk = 2, jpksed
149         DO ji = 1, jpoce
150            IF (accmask(ji) == 0) THEN
151               zlimso4(ji,jk) = pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 )
152               zlimdso4(ji,jk) = xksedso4 / ( pwcp(ji,jk,jwso4) + xksedso4 )**2
153               ! For sulfur
154               jni = (jk-1) * jpvode + isvode(jwso4)
155               jnij = jpvode + 1
156               zlimtmpso4 = ( 1.0 - zlimno3(ji,jk) ) * ( 1.0 - zlimo2(ji,jk) ) * ( 1.0 - zlimfeo(ji,jk) ) * zlimdso4(ji,jk)
157               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 0.5 * rearatpom(ji,jk) * zlimtmpso4
158               jnj = (jk-1) * jpvode + isvode(jwoxy)
159               jnij = jni - jnj + jpvode + 1
160               zlimtmpo2 = zlimso4(ji,jk) * zlimdo2(ji,jk) * ( 1.0 - zlimno3(ji,jk) )  &
161               &      * ( 1.0 - zlimfeo(ji,jk) ) 
162               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 0.5 * rearatpom(ji,jk) * zlimtmpo2
163               jnj = (jk-1) * jpvode + isvode(jwno3)
164               jnij = jni - jnj + jpvode + 1
165               zlimtmpno3 = zlimso4(ji,jk) * ( 1.0 - zlimo2(ji,jk) ) * zlimdno3(ji,jk)   &
166               &      * ( 1.0 - zlimfeo(ji,jk) )
167               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 0.5 * rearatpom(ji,jk) * zlimtmpno3 
168               jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo)
169               jnij = jni - jnj + jpvode + 1
170               zlimtmpfeo = zlimso4(ji,jk) * ( 1.0 - zlimo2(ji,jk) ) * ( 1.0 - zlimno3(ji,jk) ) * zlimdfeo(ji,jk)
171               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 0.5 * rearatpom(ji,jk) * zlimtmpfeo
172               jni = (jk-1) * jpvode + isvode(jwh2s)
173               jnj = (jk-1) * jpvode + isvode(jwso4)
174               jnij = jni - jnj + jpvode + 1
175               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 0.5 * rearatpom(ji,jk) * zlimtmpso4
176               jnj = (jk-1) * jpvode + isvode(jwoxy)
177               jnij = jni - jnj + jpvode + 1
178               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 0.5 * rearatpom(ji,jk) * zlimtmpo2
179               jnj = (jk-1) * jpvode + isvode(jwno3)
180               jnij = jni - jnj + jpvode + 1
181               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 0.5 * rearatpom(ji,jk) * zlimtmpno3
182               jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo)
183               jnij = jni - jnj + jpvode + 1
184               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 0.5 * rearatpom(ji,jk) * zlimtmpfeo
185            ENDIF
186         END DO
187      ENDDO
188
189      ! Secondary redox reactions
190      ! -------------------------
191
192      call sed_dsr_redoxbjac( NEQ, NROWPD, jacvode, accmask )
193
194      IF( ln_timing )  CALL timing_stop('sed_dsrjac')
195!     
196   END SUBROUTINE sed_dsrjac
197
198   SUBROUTINE sed_dsr_redoxbjac( NEQ, NROWPD, jacvode, accmask )
199      !!----------------------------------------------------------------------
200      !!                   ***  ROUTINE sed_dsr_redox  ***
201      !!
202      !!  ** Purpose :  computes secondary redox reactions
203      !!
204      !!   History :
205      !!        !  18-08 (O. Aumont)  Original code
206      !!----------------------------------------------------------------------
207      !! Arguments
208      INTEGER, INTENT(in) :: NEQ, NROWPD
209      REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode
210      INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask
211      ! --- local variables
212      INTEGER   ::  ji, jni, jnj, jnij, jk   ! dummy looop indices
213
214      REAL(wp)  ::  zalpha, zexcess, zh2seq, zsedfer, zdsedfer
215      !!
216      !!----------------------------------------------------------------------
217
218      IF( ln_timing )  CALL timing_start('sed_dsr_redoxbjac')
219
220      DO jk = 2, jpksed
221         DO ji = 1, jpoce
222            IF (accmask(ji) == 0) THEN
223               zalpha = ( pwcp(ji,jk,jwfe2) - pwcp(ji,jk,jwlgw) ) * 1E9
224               zsedfer = ( zalpha + SQRT( zalpha**2 + 1.E-10 ) ) /2.0 * 1E-9
225               zdsedfer = zsedfer*1E9 / SQRT( zalpha**2 +1E-10 ) 
226
227               ! First pass of the scheme. At the end, it is 1st order
228               ! -----------------------------------------------------
229               ! Fe (both adsorbed and solute) + O2
230               jni = (jk-1) * jpvode + isvode(jwfe2)
231               jnij = jpvode + 1
232               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_fe2 * pwcp(ji,jk,jwoxy) * zdsedfer
233               jnj = (jk-1) * jpvode + isvode(jwoxy)
234               jnij = jni - jnj + jpvode + 1
235               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fe2 * zsedfer
236               jni = (jk-1) * jpvode + isvode(jwoxy)
237               jnij = jpvode + 1
238               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 0.25 * reac_fe2 * zsedfer / radssol(jk,jwfe2)
239               jnj = (jk-1) * jpvode + isvode(jwfe2)
240               jnij = jni - jnj + jpvode + 1
241               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 0.25 * reac_fe2 / radssol(jk,jwfe2) * pwcp(ji,jk,jwoxy) * zdsedfer
242               jni = (jk-1) * jpvode + isvode(jpwat+jsfeo)
243               jnj = (jk-1) * jpvode + isvode(jwfe2)
244               jnij = jni - jnj + jpvode + 1
245               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fe2 / radssol(jk,jwfe2) * pwcp(ji,jk,jwoxy)   &
246               &     * zdsedfer / volc(ji,jk,jsfeo)
247               jnj = (jk-1) * jpvode + isvode(jwoxy)
248               jnij = jni - jnj + jpvode + 1
249               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fe2 / radssol(jk,jwfe2) * zsedfer  &
250               &     / volc(ji,jk,jsfeo)
251
252               ! H2S + O2
253               jni = (jk-1) * jpvode + isvode(jwh2s)
254               jnij = jpvode + 1
255               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_h2s * pwcp(ji,jk,jwoxy)
256               jnj = (jk-1) * jpvode + isvode(jwoxy)
257               jnij = jni - jnj + jpvode + 1
258               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_h2s * pwcp(ji,jk,jwh2s)
259               jni = (jk-1) * jpvode + isvode(jwoxy)
260               jnij = jpvode + 1
261               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 2.0 * reac_h2s * pwcp(ji,jk,jwh2s)
262               jnj = (jk-1) * jpvode + isvode(jwh2s)
263               jnij = jni - jnj + jpvode + 1
264               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 2.0 * reac_h2s * pwcp(ji,jk,jwoxy)
265               jni = (jk-1) * jpvode + isvode(jwso4)
266               jnj = (jk-1) * jpvode + isvode(jwh2s)
267               jnij = jni - jnj + jpvode + 1
268               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_h2s * pwcp(ji,jk,jwoxy)
269               jnj = (jk-1) * jpvode + isvode(jwoxy)
270               jnij = jni - jnj + jpvode + 1
271               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_h2s * pwcp(ji,jk,jwh2s)
272
273               ! NH4 + O2
274               jni = (jk-1) * jpvode + isvode(jwnh4)
275               jnij = jpvode + 1
276               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_nh4 * pwcp(ji,jk,jwoxy)
277               jnj = (jk-1) * jpvode + isvode(jwoxy)
278               jnij = jni - jnj + jpvode + 1
279               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_nh4 * pwcp(ji,jk,jwnh4)
280               jni = (jk-1) * jpvode + isvode(jwoxy)
281               jnij = jpvode + 1
282               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 2.0 * reac_nh4 * pwcp(ji,jk,jwnh4) / radssol(jk,jwnh4)
283               jnj = (jk-1) * jpvode + isvode(jwnh4)
284               jnij = jni - jnj + jpvode + 1
285               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 2.0 * reac_nh4 * pwcp(ji,jk,jwoxy) / radssol(jk,jwnh4)
286               jni = (jk-1) * jpvode + isvode(jwno3)
287               jnj = (jk-1) * jpvode + isvode(jwoxy)
288               jnij = jni - jnj + jpvode + 1
289               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_nh4 * pwcp(ji,jk,jwnh4) / radssol(jk,jwnh4)
290               jnj = (jk-1) * jpvode + isvode(jwnh4)
291               jnij = jni - jnj + jpvode + 1
292               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_nh4 * pwcp(ji,jk,jwoxy) / radssol(jk,jwnh4)
293
294               ! FeS - O2
295               jni = (jk-1) * jpvode + isvode(jpwat+jsfes)
296               jnij = jpvode + 1
297               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_feso * pwcp(ji,jk,jwoxy) 
298               jnj = (jk-1) * jpvode + isvode(jwoxy)
299               jnij = jni - jnj + jpvode + 1
300               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_feso * solcp(ji,jk,jsfes)
301               jni = (jk-1) * jpvode + isvode(jwoxy)
302               jnij = jpvode + 1
303               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 2.0 * reac_feso * solcp(ji,jk,jsfes)  &
304               &     * volc(ji,jk,jsfes)
305               jnj = (jk-1) * jpvode + isvode(jpwat+jsfes)
306               jnij = jni - jnj + jpvode + 1
307               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 2.0 * reac_feso * pwcp(ji,jk,jwoxy)   &
308               &     * volc(ji,jk,jsfes)
309               jni = (jk-1) * jpvode + isvode(jwfe2)
310               jnj = (jk-1) * jpvode + isvode(jwoxy)
311               jnij = jni - jnj + jpvode + 1
312               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feso * solcp(ji,jk,jsfes)  &
313               &     * volc(ji,jk,jsfes) * radssol(jk,jwfe2)
314               jnj = (jk-1) * jpvode + isvode(jpwat+jsfes)
315               jnij = jni - jnj + jpvode + 1
316               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feso * pwcp(ji,jk,jwoxy)   &
317               &     * volc(ji,jk,jsfes) * radssol(jk,jwfe2)
318               jni = (jk-1) * jpvode + isvode(jwso4)
319               jnj = (jk-1) * jpvode + isvode(jwoxy)
320               jnij = jni - jnj + jpvode + 1
321               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feso * solcp(ji,jk,jsfes)  &
322               &     * volc(ji,jk,jsfes)
323               jnj = (jk-1) * jpvode + isvode(jpwat+jsfes)
324               jnij = jni - jnj + jpvode + 1
325               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feso * pwcp(ji,jk,jwoxy)   &
326               &     * volc(ji,jk,jsfes)
327
328               ! FEOH + H2S
329               jni = (jk-1) * jpvode + isvode(jpwat+jsfeo)
330               jnij = jpvode + 1
331               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 8.0 * reac_feh2s * pwcp(ji,jk,jwh2s)
332               jnj = (jk-1) * jpvode + isvode(jwh2s)
333               jnij = jni - jnj + jpvode + 1
334               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 8.0 * reac_feh2s * solcp(ji,jk,jsfeo)
335               jni = (jk-1) * jpvode + isvode(jwh2s)
336               jnij = jpvode + 1
337               jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_feh2s * solcp(ji,jk,jsfeo)  &
338               &     * volc(ji,jk,jsfeo)
339               jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo)
340               jnij = jni - jnj + jpvode + 1
341               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_feh2s * pwcp(ji,jk,jwh2s)   &
342               &     * volc(ji,jk,jsfeo)
343               jni = (jk-1) * jpvode + isvode(jwfe2)
344               jnj = (jk-1) * jpvode + isvode(jwh2s)
345               jnij = jni - jnj + jpvode + 1
346               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 8.0 * reac_feh2s * solcp(ji,jk,jsfeo)  &
347               &     * volc(ji,jk,jsfeo) * radssol(jk,jwfe2)
348               jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo)
349               jnij = jni - jnj + jpvode + 1
350               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 8.0 * reac_feh2s * pwcp(ji,jk,jwh2s)   &
351               &     * volc(ji,jk,jsfeo) * radssol(jk,jwfe2)
352               jni = (jk-1) * jpvode + isvode(jwso4)
353               jnj = (jk-1) * jpvode + isvode(jwh2s)
354               jnij = jni - jnj + jpvode + 1
355               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feh2s * solcp(ji,jk,jsfeo)  &
356               &     * volc(ji,jk,jsfeo)
357               jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo)
358               jnij = jni - jnj + jpvode + 1
359               jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feh2s * pwcp(ji,jk,jwh2s)   &
360               &     * volc(ji,jk,jsfeo)
361
362               ! Fe + H2S
363               zh2seq     = MAX(rtrn, 2.1E-3 * hipor(ji,jk))
364               zexcess = pwcp(ji,jk,jwh2s) * zsedfer / zh2seq - 1.0
365               IF ( zexcess >= 0.0 ) THEN
366                  jni = (jk-1) * jpvode + isvode(jwfe2)
367                  jnij = jpvode + 1
368                  jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_fesp * pwcp(ji,jk,jwh2s) * zdsedfer / zh2seq * radssol(jk,jwfe2)
369                  jnj = (jk-1) * jpvode + isvode(jwh2s)
370                  jnij = jni - jnj + jpvode + 1
371                  jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesp * zsedfer / zh2seq * radssol(jk,jwfe2)
372                  jni = (jk-1) * jpvode + isvode(jpwat+jsfes)
373                  jnj = (jk-1) * jpvode + isvode(jwfe2)
374                  jnij = jni - jnj + jpvode + 1
375                  jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fesp * pwcp(ji,jk,jwh2s) / zh2seq   &
376                  &     * zdsedfer / volc(ji,jk,jsfes)
377                  jnj = (jk-1) * jpvode + isvode(jwh2s)
378                  jnij = jni - jnj + jpvode + 1
379                  jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fesp * zsedfer / zh2seq   &
380                  &     / volc(ji,jk,jsfes)
381                  jni = (jk-1) * jpvode + isvode(jwh2s)
382                  jnij = jpvode + 1
383                  jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_fesp * zsedfer / zh2seq 
384                  jnj = (jk-1) * jpvode + isvode(jwfe2)
385                  jnij = jni - jnj + jpvode + 1
386                  jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesp * pwcp(ji,jk,jwh2s) * zdsedfer / zh2seq 
387               ELSE
388                  jni = (jk-1) * jpvode + isvode(jwfe2)
389                  jnij = jpvode + 1
390                  jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_fesd * pwcp(ji,jk,jwh2s) / zh2seq   &
391                  &     * zdsedfer * radssol(jk,jwfe2) * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes)
392                  jnj = (jk-1) * jpvode + isvode(jwh2s)
393                  jnij = jni - jnj + jpvode + 1
394                  jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesd * zsedfer / zh2seq * radssol(jk,jwfe2)   &
395                  &     * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes)
396                  jnj = (jk-1) * jpvode + isvode(jpwat+jsfes)
397                  jnij = jni - jnj + jpvode + 1
398                  jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesd * zexcess * radssol(jk,jwfe2) * volc(ji,jk,jsfes)
399                  jni = (jk-1) * jpvode + isvode(jpwat+jsfes)
400                  jnij = jpvode + 1
401                  jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) + reac_fesd * zexcess 
402                  jnj = (jk-1) * jpvode + isvode(jwfe2)
403                  jnij = jni - jnj + jpvode + 1
404                  jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fesd * solcp(ji,jk,jsfes)   &
405                  &     * zdsedfer * pwcp(ji,jk,jwh2s) / zh2seq 
406                  jnj = (jk-1) * jpvode + isvode(jwh2s)
407                  jnij = jni - jnj + jpvode + 1
408                  jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fesd * solcp(ji,jk,jsfes)   &
409                  &     * zsedfer / zh2seq
410                  jni = (jk-1) * jpvode + isvode(jwh2s)
411                  jnij = jpvode + 1
412                  jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_fesd * zsedfer / zh2seq    &
413                  &     * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes)
414                  jnj = (jk-1) * jpvode + isvode(jwfe2)
415                  jnij = jni - jnj + jpvode + 1
416                  jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesd * pwcp(ji,jk,jwh2s) / zh2seq   &
417                  &     * zdsedfer * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes)
418                  jnj = (jk-1) * jpvode + isvode(jpwat+jsfes)
419                  jnij = jni - jnj + jpvode + 1
420                  jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesd * zexcess * volc(ji,jk,jsfes)
421               ENDIF
422            ENDIF
423         END DO
424     END DO
425
426      IF( ln_timing )  CALL timing_stop('sed_dsr_redoxbjac')
427
428  END SUBROUTINE sed_dsr_redoxbjac
429
430END MODULE seddsrjac
Note: See TracBrowser for help on using the repository browser.