source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmbc.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: 14.2 KB
Line 
1MODULE sedmbc
2#if defined key_sed
3   !!======================================================================
4   !!              ***  MODULE  sedmbc  ***
5   !! Sediment : mass balance calculation
6   !!=====================================================================
7
8   !!----------------------------------------------------------------------
9   !!   sed_mbc    :
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE sed     ! sediment global variable
13   USE seddsr
14
15   IMPLICIT NONE
16   PRIVATE
17
18   !! *  Routine accessibility
19   PUBLIC sed_mbc
20
21   !! * Module variables
22   REAL(wp), DIMENSION(jpsol) :: rain_tot      ! total input rain
23   REAL(wp), DIMENSION(jpsol) :: fromsed_tot   ! tota input from sediment
24   REAL(wp), DIMENSION(jpsol) :: tosed_tot     ! total output from sediment
25   REAL(wp), DIMENSION(jpsol) :: rloss_tot     ! total rain loss
26
27   REAL(wp), DIMENSION(jpwat) :: diss_in_tot   ! total input in pore water
28   REAL(wp), DIMENSION(jpwat) :: diss_out_tot  ! total output from pore water
29
30   REAL(wp)  :: cons_tot_o2                   ! cumulative o2 consomation
31   REAL(wp)  :: sour_tot_no3                  ! cumulative no3 source
32   REAL(wp)  :: cons_tot_no3                  ! cumulative no3 consomation
33   REAL(wp)  :: sour_tot_c13                  ! cumulative o2 source
34
35   REAL(wp)  :: src13p 
36   REAL(wp)  :: src13ca 
37
38CONTAINS
39
40
41   SUBROUTINE sed_mbc( kt )
42      !!----------------------------------------------------------------------
43      !!                   ***  ROUTINE sed_mbc  ***
44      !!
45      !! ** Purpose :  computation of total tracer inventories for checking
46      !!               mass conservation.
47      !!
48      !!
49      !! ** Method   : tracer inventories of each reservoir are computed and added
50      !!               subsequently.
51      !!
52      !!   History :
53      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
54      !!        !  06-07  (C. Ethe)  Re-organization
55      !!----------------------------------------------------------------------
56
57      !! Arguments
58      INTEGER, INTENT(in) :: kt     ! time step
59
60      !! local declarations
61      INTEGER  :: ji,js, jw, jk
62      REAL(wp) :: zinit, zfinal 
63      REAL(wp) :: zinput, zoutput
64      REAL(wp) :: zdsw, zvol
65      REAL, DIMENSION(jpsol) :: zsolcp_inv_i, zsolcp_inv_f
66      REAL, DIMENSION(jpwat) :: zpwcp_inv_i, zpwcp_inv_f
67      REAL(wp) ::  zdelta_sil, zdelta_clay
68      REAL(wp) ::  zdelta_co2, zdelta_oxy
69      REAL(wp) ::  zdelta_po4, zdelta_no3
70      REAL(wp) ::  zdelta_c13, zdelta_c13b
71
72      !!----------------------------------------------------------------------
73      ! Initilization
74      !---------------
75      IF( kt == nitsed000 ) THEN
76         cons_tot_o2  = 0.
77         sour_tot_no3 = 0.
78         cons_tot_no3 = 0.
79         sour_tot_c13 = 0.
80
81         DO js = 1, jpsol
82            rain_tot   (js) = 0.
83            fromsed_tot(js) = 0.
84            tosed_tot  (js) = 0.
85            rloss_tot  (js) = 0.
86         ENDDO
87
88         DO jw = 1, jpwat
89            diss_in_tot (jw) = 0.
90            diss_out_tot(jw) = 0.
91         ENDDO
92
93         src13p  = rc13P  * pdb
94         src13ca = rc13Ca * pdb
95      ENDIF
96
97
98      ! Calculation of the cumulativ input and output
99      ! for mass balance check
100      !----------------------------------------------
101
102      ! cumulativ solid
103      DO js = 1, jpsol
104         DO ji = 1, jpoce
105            ! input [mol]
106            rain_tot   (js) = rain_tot   (js) + dtsed * rainrm_dta(ji,js)
107            fromsed_tot(js) = fromsed_tot(js) + fromsed(ji,js)
108            ! output [mol]
109            tosed_tot  (js) = tosed_tot (js) + tosed(ji,js)
110            rloss_tot  (js) = rloss_tot (js) + rloss(ji,js)
111         ENDDO
112      ENDDO
113
114      ! cumulativ dissolved
115      DO jw = 1, jpwat
116         DO ji = 1, jpoce
117            ! input [mol]
118            diss_in_tot (jw) = diss_in_tot (jw) + pwcp_dta(ji,jw) * 1.e-3 * dzkbot(ji)
119            ! output [mol]
120            diss_out_tot(jw) = diss_out_tot(jw) + tokbot(ji,jw)
121         ENDDO
122      ENDDO
123
124      ! cumulativ o2 and no3 consomation
125      DO ji = 1, jpoce
126         cons_tot_o2  = cons_tot_o2  + cons_o2 (ji)
127         sour_tot_no3 = sour_tot_no3 + sour_no3(ji)
128         cons_tot_no3 = cons_tot_no3 + cons_no3(ji)
129         sour_tot_c13 = sour_tot_c13 + sour_c13(ji)
130      ENDDO
131     
132
133      ! Mass balance check
134      !---------------------
135      IF( kt == nitsedend ) THEN
136         ! initial and final inventories for solid component (mole/dx.dy) in sediment
137         zsolcp_inv_i(:) = 0.
138         zsolcp_inv_f(:) = 0.
139         zpwcp_inv_i (:) = 0.       
140         zpwcp_inv_f (:) = 0.       
141         DO js = 1, jpsol
142            zdsw = dens / mol_wgt(js)
143            DO jk = 2, jpksed
144               DO ji = 1, jpoce
145                  zvol = vols3d(ji,jk) * zdsw
146                  zsolcp_inv_i(js) = zsolcp_inv_i(js) + solcp0(ji,jk,js) * zvol 
147                  zsolcp_inv_f(js) = zsolcp_inv_f(js) + solcp (ji,jk,js) * zvol
148               ENDDO
149            END DO
150         ENDDO
151
152         ! initial and final inventories for dissolved component (mole/dx.dy) in sediment
153         DO jw = 1, jpwat
154            DO jk = 2, jpksed
155               DO ji = 1, jpoce 
156                  zvol = volw3d(ji,jk) * 1.e-3
157                  zpwcp_inv_i(jw) = zpwcp_inv_i(jw) + pwcp0(ji,jk,jw) * zvol
158                  zpwcp_inv_f(jw) = zpwcp_inv_f(jw) + pwcp (ji,jk,jw) * zvol
159               ENDDO
160            END DO
161         ENDDO
162
163         ! mass balance for Silica/opal
164         zinit      = zsolcp_inv_i(jsopal) + zpwcp_inv_i(jwsil)
165         zfinal     = zsolcp_inv_f(jsopal) + zpwcp_inv_f(jwsil)
166         zinput     = rain_tot    (jsopal) + diss_in_tot (jwsil)
167         zoutput    = tosed_tot   (jsopal) + rloss_tot  (jsopal) + diss_out_tot(jwsil)
168         zdelta_sil = ( zfinal + zoutput ) - ( zinit + zinput )
169
170
171         ! mass balance for Clay
172         zinit      = zsolcp_inv_i(jsclay) 
173         zfinal     = zsolcp_inv_f(jsclay) 
174         zinput     = rain_tot   (jsclay) + fromsed_tot(jsclay) 
175         zoutput    = tosed_tot  (jsclay) + rloss_tot  (jsclay) 
176         zdelta_clay= ( zfinal + zoutput ) - ( zinit + zinput )
177
178         ! mass balance for carbon ( carbon in POC, CaCo3, DIC )
179         zinit      = zsolcp_inv_i(jspoc) + zsolcp_inv_i(jscal) + zpwcp_inv_i(jwdic)
180         zfinal     = zsolcp_inv_f(jspoc) + zsolcp_inv_f(jscal) + zpwcp_inv_f(jwdic)
181         zinput     = rain_tot   (jspoc) + rain_tot   (jscal) + diss_in_tot(jwdic)
182         zoutput    = tosed_tot  (jspoc) + tosed_tot  (jscal) + diss_out_tot(jwdic) &
183            &       + rloss_tot  (jspoc) + rloss_tot  (jscal) 
184         zdelta_co2 = ( zfinal + zoutput ) - ( zinit + zinput )
185
186         ! mass balance for oxygen : O2 is in POC remineralization
187         zinit      = zpwcp_inv_i(jwoxy)
188         zfinal     = zpwcp_inv_f(jwoxy)
189         zinput     = diss_in_tot(jwoxy)
190         zoutput    = diss_out_tot(jwoxy) + cons_tot_o2
191         zdelta_oxy = ( zfinal + zoutput ) - ( zinit + zinput )
192
193         ! mass balance for phosphate ( PO4 in POC and dissolved phosphates )
194         zinit      = zsolcp_inv_i(jspoc) * spo4r + zpwcp_inv_i(jwpo4)
195         zfinal     = zsolcp_inv_f(jspoc) * spo4r + zpwcp_inv_f(jwpo4)
196         zinput     = rain_tot   (jspoc) * spo4r + diss_in_tot(jwpo4)
197         zoutput    = tosed_tot  (jspoc) * spo4r + diss_out_tot(jwpo4) &
198            &       + rloss_tot  (jspoc) * spo4r
199         zdelta_po4 = ( zfinal + zoutput ) - ( zinit + zinput )
200
201
202         ! mass balance for Nitrate
203         zinit      = zpwcp_inv_i  (jwno3)
204         zfinal     = zpwcp_inv_f  (jwno3)
205         zinput     = diss_in_tot (jwno3) + sour_tot_no3
206         zoutput    = diss_out_tot(jwno3) + cons_tot_no3
207         zdelta_no3 = ( zfinal + zoutput ) - ( zinit + zinput )
208
209         ! mass balance for DIC13
210         zinit      =  zpwcp_inv_i(jwc13)   &
211            &        + src13p * zsolcp_inv_i(jspoc) + src13Ca * zsolcp_inv_i(jscal) 
212         zfinal     =  zpwcp_inv_f(jwc13)   &
213            &        + src13p * zsolcp_inv_f(jspoc) + src13Ca * zsolcp_inv_f(jscal)
214         zinput     =  diss_in_tot (jwc13)  &
215            &        + src13p * rain_tot(jspoc) + src13Ca * rain_tot(jscal)
216         zoutput    =  diss_out_tot(jwc13)  &
217            &        + src13p * tosed_tot(jspoc) + src13Ca * tosed_tot(jscal) &   
218            &        + src13p * rloss_tot(jspoc) + src13Ca * rloss_tot(jscal)
219         zdelta_c13 = ( zfinal + zoutput ) - ( zinit + zinput )
220
221         ! other mass balance for DIC13
222         zinit      = zpwcp_inv_i  (jwc13)
223         zfinal     = zpwcp_inv_f  (jwc13)
224         zinput     = diss_in_tot (jwc13) + sour_tot_c13
225         zoutput    = diss_out_tot(jwc13)
226         zdelta_c13b= ( zfinal + zoutput ) - ( zinit + zinput )   
227
228      END IF
229
230      IF( kt == nitsedend) THEN
231
232         WRITE(numsed,*)
233         WRITE(numsed,*)'==================    General mass balance   ==================  '
234         WRITE(numsed,*)' '
235         WRITE(numsed,*)' '
236         WRITE(numsed,*)' Initial total solid Masses (mole/dx.dy) (k=2-11) '
237         WRITE(numsed,*)' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
238         WRITE(numsed,*)'    Opale,      Clay,       POC,        CaCO3,         C13'
239         WRITE(numsed,'(4x,5(1PE10.3,2X))')zsolcp_inv_i(jsopal),zsolcp_inv_i(jsclay),zsolcp_inv_i(jspoc), &
240            & zsolcp_inv_i(jscal),( src13P * zsolcp_inv_i(jspoc) + src13Ca * zsolcp_inv_i(jscal) )
241         WRITE(numsed,*)' '
242         WRITE(numsed,*)' Initial total dissolved Masses (mole/dx.dy) (k=2-11) '
243         WRITE(numsed,*)' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
244         WRITE(numsed,*)'    Si,         O2,         DIC,        Nit         Phos,       DIC13'
245         WRITE(numsed,'(4x,6(1PE10.3,2X))') zpwcp_inv_i(jwsil), zpwcp_inv_i(jwoxy), &
246            & zpwcp_inv_i(jwdic), zpwcp_inv_i(jwno3), zpwcp_inv_i(jwpo4), zpwcp_inv_i(jwc13)
247         WRITE(numsed,*)' '
248         WRITE(numsed,*)'  Solid inputs :  Opale,      Clay,       POC,        CaCO3,         C13'
249         WRITE(numsed,'(A4,10X,5(1PE10.3,2X))')'Rain : ',rain_tot(jsopal),rain_tot(jsclay),rain_tot(jspoc),&
250            & rain_tot(jscal),( src13P * rain_tot(jspoc) + src13Ca * rain_tot(jscal) )
251         WRITE(numsed,'(A12,6x,4(1PE10.3,2X))')' From Sed : ',fromsed_tot(jsopal), fromsed_tot(jsclay), &
252            & fromsed_tot(jspoc), fromsed_tot(jscal)
253         WRITE(numsed,*)'Diss. inputs : Si,    O2,         DIC,         Nit,       Phos,       DIC13'
254         WRITE(numsed,'(A9,1x,6(1PE10.3,2X))')' From Pisc : ', diss_in_tot(jwsil), &
255            & diss_in_tot(jwoxy), diss_in_tot(jwdic), diss_in_tot(jwno3), diss_in_tot(jwpo4), &
256            & diss_in_tot(jwc13)
257         WRITE(numsed,*)' '
258         WRITE(numsed,*)'Solid output : Opale,      Clay,       POC,        CaCO3,          C13'
259         WRITE(numsed,'(A6,8x,5(1PE10.3,2X))')'To sed', tosed_tot(jsopal),tosed_tot(jsclay),tosed_tot(jspoc),&
260            & tosed_tot(jscal),( src13P * tosed_tot(jspoc) + src13Ca * tosed_tot(jscal) )
261         WRITE(numsed,'(A5,9x,5(1PE10.3,2X))')'Perdu', rloss_tot(jsopal),rloss_tot(jsclay),rloss_tot(jspoc),&
262            & rloss_tot(jscal),( src13P * rloss_tot(jspoc) + src13Ca * rloss_tot(jscal) )
263         WRITE(numsed,*)'Diss. output : Si,     O2,        DIC,          Nit,       Phos,       DIC13 ' 
264         WRITE(numsed,'(A7,2x,6(1PE10.3,2X))')'To kbot', diss_out_tot(jwsil), &
265            & diss_out_tot(jwoxy), diss_out_tot(jwdic), diss_out_tot(jwno3), diss_out_tot(jwpo4), &
266            & diss_out_tot(jwc13)
267         WRITE(numsed,*)' '
268         WRITE(numsed,*)' Total consomation in POC remineralization [mol]:  O2,         NO3'
269         WRITE(numsed,'(51x,2(1PE10.3,2X))') cons_tot_o2,cons_tot_no3
270         WRITE(numsed,*)' '
271         WRITE(numsed,*)'Final solid  Masses (mole/dx.dy) (k=2-11)'
272         WRITE(numsed,*)'    Opale,      Clay,       POC,        CaCO3,          C13'
273         WRITE(numsed,'(4x,5(1PE10.3,2X))')zsolcp_inv_f(jsopal),zsolcp_inv_f(jsclay),zsolcp_inv_f(jspoc), &
274            & zsolcp_inv_f(jscal),( src13P * zsolcp_inv_f(jspoc) + src13Ca * zsolcp_inv_f(jscal) )
275         WRITE(numsed,*)' '
276         WRITE(numsed,*)'Final dissolved  Masses (mole/dx.dy) (k=2-11)'
277         WRITE(numsed,*)'    Si,        O2,         DIC,        Nit,        Phos,    DIC13'
278         WRITE(numsed,'(4x,6(1PE10.3,2X))') zpwcp_inv_f(jwsil), zpwcp_inv_f(jwoxy), &
279            & zpwcp_inv_f(jwdic), zpwcp_inv_f(jwno3), zpwcp_inv_f(jwpo4), zpwcp_inv_f(jwc13)
280         WRITE(numsed,*)' '     
281         WRITE(numsed,*)'Delta : Opale,      Clay,       C,          O,          N,          P,        C13'
282         WRITE(numsed,'(7x,7(1PE11.3,1X))') zdelta_sil, zdelta_clay, zdelta_co2, zdelta_oxy, zdelta_no3,&
283            &                          zdelta_po4, zdelta_c13
284         WRITE(numsed,*)' ' 
285         WRITE(numsed,*)'deltaC13bis : ',zdelta_c13b     
286
287         WRITE(numsed,*)'=========================================================================='
288         WRITE(numsed,*)' Composition of final sediment for point jpoce'
289         WRITE(numsed,*)' ========================================='
290         WRITE(numsed,*)'Opale,      Clay,       POC,        CaCo3,      hipor,      pH,         co3por'
291         DO jk = 1,jpksed
292            WRITE(numsed,'(4(F8.4,4X),3(1PE10.3,2X))') solcp(jpoce,jk,jsopal)*100.,solcp(jpoce,jk,jsclay)*100.,&
293               &       solcp(jpoce,jk,jspoc)*100.,solcp(jpoce,jk,jscal)*100.,&
294               &       hipor(jpoce,jk),-LOG10(hipor(jpoce,jk)/densSW(jpoce)),co3por(jpoce,jk)
295         ENDDO
296         WRITE(numsed,*)'Silicic A.,  Oxygen,     DIC,        Nitrats,    Phosphats,  Alkal.,     DIC13'
297         DO jk = 1, jpksed
298            WRITE(numsed,'(8(1PE10.3,2X))')pwcp(jpoce,jk,jwsil),pwcp(jpoce,jk,jwoxy),&
299               & pwcp(jpoce,jk,jwdic),pwcp(jpoce,jk,jwno3),pwcp(jpoce,jk,jwpo4),pwcp(jpoce,jk,jwalk),pwcp(jpoce,jk,jwc13)
300         ENDDO
301         WRITE(numsed,*)'densSW at the end : ',densSW(jpoce)
302         WRITE(numsed,*)'=========================================================================='
303
304      ENDIF
305 
306   END SUBROUTINE sed_mbc
307
308
309#else
310   !!======================================================================
311   !! MODULE sedmbc :   Dummy module
312   !!======================================================================
313CONTAINS
314   SUBROUTINE sed_mbc( kt )         ! Empty routine
315      INTEGER, INTENT(in) :: kt
316      WRITE(*,*) 'sed_mbc: You should not have seen this print! error?', kt
317   END SUBROUTINE sed_mbc
318#endif
319END MODULE sedmbc
Note: See TracBrowser for help on using the repository browser.