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.
sedmat.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/sedmat.F90 @ 15352

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

dev_PISCO : minor corrections

  • Property svn:keywords set to Id
File size: 24.7 KB
Line 
1MODULE sedmat
2   !!======================================================================
3   !!              ***  MODULE  sedmat  ***
4   !!    Sediment : linear system of equations
5   !!=====================================================================
6   !! * Modules used
7   !!----------------------------------------------------------------------
8
9   USE sed     ! sediment global variable
10   USE lib_mpp         ! distribued memory computing library
11
12
13   IMPLICIT NONE
14   PRIVATE
15
16   PUBLIC sed_mat_dsr 
17   PUBLIC sed_mat_dsrjac
18   PUBLIC sed_mat_dsri
19   PUBLIC sed_mat_btb
20   PUBLIC sed_mat_btbjac
21   PUBLIC sed_mat_btbi
22   PUBLIC sed_mat_coef
23
24   !! $Id$
25 CONTAINS
26
27    SUBROUTINE sed_mat_coef( nksed )
28       !!---------------------------------------------------------------------
29       !!                  ***  ROUTINE sed_mat_coef  ***
30       !!
31       !! ** Purpose :  solves tridiagonal system of linear equations
32       !!
33       !! ** Method  :
34       !!        1 - computes left hand side of linear system of equations
35       !!            for dissolution reaction
36       !!         For mass balance in kbot+sediment :
37       !!              dz3d  (:,1) = dz(1) = 0.5 cm
38       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )
39       !!              dz(2)       = 0.3 cm
40       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )     
41       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see seddsr.F90 )
42       !!
43       !!         2 - forward/backward substitution.
44       !!
45       !!   History :
46       !!        !  04-10 (N. Emprin, M. Gehlen ) original
47       !!        !  06-04 (C. Ethe)  Module Re-organization
48       !!----------------------------------------------------------------------
49       !! * Arguments
50       INTEGER, INTENT(in) :: nksed
51
52       !---Local declarations
53       INTEGER  ::  ji, jk
54       REAL(wp) ::  aplus, aminus, dxplus, dxminus
55       !----------------------------------------------------------------------
56
57       IF( ln_timing )  CALL timing_start('sed_mat_coef')
58
59       ! Computation left hand side of linear system of
60       ! equations for dissolution reaction
61       !---------------------------------------------
62       ! first sediment level         
63       DO ji = 1, jpoce
64          aplus  = ( por(1) + por(2) ) * 0.5
65          dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2.
66          apluss(ji,1) = ( 1.0 / ( volw3d(ji,1) ) ) * aplus / dxplus
67
68          DO jk = 2, nksed - 1
69             aminus  = ( por(jk-1) + por(jk) ) * 0.5
70             dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2.
71
72             aplus   = ( por(jk+1) + por(jk) ) * 0.5
73             dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2
74             !
75             aminuss(ji,jk) = ( 1.0 / volw3d(ji,jk) ) * aminus / dxminus
76             apluss (ji,jk) = ( 1.0 / volw3d(ji,jk) ) * aplus / dxplus
77          END DO
78
79          aminus  = ( por(nksed-1) + por(nksed) ) * 0.5
80          dxminus = ( dz3d(ji,nksed-1) + dz3d(ji,nksed) ) / 2.
81          aminuss(ji,nksed)  = ( 1.0 / volw3d(ji,nksed) ) * aminus / dxminus
82
83       END DO
84
85       IF( ln_timing )  CALL timing_stop('sed_mat_coef')
86
87    END SUBROUTINE sed_mat_coef
88
89    SUBROUTINE sed_mat_dsr( nksed, nvar, accmask )
90       !!---------------------------------------------------------------------
91       !!                  ***  ROUTINE sed_mat_dsr  ***
92       !!
93       !! ** Purpose :  solves tridiagonal system of linear equations
94       !!
95       !! ** Method  :
96       !!        1 - computes left hand side of linear system of equations
97       !!            for dissolution reaction
98       !!         For mass balance in kbot+sediment :
99       !!              dz3d  (:,1) = dz(1) = 0.5 cm
100       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )
101       !!              dz(2)       = 0.3 cm
102       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )     
103       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see seddsr.F90 )
104       !!
105       !!         2 - forward/backward substitution.
106       !!
107       !!   History :
108       !!        !  04-10 (N. Emprin, M. Gehlen ) original
109       !!        !  06-04 (C. Ethe)  Module Re-organization
110       !!----------------------------------------------------------------------
111       !! * Arguments
112       INTEGER , INTENT(in) ::  nvar, nksed  ! number of variable
113       INTEGER, DIMENSION(jpoce) :: accmask
114       INTEGER :: ji
115
116       !---Local declarations
117       INTEGER  ::  jk, jn
118       REAL(wp), DIMENSION(nksed) :: za, zb, zc
119
120       REAL(wp) ::  rplus,rminus   
121       !----------------------------------------------------------------------
122
123       IF( ln_timing )  CALL timing_start('sed_mat_dsr')
124
125       ! Computation left hand side of linear system of
126       ! equations for dissolution reaction
127       !---------------------------------------------
128       jn = nvar
129       ! first sediment level         
130       DO ji = 1, jpoce
131          IF (accmask(ji) == 0) THEN
132             rplus  = apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn)
133
134             za(1) = 0.
135             zb(1) = rplus
136             zc(1) = -rplus
137 
138             DO jk = 2, nksed - 1
139                rminus  = aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn)
140                rplus   = apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn)
141                !     
142                za(jk) = -rminus
143                zb(jk) = rminus + rplus 
144                zc(jk) = -rplus
145             END DO
146
147             rminus  = aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn)
148             !
149             za(nksed) = -rminus
150             zb(nksed) = rminus
151             zc(nksed) = 0.
152
153             ! solves tridiagonal system of linear equations
154             ! -----------------------------------------------
155
156             pwcpa(ji,1,jn) = pwcpa(ji,1,jn) - ( zc(1) * pwcp(ji,2,jn) + zb(1) * pwcp(ji,1,jn) )
157             DO jk = 2, nksed - 1
158                pwcpa(ji,jk,jn) =  pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn)    &
159                &                  + zb(jk) * pwcp(ji,jk,jn) )
160             ENDDO
161             pwcpa(ji,nksed,jn) = pwcpa(ji,nksed,jn) - ( za(nksed) * pwcp(ji,nksed-1,jn)    &
162             &                     + zb(nksed) * pwcp(ji,nksed,jn) )
163
164          ENDIF
165       END DO
166
167       IF( ln_timing )  CALL timing_stop('sed_mat_dsr')
168
169    END SUBROUTINE sed_mat_dsr
170
171    SUBROUTINE sed_mat_dsrjac( nksed, nvar, NEQ, NROWPD, jacvode, accmask )
172       !!---------------------------------------------------------------------
173       !!                  ***  ROUTINE sed_mat_dsrjac  ***
174       !!
175       !! ** Purpose :  solves tridiagonal system of linear equations
176       !!
177       !! ** Method  :
178       !!        1 - computes left hand side of linear system of equations
179       !!            for dissolution reaction
180       !!         For mass balance in kbot+sediment :
181       !!              dz3d  (:,1) = dz(1) = 0.5 cm
182       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )
183       !!              dz(2)       = 0.3 cm
184       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )     
185       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see
186       !seddsr.F90 )
187       !!
188       !!         2 - forward/backward substitution.
189       !!
190       !!   History :
191       !!        !  04-10 (N. Emprin, M. Gehlen ) original
192       !!        !  06-04 (C. Ethe)  Module Re-organization
193       !!----------------------------------------------------------------------
194       !! * Arguments
195       INTEGER , INTENT(in) ::  nvar, nksed, NEQ, NROWPD  ! number of variable
196       REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode
197       INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask
198
199       !---Local declarations
200       INTEGER  ::  ji,jk, jn, jnn, jni, jnj ,jnij
201       REAL(wp), DIMENSION(nksed) :: za, zb, zc
202
203       REAL(wp) ::  rplus,rminus
204       !----------------------------------------------------------------------
205
206       IF( ln_timing )  CALL timing_start('sed_mat_dsrjac')
207
208       ! Computation left hand side of linear system of
209       ! equations for dissolution reaction
210       !---------------------------------------------
211       jn = nvar
212       ! first sediment level         
213       DO ji = 1, jpoce
214          IF (accmask(ji) == 0 ) THEN
215             rplus  = apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn)
216
217             za(1) = 0.
218             zb(1) = rplus
219             zc(1) = -rplus
220
221             DO jk = 2, nksed - 1
222                rminus  = aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn)
223                rplus   = apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn)
224                !     
225                za(jk) = -rminus
226                zb(jk) = rminus + rplus
227                zc(jk) = -rplus
228             END DO
229
230             rminus  = aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn)
231             !
232             za(nksed) = -rminus
233             zb(nksed) = rminus
234             zc(nksed) = 0.
235
236             ! solves tridiagonal system of linear equations
237
238             jnn = isvode(jn)
239             jnij = jpvode + 1
240             jacvode(ji, jnij, jnn) = jacvode(ji,jnij,jnn) - zb(1)
241             jnj = jpvode + jnn
242             jnij = jnn - jnj + jpvode + 1
243             jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) -zc(1)
244             DO jk = 2, nksed - 1
245                jni = (jk-1) * jpvode + jnn
246                jnj = (jk-2) * jpvode + jnn
247                jnij = jni - jnj + jpvode + 1
248                jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jk)
249                jnj = (jk-1) * jpvode + jnn
250                jnij = jni - jnj + jpvode + 1
251                jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zb(jk)
252                jnj = (jk) * jpvode + jnn
253                jnij = jni - jnj + jpvode + 1
254                jacvode(ji, jnij, jnj)   = jacvode(ji, jnij, jnj) - zc(jk)
255             END DO
256             jni = (nksed-1) * jpvode + jnn
257             jnj = (nksed-2) * jpvode + jnn
258             jnij = jni - jnj + jpvode + 1
259             jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(nksed)
260             jnij = jpvode + 1
261             jacvode(ji, jnij, jni) = jacvode(ji, jnij, jni) - zb(nksed)
262          ENDIF
263       END DO
264
265       IF( ln_timing )  CALL timing_stop('sed_mat_dsrjac')
266
267    END SUBROUTINE sed_mat_dsrjac
268
269    SUBROUTINE sed_mat_btbi( nksed, nvar, psol, preac, dtsed_in )
270       !!---------------------------------------------------------------------
271       !!                  ***  ROUTINE sed_mat_btb  ***
272       !!
273       !! ** Purpose :  solves tridiagonal system of linear equations
274       !!
275       !! ** Method  :
276       !!        1 - computes left hand side of linear system of equations
277       !!            for dissolution reaction
278       !!
279       !!
280       !!         2 - forward/backward substitution.
281       !!
282       !!   History :
283       !!        !  04-10 (N. Emprin, M. Gehlen ) original
284       !!        !  06-04 (C. Ethe)  Module Re-organization
285       !!----------------------------------------------------------------------
286       !! * Arguments
287       INTEGER , INTENT(in) :: nksed, nvar      ! number of sediment levels
288
289      REAL(wp), DIMENSION(jpoce,nksed,nvar), INTENT(inout) :: &
290          psol, preac
291
292      REAL(wp), INTENT(in) :: dtsed_in
293
294       !---Local declarations
295       INTEGER  ::  &
296          ji, jk, jn
297
298       REAL(wp) ::  &
299          aplus,aminus   ,  &
300          rplus,rminus   ,  &
301          dxplus,dxminus
302
303       REAL(wp), DIMENSION(nksed)    :: za, zb, zc
304       REAL(wp), DIMENSION(nksed)    :: zr, zgamm
305       REAL(wp) ::  zbet
306
307       !----------------------------------------------------------------------
308
309      ! Computation left hand side of linear system of
310      ! equations for dissolution reaction
311      !---------------------------------------------
312      IF( ln_timing )  CALL timing_start('sed_mat_btbi')
313
314      ! first sediment level         
315      DO ji = 1, jpoce
316         aplus  = ( por1(2) + por1(3) ) / 2.0
317         dxplus = ( dz(2) + dz(3) ) / 2.
318         rplus  = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus
319         za(2) = 0.
320         zb(2) = 1. + rplus
321         zc(2) = -rplus
322
323         DO jk = 3, nksed - 1
324            aminus  = ( por1(jk-1) + por1(jk) ) * 0.5
325            aminus  = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2.
326            dxminus = ( dz(jk-1) + dz(jk) ) / 2.
327            rminus  = ( dtsed_in / vols(jk) ) * db(ji,jk-1) * aminus / dxminus
328            !
329            aplus   = ( por1(jk) + por1(jk+1) ) * 0.5
330            dxplus  = ( dz(jk) + dz(jk+1) ) / 2.
331            rplus   = ( dtsed_in / vols(jk) ) * db(ji,jk) * aplus / dxplus
332            !     
333            za(jk) = -rminus
334            zb(jk) = 1. + rminus + rplus
335            zc(jk) = -rplus
336
337         ENDDO
338
339         aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5
340         dxminus = ( dz(nksed-1) + dz(nksed) ) / 2.
341         rminus  = ( dtsed_in / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus
342         !
343         za(nksed) = -rminus
344         zb(nksed) = 1. + rminus
345         zc(nksed) = 0.
346
347         ! solves tridiagonal system of linear equations
348         ! -----------------------------------------------   
349         DO jn = 1, nvar
350            zr(:) = psol(ji,:,jn)
351            zbet     = zb(2) - preac(ji,2,jn) * dtsed_in
352            psol(ji,2,jn) = zr(2) / zbet
353            !
354            DO jk = 3, nksed
355               zgamm(jk) =  zc(jk-1) / zbet
356               zbet      =  zb(jk) - preac(ji,jk,jn) * dtsed_in - za(jk) * zgamm(jk)
357               psol(ji,jk,jn) = ( zr(jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet
358            ENDDO
359            !
360            DO jk = nksed - 1, 2, -1
361               psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn)
362            ENDDO
363         END DO
364      END DO
365      !
366      IF( ln_timing )  CALL timing_stop('sed_mat_btbi')
367
368    END SUBROUTINE sed_mat_btbi
369
370
371    SUBROUTINE sed_mat_btb( nksed, nvar, accmask )
372       !!---------------------------------------------------------------------
373       !!                  ***  ROUTINE sed_mat_btb  ***
374       !!
375       !! ** Purpose :  solves tridiagonal system of linear equations
376       !!
377       !! ** Method  :
378       !!        1 - computes left hand side of linear system of equations
379       !!            for dissolution reaction
380       !!
381       !!         2 - forward/backward substitution.
382       !!
383       !!   History :
384       !!        !  04-10 (N. Emprin, M. Gehlen ) original
385       !!        !  06-04 (C. Ethe)  Module Re-organization
386       !!----------------------------------------------------------------------
387       !! * Arguments
388       INTEGER , INTENT(in) :: &
389          nvar, nksed     ! number of sediment levels
390       INTEGER, DIMENSION(jpoce) :: accmask
391
392       !---Local declarations
393       INTEGER  ::  ji, jk, jn
394
395       REAL(wp) ::  &
396          aplus,aminus   ,  &
397          rplus,rminus   ,  &
398          dxplus,dxminus
399
400       REAL(wp), DIMENSION(nksed)      :: za, zb, zc
401
402       !----------------------------------------------------------------------
403
404      ! Computation left hand side of linear system of
405      ! equations for dissolution reaction
406      !---------------------------------------------
407      IF( ln_timing )  CALL timing_start('sed_mat_btb')
408
409      ! first sediment level         
410      jn = nvar
411      DO ji = 1, jpoce
412         IF (accmask(ji) == 0) THEN
413         aplus  = ( por1(2) + por1(3) ) / 2.0
414         dxplus = ( dz(2) + dz(3) ) / 2.
415         rplus  = ( 1.0 / vols(2) ) * db(ji,2) * aplus / dxplus * rads1sol(2,jn)
416
417         za(2) = 0.
418         zb(2) = rplus 
419         zc(2) = -rplus
420
421         DO jk = 3, nksed - 1
422            aminus  = ( por1(jk-1) + por1(jk) ) * 0.5
423            aminus  = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2.
424            dxminus = ( dz(jk-1) + dz(jk) ) / 2.
425            rminus  = ( 1.0 / vols(jk) ) * db(ji,jk-1) * aminus / dxminus * rads1sol(jk,jn)
426            !
427            aplus   = ( por1(jk) + por1(jk+1) ) * 0.5
428            dxplus  = ( dz(jk) + dz(jk+1) ) / 2.
429            rplus   = ( 1.0 / vols(jk) ) * db(ji,jk) * aplus / dxplus * rads1sol(jk,jn)
430            !     
431            za(jk) = -rminus
432            zb(jk) = rminus + rplus
433            zc(jk) = -rplus
434
435         ENDDO
436
437         aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5
438         dxminus = ( dz(nksed-1) + dz(nksed) ) / 2.
439         rminus  = ( 1.0 / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus * rads1sol(nksed,jn)
440         !
441         za(nksed) = -rminus
442         zb(nksed) = rminus
443         zc(nksed) = 0.
444
445         ! solves tridiagonal system of linear equations
446         ! -----------------------------------------------   
447         pwcpa(ji,2,jn) = pwcpa(ji,2,jn) - ( zc(2) * pwcp(ji,3,jn) + zb(2) * pwcp(ji,2,jn) )
448         DO jk = 3, nksed-1
449            pwcpa(ji,jk,jn) =  pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn)    &
450            &                  + zb(jk) * pwcp(ji,jk,jn) )
451         ENDDO
452         pwcpa(ji,nksed,jn) = pwcpa(ji,nksed,jn) - ( za(nksed) * pwcp(ji,nksed-1,jn)    &
453         &                     + zb(nksed) * pwcp(ji,nksed,jn) )
454         !
455         ENDIF
456      END DO
457      !
458      IF( ln_timing )  CALL timing_stop('sed_mat_btb')
459       
460    END SUBROUTINE sed_mat_btb
461
462    SUBROUTINE sed_mat_btbjac( nksed, nvar, NEQ, NROWPD, jacvode, accmask )
463       !!---------------------------------------------------------------------
464       !!                  ***  ROUTINE sed_mat_btb  ***
465       !!
466       !! ** Purpose :  solves tridiagonal system of linear equations
467       !!
468       !! ** Method  :
469       !!        1 - computes left hand side of linear system of equations
470       !!            for dissolution reaction
471       !!
472       !!         2 - forward/backward substitution.
473       !!
474       !!   History :
475       !!        !  04-10 (N. Emprin, M. Gehlen ) original
476       !!        !  06-04 (C. Ethe)  Module Re-organization
477       !!----------------------------------------------------------------------
478       !! * Arguments
479       INTEGER , INTENT(in) ::  nvar, nksed, NEQ, NROWPD  ! number of variable
480       REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode
481       INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask
482
483       !---Local declarations
484       INTEGER  ::  ji, jk, jn, jnn, jni, jnj ,jnij
485
486       REAL(wp) ::  &
487          aplus,aminus   ,  &
488          rplus,rminus   ,  &
489          dxplus,dxminus
490
491       REAL(wp), DIMENSION(nksed)      :: za, zb, zc
492
493       !----------------------------------------------------------------------
494
495      ! Computation left hand side of linear system of
496      ! equations for dissolution reaction
497      !---------------------------------------------
498      IF( ln_timing )  CALL timing_start('sed_mat_btbjac')
499
500      ! first sediment level         
501      jn = nvar
502      DO ji = 1, jpoce
503         IF (accmask(ji) == 0) THEN
504         aplus  = ( por1(2) + por1(3) ) / 2.0
505         dxplus = ( dz(2) + dz(3) ) / 2.
506         rplus  = ( 1.0 / vols(2) ) * db(ji,2) * aplus / dxplus * rads1sol(2,jn)
507
508         za(2) = 0.
509         zb(2) = rplus
510         zc(2) = -rplus
511
512         DO jk = 3, nksed - 1
513            aminus  = ( por1(jk-1) + por1(jk) ) * 0.5
514            aminus  = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2.
515            dxminus = ( dz(jk-1) + dz(jk) ) / 2.
516            rminus  = ( 1.0 / vols(jk) ) * db(ji,jk-1) * aminus / dxminus * rads1sol(jk,jn)
517            !
518            aplus   = ( por1(jk) + por1(jk+1) ) * 0.5
519            dxplus  = ( dz(jk) + dz(jk+1) ) / 2.
520            rplus   = ( 1.0 / vols(jk) ) * db(ji,jk) * aplus / dxplus * rads1sol(jk,jn)
521            !     
522            za(jk) = -rminus
523            zb(jk) = rminus + rplus
524            zc(jk) = -rplus
525
526         ENDDO
527
528         aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5
529         dxminus = ( dz(nksed-1) + dz(nksed) ) / 2.
530         rminus  = ( 1.0 / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus * rads1sol(nksed,jn)
531         !
532         za(nksed) = -rminus
533         zb(nksed) = rminus
534         zc(nksed) = 0.
535
536         ! solves tridiagonal system of linear equations
537         ! -----------------------------------------------   
538         jnn = isvode(jn)
539         jni = jpvode + jnn
540         jnij = jpvode + 1
541         jacvode(ji, jnij, jni) = jacvode(ji,jnij,jni) - zb(2)
542         jnj = 2 * jpvode + jnn
543         jnij = jni - jnj + jpvode + 1
544         jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) -zc(2)
545         DO jk = 3, nksed-1
546            jni = (jk-1) * jpvode + jnn
547            jnj = (jk-2) * jpvode + jnn
548            jnij = jni - jnj + jpvode + 1
549            jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jk)
550            jnj = (jk-1) * jpvode + jnn
551            jnij = jni - jnj + jpvode + 1
552            jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zb(jk)
553            jnj = (jk) * jpvode + jnn
554            jnij = jni - jnj + jpvode + 1
555            jacvode(ji, jnij, jnj)   = jacvode(ji, jnij, jnj) - zc(jk)
556         ENDDO
557         jni = (nksed-1) * jpvode + jnn
558         jnj = (nksed-2) * jpvode + jnn
559         jnij = jni - jnj + jpvode + 1
560         jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(nksed)
561         jnij = jpvode + 1
562         jacvode(ji, jnij, jni) = jacvode(ji, jnij, jni) - zb(nksed)
563         !
564         ENDIF
565      END DO
566      !
567      IF( ln_timing )  CALL timing_stop('sed_mat_btbjac')
568
569    END SUBROUTINE sed_mat_btbjac
570
571
572    SUBROUTINE sed_mat_dsri( nksed, nvar, preac, psms, dtsed_in, psol )
573       !!---------------------------------------------------------------------
574       !!                  ***  ROUTINE sed_mat_dsr  ***
575       !!
576       !! ** Purpose :  solves tridiagonal system of linear equations
577       !!
578       !! ** Method  :
579       !!        1 - computes left hand side of linear system of equations
580       !!            for dissolution reaction
581       !!         For mass balance in kbot+sediment :
582       !!              dz3d  (:,1) = dz(1) = 0.5 cm
583       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )
584       !!              dz(2)       = 0.3 cm
585       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )     
586       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see
587       !seddsr.F90 )
588       !!
589       !!         2 - forward/backward substitution.
590       !!
591       !!   History :
592       !!        !  04-10 (N. Emprin, M. Gehlen ) original
593       !!        !  06-04 (C. Ethe)  Module Re-organization
594       !!----------------------------------------------------------------------
595       !! * Arguments
596       INTEGER , INTENT(in) ::  nksed, nvar  ! number of variable
597
598       REAL(wp), DIMENSION(jpoce,nksed), INTENT(in   ) :: preac  ! reaction rates
599       REAL(wp), DIMENSION(jpoce,nksed), INTENT(in   ) :: psms  ! reaction rates
600       REAL(wp), DIMENSION(jpoce,nksed), INTENT(inout) :: psol  ! reaction rates
601       REAL(wp), INTENT(in) ::  dtsed_in
602
603       !---Local declarations
604       INTEGER  ::  ji, jk, jn
605       REAL(wp), DIMENSION(jpoce,nksed) :: za, zb, zc, zr
606       REAL(wp), DIMENSION(jpoce)        :: zbet
607       REAL(wp), DIMENSION(jpoce,nksed) :: zgamm
608
609       REAL(wp) ::  rplus,rminus
610       !----------------------------------------------------------------------
611
612       IF( ln_timing )  CALL timing_start('sed_mat_dsri')
613
614       ! Computation left hand side of linear system of
615       ! equations for dissolution reaction
616       !---------------------------------------------
617       jn = nvar
618       ! first sediment level         
619       DO ji = 1, jpoce
620          rplus  = dtsed_in * apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn)
621
622          za(ji,1) = 0.
623          zb(ji,1) = 1. + rplus
624          zc(ji,1) = -rplus
625       ENDDO
626
627       DO jk = 2, nksed - 1
628          DO ji = 1, jpoce
629             rminus  = dtsed_in * aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn)
630             rplus   = dtsed_in * apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn)
631                !     
632             za(ji,jk) = -rminus
633             zb(ji,jk) = 1. + rminus + rplus
634             zc(ji,jk) = -rplus
635          END DO
636       END DO
637
638       DO ji = 1, jpoce
639          rminus  = dtsed_in * aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn)
640          !
641          za(ji,nksed) = -rminus
642          zb(ji,nksed) = 1. + rminus
643          zc(ji,nksed) = 0.
644       END DO
645
646
647       ! solves tridiagonal system of linear equations
648       ! -----------------------------------------------
649
650       zr  (:,:) = psol(:,:) + psms(:,:) * dtsed_in
651       zb  (:,:) = zb(:,:) - preac(:,:) * dtsed_in
652       zbet(:  ) = zb(:,1)
653       psol(:,1) = zr(:,1) / zbet(:)
654
655          !
656       DO jk = 2, nksed
657          DO ji = 1, jpoce
658             zgamm(ji,jk) =  zc(ji,jk-1) / zbet(ji)
659             zbet(ji)     =  zb(ji,jk) - za(ji,jk) * zgamm(ji,jk)
660             psol(ji,jk)  = ( zr(ji,jk) - za(ji,jk) * psol(ji,jk-1) ) / zbet(ji)
661          END DO
662       ENDDO
663          !
664       DO jk = nksed - 1, 1, -1
665          DO ji = 1, jpoce
666             psol(ji,jk) = psol(ji,jk) - zgamm(ji,jk+1) * psol(ji,jk+1)
667          END DO
668       ENDDO
669
670       IF( ln_timing )  CALL timing_stop('sed_mat_dsri')
671
672
673    END SUBROUTINE sed_mat_dsri
674
675
676 END MODULE sedmat
Note: See TracBrowser for help on using the repository browser.