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 @ 15127

Last change on this file since 15127 was 15127, checked in by cetlod, 2 years ago

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