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

Last change on this file since 15075 was 15075, checked in by aumont, 3 years ago

major update of the sediment module
  • Property svn:keywords set to Id
File size: 16.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
21   INTEGER, PARAMETER :: nmax = 60 
22
23
24   !! $Id$
25 CONTAINS
26
27    SUBROUTINE sed_mat_dsr( ji, nvar, dtsed_in )
28       !!---------------------------------------------------------------------
29       !!                  ***  ROUTINE sed_mat_dsr  ***
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) ::  ji, nvar  ! number of variable
51
52       REAL(wp), INTENT(in) ::  dtsed_in
53 
54       !---Local declarations
55       INTEGER  ::  jk, jn
56       REAL(wp), DIMENSION(jpksed) :: za, zb, zc
57
58       REAL(wp) ::  aplus,aminus 
59       REAL(wp) ::  rplus,rminus   
60       REAL(wp) ::  dxplus,dxminus
61
62       !----------------------------------------------------------------------
63
64       IF( ln_timing )  CALL timing_start('sed_mat_dsr')
65
66       ! Computation left hand side of linear system of
67       ! equations for dissolution reaction
68       !---------------------------------------------
69
70
71       jn = nvar
72       ! first sediment level         
73       aplus  = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + &
74                        ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2.
75       dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2.
76       rplus  = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 
77
78       za(1) = 0.
79       zb(1) = rplus
80       zc(1) = -rplus
81 
82       DO jk = 2, jpksed - 1
83          aminus  = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + &
84          &        ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) ) / 2.
85          dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2.
86
87          aplus   = ( ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) + &
88          &        ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2.
89          dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2
90          !
91          rminus  = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus
92          rplus   = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn)   * aplus / dxplus
93          !     
94          za(jk) = -rminus
95          zb(jk) = rminus + rplus 
96          zc(jk) = -rplus
97       END DO
98
99       aminus  = ( ( volw3d(ji,jpksed-1) / dz3d(ji,jpksed-1)  ) + &
100       &        ( volw3d(ji,jpksed)  / dz3d(ji,jpksed) ) ) / 2.
101       dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2.
102       rminus  = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus / dxminus
103       !
104       za(jpksed) = -rminus
105       zb(jpksed) = rminus
106       zc(jpksed) = 0.
107
108       ! solves tridiagonal system of linear equations
109       ! -----------------------------------------------
110
111       pwcpa(ji,1,jn) = pwcpa(ji,1,jn) - ( zc(1) * pwcp(ji,2,jn) + zb(1) * pwcp(ji,1,jn) ) 
112       DO jk = 2, jpksed - 1
113          pwcpa(ji,jk,jn) =  pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn)    &
114          &                  + zb(jk) * pwcp(ji,jk,jn) )
115       ENDDO
116       pwcpa(ji,jpksed,jn) = pwcpa(ji,jpksed,jn) - ( za(jk) * pwcp(ji,jk-1,jn) + zb(jk) * pwcp(ji,jk,jn) )
117
118       IF( ln_timing )  CALL timing_stop('sed_mat_dsr')
119
120    END SUBROUTINE sed_mat_dsr
121
122    SUBROUTINE sed_mat_dsrjac( ji, nvar, dtsed_in )
123       !!---------------------------------------------------------------------
124       !!                  ***  ROUTINE sed_mat_dsrjac  ***
125       !!
126       !! ** Purpose :  solves tridiagonal system of linear equations
127       !!
128       !! ** Method  :
129       !!        1 - computes left hand side of linear system of equations
130       !!            for dissolution reaction
131       !!         For mass balance in kbot+sediment :
132       !!              dz3d  (:,1) = dz(1) = 0.5 cm
133       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )
134       !!              dz(2)       = 0.3 cm
135       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )     
136       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see
137       !seddsr.F90 )
138       !!
139       !!         2 - forward/backward substitution.
140       !!
141       !!   History :
142       !!        !  04-10 (N. Emprin, M. Gehlen ) original
143       !!        !  06-04 (C. Ethe)  Module Re-organization
144       !!----------------------------------------------------------------------
145       !! * Arguments
146       INTEGER , INTENT(in) ::  ji, nvar  ! number of variable
147
148       REAL(wp), INTENT(in) ::  dtsed_in
149
150       !---Local declarations
151       INTEGER  ::  jk, jn, jnn
152       REAL(wp), DIMENSION(jpksed) :: za, zb, zc
153
154       REAL(wp) ::  aplus,aminus
155       REAL(wp) ::  rplus,rminus
156       REAL(wp) ::  dxplus,dxminus
157
158       !----------------------------------------------------------------------
159
160       IF( ln_timing )  CALL timing_start('sed_mat_dsrjac')
161
162       ! Computation left hand side of linear system of
163       ! equations for dissolution reaction
164       !---------------------------------------------
165
166
167       jn = nvar
168       ! first sediment level         
169       aplus  = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + &
170                        ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2.
171       dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2.
172       rplus  = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus
173
174       za(1) = 0.
175       zb(1) = rplus
176       zc(1) = -rplus
177
178       DO jk = 2, jpksed - 1
179          aminus  = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + &
180          &        ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) ) / 2.
181          dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2.
182
183          aplus   = ( ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) + &
184          &        ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2.
185          dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2
186          !
187          rminus  = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus
188          rplus   = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn)   * aplus / dxplus
189          !     
190          za(jk) = -rminus
191          zb(jk) = rminus + rplus
192          zc(jk) = -rplus
193       END DO
194
195       aminus  = ( ( volw3d(ji,jpksed-1) / dz3d(ji,jpksed-1)  ) + &
196       &        ( volw3d(ji,jpksed)  / dz3d(ji,jpksed) ) ) / 2.
197       dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2.
198       rminus  = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus / dxminus
199       !
200       za(jpksed) = -rminus
201       zb(jpksed) = rminus
202       zc(jpksed) = 0.
203
204       ! solves tridiagonal system of linear equations
205
206       IF (isvode(jn) > 0) THEN
207          jnn = isvode(jn)
208          Jacobian(ji, jnn, jnn) = Jacobian(ji,jnn,jnn) - zb(1)
209          Jacobian(ji, jnn, jpvode + jnn) = Jacobian(ji, jnn, jpvode + jnn) -zc(1)
210          DO jk = 2, jpksed - 1
211             Jacobian(ji, (jk-1) * jpvode + jnn, (jk-2) * jpvode + jnn) = Jacobian(ji, (jk-1) * jpvode + jnn, (jk-2) * jpvode + jnn) - za(jk)
212             Jacobian(ji, (jk-1) * jpvode + jnn, (jk-1) * jpvode + jnn) = Jacobian(ji, (jk-1) * jpvode + jnn, (jk-1) * jpvode + jnn) - zb(jk)
213             Jacobian(ji, (jk-1) * jpvode + jnn, (jk) * jpvode + jnn)   = Jacobian(ji, (jk-1) * jpvode + jnn, (jk) * jpvode + jnn) - zc(jk)
214          END DO
215          Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-2) * jpvode + jnn) = Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-2) * jpvode + jnn) - za(jpksed)
216          Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-1) * jpvode + jnn) = Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-1) * jpvode + jnn) - zb(jpksed)
217       ENDIF
218
219       IF( ln_timing )  CALL timing_stop('sed_mat_dsrjac')
220
221    END SUBROUTINE sed_mat_dsrjac
222
223
224 
225
226    SUBROUTINE sed_mat_btb( nlev, nvar, psol, preac, dtsed_in )
227       !!---------------------------------------------------------------------
228       !!                  ***  ROUTINE sed_mat_btb  ***
229       !!
230       !! ** Purpose :  solves tridiagonal system of linear equations
231       !!
232       !! ** Method  :
233       !!        1 - computes left hand side of linear system of equations
234       !!            for dissolution reaction
235       !!
236       !!         2 - forward/backward substitution.
237       !!
238       !!   History :
239       !!        !  04-10 (N. Emprin, M. Gehlen ) original
240       !!        !  06-04 (C. Ethe)  Module Re-organization
241       !!----------------------------------------------------------------------
242       !! * Arguments
243       INTEGER , INTENT(in) :: &
244          nlev, nvar      ! number of sediment levels
245
246      REAL(wp), DIMENSION(jpoce,nlev,nvar), INTENT(inout) :: &
247          psol, preac
248
249      REAL(wp), INTENT(in) :: dtsed_in
250
251       !---Local declarations
252       INTEGER  ::  &
253          ji, jk, jn
254
255       REAL(wp) ::  &
256          aplus,aminus   ,  &
257          rplus,rminus   ,  &
258          dxplus,dxminus
259
260       REAL(wp), DIMENSION(nlev)      :: za, zb, zc
261       REAL(wp), DIMENSION(jpoce,nlev) :: zr
262       REAL(wp), DIMENSION(nmax)      :: zgamm
263       REAL(wp) ::  zbet
264
265       !----------------------------------------------------------------------
266
267       ! Computation left hand side of linear system of
268       ! equations for dissolution reaction
269       !---------------------------------------------
270
271
272      IF( ln_timing )  CALL timing_start('sed_mat_btb')
273
274       ! first sediment level         
275      DO ji = 1, jpoce
276         aplus  = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2.
277         dxplus = ( dz(2) + dz(3) ) / 2.
278         rplus  = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus
279
280         za(1) = 0.
281         zb(1) = 1. + rplus 
282         zc(1) = -rplus
283
284         DO jk = 2, nlev - 1
285            aminus  = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2.
286            dxminus = ( dz(jk) + dz(jk+1) ) / 2.
287            rminus  = ( dtsed_in / vols(jk+1) ) * db(ji,jk) * aminus / dxminus
288            !
289            aplus   = ( ( vols(jk+1) / dz(jk+1  ) ) + ( vols(jk+2) / dz(jk+2) ) ) / 2.
290            dxplus  = ( dz(jk+1) + dz(jk+2) ) / 2.
291            rplus   = ( dtsed_in / vols(jk+1) ) * db(ji,jk+1) * aplus / dxplus
292            !     
293            za(jk) = -rminus
294            zb(jk) = 1. + rminus + rplus
295            zc(jk) = -rplus
296
297         ENDDO
298
299         aminus  = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2.
300         dxminus = ( dz(nlev) + dz(nlev+1) ) / 2.
301         rminus  = ( dtsed_in / vols(nlev+1) ) * db(ji,nlev) * aminus / dxminus
302         !
303         za(nlev) = -rminus
304         zb(nlev) = 1. + rminus
305         zc(nlev) = 0.
306
307         ! solves tridiagonal system of linear equations
308         ! -----------------------------------------------   
309         DO jn = 1, nvar
310            DO jk = 1, nlev
311               zr(ji,jk) = psol(ji,jk,jn)
312            END DO
313            zbet          = zb(1) - preac(ji,1,jn) * dtsed_in
314            psol(ji,1,jn) = zr(ji,1) / zbet
315            !
316            DO jk = 2, nlev
317               zgamm(jk) =  zc(jk-1) / zbet
318               zbet      =  zb(jk) - preac(ji,jk,jn) * dtsed_in - za(jk) * zgamm(jk)
319               psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet
320            ENDDO
321            !
322            DO jk = nlev - 1, 1, -1
323               psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn)
324            ENDDO
325         END DO
326      END DO
327      !
328      IF( ln_timing )  CALL timing_stop('sed_mat_btb')
329
330       
331    END SUBROUTINE sed_mat_btb
332
333    SUBROUTINE sed_mat_dsri( nvar, preac, psms, dtsed_in )
334       !!---------------------------------------------------------------------
335       !!                  ***  ROUTINE sed_mat_dsr  ***
336       !!
337       !! ** Purpose :  solves tridiagonal system of linear equations
338       !!
339       !! ** Method  :
340       !!        1 - computes left hand side of linear system of equations
341       !!            for dissolution reaction
342       !!         For mass balance in kbot+sediment :
343       !!              dz3d  (:,1) = dz(1) = 0.5 cm
344       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )
345       !!              dz(2)       = 0.3 cm
346       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )     
347       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see
348       !seddsr.F90 )
349       !!
350       !!         2 - forward/backward substitution.
351       !!
352       !!   History :
353       !!        !  04-10 (N. Emprin, M. Gehlen ) original
354       !!        !  06-04 (C. Ethe)  Module Re-organization
355       !!----------------------------------------------------------------------
356       !! * Arguments
357       INTEGER , INTENT(in) ::  nvar  ! number of variable
358
359       REAL(wp), DIMENSION(jpoce,jpksed), INTENT(in   ) :: preac  ! reaction rates
360       REAL(wp), DIMENSION(jpoce,jpksed), INTENT(in   ) :: psms  ! reaction rates
361       REAL(wp), INTENT(in) ::  dtsed_in
362
363       !---Local declarations
364       INTEGER  ::  ji, jk, jn
365       REAL(wp), DIMENSION(jpoce,jpksed) :: za, zb, zc, zr
366       REAL(wp), DIMENSION(jpoce)      :: zbet
367       REAL(wp), DIMENSION(jpoce,nmax) :: zgamm
368
369       REAL(wp) ::  aplus,aminus
370       REAL(wp) ::  rplus,rminus
371       REAL(wp) ::  dxplus,dxminus
372
373       !----------------------------------------------------------------------
374
375       IF( ln_timing )  CALL timing_start('sed_mat_dsri')
376
377       ! Computation left hand side of linear system of
378       ! equations for dissolution reaction
379       !---------------------------------------------
380
381
382       jn = nvar
383       ! first sediment level         
384       DO ji = 1, jpoce
385          aplus  = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + &
386                        ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2.
387          dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2.
388          rplus  = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus
389
390          za(ji,1) = 0.
391          zb(ji,1) = 1. + rplus
392          zc(ji,1) = -rplus
393       ENDDO
394
395       DO jk = 2, jpksed - 1
396          DO ji = 1, jpoce
397             aminus  = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + &
398             &        ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) ) / 2.
399             dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2.
400
401             aplus   = ( ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) + &
402             &        ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2.
403             dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2
404                !
405             rminus  = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus
406             rplus   = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn)   * aplus /dxplus
407                !     
408             za(ji,jk) = -rminus
409             zb(ji,jk) = 1. + rminus + rplus
410             zc(ji,jk) = -rplus
411          END DO
412       END DO
413
414       DO ji = 1, jpoce
415          aminus  = ( ( volw3d(ji,jpksed-1) / dz3d(ji,jpksed-1)  ) + &
416          &        ( volw3d(ji,jpksed)  / dz3d(ji,jpksed) ) ) / 2.
417          dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2.
418          rminus  = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus/ dxminus
419          !
420          za(ji,jpksed) = -rminus
421          zb(ji,jpksed) = 1. + rminus
422          zc(ji,jpksed) = 0.
423       END DO
424
425
426       ! solves tridiagonal system of linear equations
427       ! -----------------------------------------------
428
429       zr  (:,:) = pwcp(:,:,jn) + psms(:,:)
430       zb  (:,:) = zb(:,:) - preac(:,:)
431       zbet(:  ) = zb(:,1)
432       pwcp(:,1,jn) = zr(:,1) / zbet(:)
433
434          !
435       DO jk = 2, jpksed
436          DO ji = 1, jpoce
437             zgamm(ji,jk) =  zc(ji,jk-1) / zbet(ji)
438             zbet(ji)     =  zb(ji,jk) - za(ji,jk) * zgamm(ji,jk)
439             pwcp(ji,jk,jn)  = ( zr(ji,jk) - za(ji,jk) * pwcp(ji,jk-1,jn) ) / zbet(ji)
440          END DO
441       ENDDO
442          !
443       DO jk = jpksed - 1, 1, -1
444          DO ji = 1, jpoce
445             pwcp(ji,jk,jn) = pwcp(ji,jk,jn) - zgamm(ji,jk+1) * pwcp(ji,jk+1,jn)
446          END DO
447       ENDDO
448
449       IF( ln_timing )  CALL timing_stop('sed_mat_dsri')
450
451
452    END SUBROUTINE sed_mat_dsri
453
454
455 END MODULE sedmat
Note: See TracBrowser for help on using the repository browser.