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 branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/TOP_SRC/PISCES/SED – NEMO

source: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmat.F90 @ 5965

Last change on this file since 5965 was 5965, checked in by timgraham, 8 years ago

Upgraded branch to r5518 of trunk (v3.6 stable revision)

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