source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmat.F90 @ 4528

Last change on this file since 4528 was 3443, checked in by cetlod, 9 years ago

branch:2012/dev_r3438_LOCEAN15_PISLOB : 1st step of the merge, see ticket #972

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