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/UKMO/NEMO_4.0_mirror/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror/src/TOP/PISCES/SED/sedmat.F90 @ 10888

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

branches/UKMO/NEMO_4.0_mirror : clear SVN keywords

File size: 9.0 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 
17
18   INTERFACE sed_mat
19      MODULE PROCEDURE sed_mat_dsr, sed_mat_btb
20   END INTERFACE
21
22   INTEGER, PARAMETER :: nmax = 30 
23
24
25   !! $Id$
26 CONTAINS
27
28    SUBROUTINE sed_mat_dsr( nvar, ndim, nlev, preac, psms, psol, dtsed_in )
29       !!---------------------------------------------------------------------
30       !!                  ***  ROUTINE sed_mat_dsr  ***
31       !!
32       !! ** Purpose :  solves tridiagonal system of linear equations
33       !!
34       !! ** Method  :
35       !!        1 - computes left hand side of linear system of equations
36       !!            for dissolution reaction
37       !!         For mass balance in kbot+sediment :
38       !!              dz3d  (:,1) = dz(1) = 0.5 cm
39       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )
40       !!              dz(2)       = 0.3 cm
41       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )     
42       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see seddsr.F90 )
43       !!
44       !!         2 - forward/backward substitution.
45       !!
46       !!   History :
47       !!        !  04-10 (N. Emprin, M. Gehlen ) original
48       !!        !  06-04 (C. Ethe)  Module Re-organization
49       !!----------------------------------------------------------------------
50       !! * Arguments
51       INTEGER , INTENT(in) ::  nvar  ! number of variable
52       INTEGER , INTENT(in) ::  ndim  ! number of points
53       INTEGER , INTENT(in) ::  nlev  ! number of sediment levels
54
55       REAL(wp), DIMENSION(ndim,nlev), INTENT(in   ) :: preac  ! reaction rates
56       REAL(wp), DIMENSION(ndim,nlev), INTENT(in   ) :: psms  ! reaction rates
57       REAL(wp), DIMENSION(ndim,nlev), INTENT(inout) :: psol   ! solution ( undersaturation values )
58       REAL(wp), INTENT(in) ::  dtsed_in
59 
60       !---Local declarations
61       INTEGER  ::  ji, jk, jn
62       REAL(wp), DIMENSION(ndim,nlev) :: za, zb, zc, zr
63       REAL(wp), DIMENSION(ndim)      :: zbet
64       REAL(wp), DIMENSION(ndim,nmax) :: zgamm
65
66       REAL(wp) ::  aplus,aminus 
67       REAL(wp) ::  rplus,rminus   
68       REAL(wp) ::  dxplus,dxminus
69
70       !----------------------------------------------------------------------
71
72       IF( ln_timing )  CALL timing_start('sed_mat_dsr')
73
74       ! Computation left hand side of linear system of
75       ! equations for dissolution reaction
76       !---------------------------------------------
77
78
79       jn = nvar
80       ! first sediment level         
81       DO ji = 1, ndim
82          aplus  = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + &
83                        ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2.
84          dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2.
85          rplus  = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 
86
87          za(ji,1) = 0.
88          zb(ji,1) = 1. + rplus
89          zc(ji,1) = -rplus
90       ENDDO
91 
92       DO jk = 2, nlev - 1
93          DO ji = 1, ndim
94             aminus  = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + &
95             &        ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) ) / 2.
96             dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2.
97
98             aplus   = ( ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) + &
99             &        ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2.
100             dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2
101                !
102             rminus  = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus
103             rplus   = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn)   * aplus / dxplus
104                !     
105             za(ji,jk) = -rminus
106             zb(ji,jk) = 1. + rminus + rplus 
107             zc(ji,jk) = -rplus
108          END DO
109       END DO
110
111       DO ji = 1, ndim
112          aminus  = ( ( volw3d(ji,nlev-1) / dz3d(ji,nlev-1)  ) + &
113          &        ( volw3d(ji,nlev)  / dz3d(ji,nlev) ) ) / 2.
114          dxminus = ( dz3d(ji,nlev-1) + dz3d(ji,nlev) ) / 2.
115          rminus  = ( dtsed_in / volw3d(ji,nlev) ) * diff(ji,nlev-1,jn) * aminus / dxminus
116          !
117          za(ji,nlev) = -rminus
118          zb(ji,nlev) = 1. + rminus
119          zc(ji,nlev) = 0.
120       END DO
121
122
123       ! solves tridiagonal system of linear equations
124       ! -----------------------------------------------
125
126       zr  (:,:) = psol(:,:) + psms(:,:)
127       zb  (:,:) = zb(:,:) + preac(:,:)
128       zbet(:  ) = zb(:,1)
129       psol(:,1) = zr(:,1) / zbet(:)
130
131          !
132       DO jk = 2, nlev
133          DO ji = 1, ndim
134             zgamm(ji,jk) =  zc(ji,jk-1) / zbet(ji)
135             zbet(ji)     =  zb(ji,jk) - za(ji,jk) * zgamm(ji,jk)
136             psol(ji,jk)  = ( zr(ji,jk) - za(ji,jk) * psol(ji,jk-1) ) / zbet(ji)
137          END DO
138       ENDDO
139          !
140       DO jk = nlev - 1, 1, -1
141          DO ji = 1,ndim
142             psol(ji,jk) = psol(ji,jk) - zgamm(ji,jk+1) * psol(ji,jk+1)
143          END DO
144       ENDDO
145
146       IF( ln_timing )  CALL timing_stop('sed_mat_dsr')
147
148
149    END SUBROUTINE sed_mat_dsr
150
151    SUBROUTINE sed_mat_btb( nvar, ndim, nlev, psol, dtsed_in )
152       !!---------------------------------------------------------------------
153       !!                  ***  ROUTINE sed_mat_btb  ***
154       !!
155       !! ** Purpose :  solves tridiagonal system of linear equations
156       !!
157       !! ** Method  :
158       !!        1 - computes left hand side of linear system of equations
159       !!            for dissolution reaction
160       !!
161       !!         2 - forward/backward substitution.
162       !!
163       !!   History :
164       !!        !  04-10 (N. Emprin, M. Gehlen ) original
165       !!        !  06-04 (C. Ethe)  Module Re-organization
166       !!----------------------------------------------------------------------
167       !! * Arguments
168       INTEGER , INTENT(in) :: &
169          nvar , &  ! number of variables
170          ndim , &  ! number of points
171          nlev      ! number of sediment levels
172
173      REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(inout) :: &
174          psol      ! solution
175
176      REAL(wp), INTENT(in) :: dtsed_in
177
178       !---Local declarations
179       INTEGER  ::  &
180          ji, jk, jn
181
182       REAL(wp) ::  &
183          aplus,aminus   ,  & 
184          rplus,rminus   ,  &
185          dxplus,dxminus 
186
187       REAL(wp), DIMENSION(nlev)      :: za, zb, zc
188       REAL(wp), DIMENSION(ndim,nlev) :: zr
189       REAL(wp), DIMENSION(nmax)      :: zgamm 
190       REAL(wp) ::  zbet
191
192 
193       !----------------------------------------------------------------------
194
195       ! Computation left hand side of linear system of
196       ! equations for dissolution reaction
197       !---------------------------------------------
198
199
200      IF( ln_timing )  CALL timing_start('sed_mat_btb')
201
202       ! first sediment level         
203      DO ji = 1, ndim
204         aplus  = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2.
205         dxplus = ( dz(2) + dz(3) ) / 2.
206         rplus  = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus 
207
208         za(1) = 0.
209         zb(1) = 1. + rplus
210         zc(1) = -rplus
211
212             
213         DO jk = 2, nlev - 1
214            aminus  = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2.
215            dxminus = ( dz(jk) + dz(jk+1) ) / 2.
216            rminus  = ( dtsed_in / vols(jk+1) ) * db(ji,jk) * aminus / dxminus
217            !
218            aplus   = ( ( vols(jk+1) / dz(jk+1  ) ) + ( vols(jk+2) / dz(jk+2) ) ) / 2.
219            dxplus  = ( dz(jk+1) + dz(jk+2) ) / 2.
220            rplus   = ( dtsed_in / vols(jk+1) ) * db(ji,jk+1) * aplus / dxplus
221            !     
222            za(jk) = -rminus
223            zb(jk) = 1. + rminus + rplus 
224            zc(jk) = -rplus
225         ENDDO
226 
227         aminus  = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2.
228         dxminus = ( dz(nlev) + dz(nlev+1) ) / 2.
229         rminus  = ( dtsed_in / vols(nlev+1) ) * db(ji,nlev) * aminus / dxminus
230         !
231         za(nlev) = -rminus
232         zb(nlev) = 1. + rminus
233         zc(nlev) = 0.
234
235
236         ! solves tridiagonal system of linear equations
237         ! -----------------------------------------------   
238         DO jn = 1, nvar
239         
240            DO jk = 1, nlev
241               zr  (ji,jk)    = psol(ji,jk,jn)
242            END DO
243            zbet          = zb(1)
244            psol(ji,1,jn) = zr(ji,1) / zbet
245            !
246            DO jk = 2, nlev
247               zgamm(jk) =  zc(jk-1) / zbet
248               zbet      =  zb(jk) - za(jk) * zgamm(jk)
249               psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet
250            ENDDO
251            !
252            DO jk = nlev - 1, 1, -1
253                psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn)
254            ENDDO
255
256         ENDDO
257
258       END DO
259       !
260       IF( ln_timing )  CALL timing_stop('sed_mat_btb')
261
262       
263    END SUBROUTINE sed_mat_btb
264
265 END MODULE sedmat
Note: See TracBrowser for help on using the repository browser.