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.
sedmbc.F90 in trunk/NEMO/TOP_SRC/SED – NEMO

source: trunk/NEMO/TOP_SRC/SED/sedmbc.F90 @ 1179

Last change on this file since 1179 was 1179, checked in by cetlod, 16 years ago

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