source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmbc.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

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