1 | MODULE sedmbc |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE sedmbc *** |
---|
4 | !! Sediment : mass balance calculation |
---|
5 | !!===================================================================== |
---|
6 | |
---|
7 | !!---------------------------------------------------------------------- |
---|
8 | !! sed_mbc : |
---|
9 | !!---------------------------------------------------------------------- |
---|
10 | !! * Modules used |
---|
11 | USE sed ! sediment global variable |
---|
12 | USE seddsr |
---|
13 | USE lib_mpp ! distribued memory computing library |
---|
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 | !! $Id$ |
---|
31 | CONTAINS |
---|
32 | |
---|
33 | |
---|
34 | SUBROUTINE sed_mbc( kt ) |
---|
35 | !!---------------------------------------------------------------------- |
---|
36 | !! *** ROUTINE sed_mbc *** |
---|
37 | !! |
---|
38 | !! ** Purpose : computation of total tracer inventories for checking |
---|
39 | !! mass conservation. |
---|
40 | !! |
---|
41 | !! |
---|
42 | !! ** Method : tracer inventories of each reservoir are computed and added |
---|
43 | !! subsequently. |
---|
44 | !! |
---|
45 | !! History : |
---|
46 | !! ! 04-10 (N. Emprin, M. Gehlen ) Original code |
---|
47 | !! ! 06-07 (C. Ethe) Re-organization |
---|
48 | !!---------------------------------------------------------------------- |
---|
49 | |
---|
50 | !! Arguments |
---|
51 | INTEGER, INTENT(in) :: kt ! time step |
---|
52 | |
---|
53 | !! local declarations |
---|
54 | INTEGER :: ji,js, jw, jk |
---|
55 | REAL(wp) :: zinit, zfinal |
---|
56 | REAL(wp) :: zinput, zoutput |
---|
57 | REAL(wp) :: zdsw, zvol |
---|
58 | REAL, DIMENSION(jpsol) :: zsolcp_inv_i, zsolcp_inv_f |
---|
59 | REAL, DIMENSION(jpwat) :: zpwcp_inv_i, zpwcp_inv_f |
---|
60 | REAL(wp) :: zdelta_sil, zdelta_clay |
---|
61 | REAL(wp) :: zdelta_co2, zdelta_fe |
---|
62 | REAL(wp) :: zdelta_po4, zdelta_no3 |
---|
63 | |
---|
64 | !!---------------------------------------------------------------------- |
---|
65 | ! Initilization |
---|
66 | !--------------- |
---|
67 | IF( ln_timing ) CALL timing_start('sed_mbc') |
---|
68 | ! |
---|
69 | IF( kt == nitsed000 ) THEN |
---|
70 | |
---|
71 | DO js = 1, jpsol |
---|
72 | rain_tot (js) = 0. |
---|
73 | fromsed_tot(js) = 0. |
---|
74 | tosed_tot (js) = 0. |
---|
75 | rloss_tot (js) = 0. |
---|
76 | ENDDO |
---|
77 | |
---|
78 | DO jw = 1, jpwat |
---|
79 | diss_in_tot (jw) = 0. |
---|
80 | diss_out_tot(jw) = 0. |
---|
81 | ENDDO |
---|
82 | |
---|
83 | ENDIF |
---|
84 | |
---|
85 | |
---|
86 | ! Calculation of the cumulativ input and output |
---|
87 | ! for mass balance check |
---|
88 | !---------------------------------------------- |
---|
89 | |
---|
90 | ! cumulativ solid |
---|
91 | DO js = 1, jpsol |
---|
92 | DO ji = 1, jpoce |
---|
93 | ! input [mol] |
---|
94 | rain_tot (js) = rain_tot (js) + dtsed * rainrm_dta(ji,js) |
---|
95 | fromsed_tot(js) = fromsed_tot(js) + fromsed(ji,js) / mol_wgt(js) |
---|
96 | ! output [mol] |
---|
97 | tosed_tot (js) = tosed_tot (js) + tosed(ji,js) / mol_wgt(js) |
---|
98 | rloss_tot (js) = rloss_tot (js) + rloss(ji,js) / mol_wgt(js) |
---|
99 | ENDDO |
---|
100 | ENDDO |
---|
101 | |
---|
102 | ! cumulativ dissolved |
---|
103 | DO jw = 1, jpwat |
---|
104 | DO ji = 1, jpoce |
---|
105 | ! input [mol] |
---|
106 | diss_in_tot (jw) = diss_in_tot (jw) + pwcp_dta(ji,jw) * 1.e-3 * dzkbot(ji) |
---|
107 | ! output [mol] |
---|
108 | diss_out_tot(jw) = diss_out_tot(jw) + tokbot(ji,jw) |
---|
109 | ENDDO |
---|
110 | ENDDO |
---|
111 | |
---|
112 | ! Mass balance check |
---|
113 | !--------------------- |
---|
114 | IF( kt == nitsedend ) THEN |
---|
115 | ! initial and final inventories for solid component (mole/dx.dy) in sediment |
---|
116 | zsolcp_inv_i(:) = 0. |
---|
117 | zsolcp_inv_f(:) = 0. |
---|
118 | zpwcp_inv_i (:) = 0. |
---|
119 | zpwcp_inv_f (:) = 0. |
---|
120 | DO js = 1, jpsol |
---|
121 | zdsw = denssol / mol_wgt(js) |
---|
122 | DO jk = 2, jpksed |
---|
123 | DO ji = 1, jpoce |
---|
124 | zvol = vols3d(ji,jk) * zdsw |
---|
125 | zsolcp_inv_i(js) = zsolcp_inv_i(js) + solcp0(ji,jk,js) * zvol |
---|
126 | zsolcp_inv_f(js) = zsolcp_inv_f(js) + solcp (ji,jk,js) * zvol |
---|
127 | ENDDO |
---|
128 | END DO |
---|
129 | ENDDO |
---|
130 | |
---|
131 | ! initial and final inventories for dissolved component (mole/dx.dy) in sediment |
---|
132 | DO jw = 1, jpwat |
---|
133 | DO jk = 2, jpksed |
---|
134 | DO ji = 1, jpoce |
---|
135 | zvol = volw3d(ji,jk) * 1.e-3 |
---|
136 | zpwcp_inv_i(jw) = zpwcp_inv_i(jw) + pwcp0(ji,jk,jw) * zvol |
---|
137 | zpwcp_inv_f(jw) = zpwcp_inv_f(jw) + pwcp (ji,jk,jw) * zvol |
---|
138 | ENDDO |
---|
139 | END DO |
---|
140 | ENDDO |
---|
141 | |
---|
142 | ! mass balance for Silica/opal |
---|
143 | zinit = zsolcp_inv_i(jsopal) + zpwcp_inv_i(jwsil) |
---|
144 | zfinal = zsolcp_inv_f(jsopal) + zpwcp_inv_f(jwsil) |
---|
145 | zinput = rain_tot (jsopal) + diss_in_tot (jwsil) |
---|
146 | zoutput = tosed_tot (jsopal) + rloss_tot (jsopal) + diss_out_tot(jwsil) |
---|
147 | zdelta_sil = ( zfinal + zoutput ) - ( zinit + zinput ) |
---|
148 | |
---|
149 | |
---|
150 | ! mass balance for Clay |
---|
151 | zinit = zsolcp_inv_i(jsclay) |
---|
152 | zfinal = zsolcp_inv_f(jsclay) |
---|
153 | zinput = rain_tot (jsclay) + fromsed_tot(jsclay) |
---|
154 | zoutput = tosed_tot (jsclay) + rloss_tot (jsclay) |
---|
155 | zdelta_clay= ( zfinal + zoutput ) - ( zinit + zinput ) |
---|
156 | |
---|
157 | ! mass balance for carbon ( carbon in POC, CaCo3, DIC ) |
---|
158 | zinit = zsolcp_inv_i(jspoc) + zsolcp_inv_i(jspos) + zsolcp_inv_i(jspor) & |
---|
159 | & + zsolcp_inv_i(jscal) + zpwcp_inv_i(jwdic) |
---|
160 | zfinal = zsolcp_inv_f(jspoc) + zsolcp_inv_f(jspos) + zsolcp_inv_f(jspor) & |
---|
161 | & + zsolcp_inv_f(jscal) + zpwcp_inv_f(jwdic) |
---|
162 | zinput = rain_tot (jspoc) + rain_tot (jspos) + rain_tot (jspor) & |
---|
163 | & + rain_tot (jscal) + diss_in_tot(jwdic) |
---|
164 | zoutput = tosed_tot(jspoc) + tosed_tot(jspos) + tosed_tot(jspor) + tosed_tot(jscal) + diss_out_tot(jwdic) & |
---|
165 | & + rloss_tot(jspoc) + rloss_tot(jspos) + rloss_tot(jspor) + rloss_tot(jscal) |
---|
166 | zdelta_co2 = ( zfinal + zoutput ) - ( zinit + zinput ) |
---|
167 | |
---|
168 | ! mass balance for Sulfur |
---|
169 | zinit = zpwcp_inv_i(jwso4) + zpwcp_inv_i(jwh2s) & |
---|
170 | & + zsolcp_inv_i(jsfes) |
---|
171 | zfinal = zpwcp_inv_f(jwso4) + zpwcp_inv_f(jwh2s) & |
---|
172 | & + zsolcp_inv_f(jsfes) |
---|
173 | zinput = diss_in_tot (jwso4) + diss_in_tot (jwh2s) & |
---|
174 | & + rain_tot (jsfes) |
---|
175 | zoutput = diss_out_tot(jwso4) + diss_out_tot(jwh2s) & |
---|
176 | & + tosed_tot(jsfes) + rloss_tot(jsfes) |
---|
177 | zdelta_no3 = ( zfinal + zoutput ) - ( zinit + zinput ) |
---|
178 | |
---|
179 | ! mass balance for iron |
---|
180 | zinit = zpwcp_inv_i(jwfe2) + zsolcp_inv_i(jsfeo) & |
---|
181 | & + zsolcp_inv_i(jsfes) |
---|
182 | zfinal = zpwcp_inv_f(jwfe2) + zsolcp_inv_f(jsfeo) & |
---|
183 | & + zsolcp_inv_f(jsfes) |
---|
184 | zinput = diss_in_tot (jwfe2) + rain_tot (jsfeo) & |
---|
185 | & + rain_tot (jsfes) |
---|
186 | zoutput = diss_out_tot(jwfe2) + tosed_tot(jsfeo) & |
---|
187 | & + tosed_tot(jsfes) + rloss_tot(jsfes) + rloss_tot(jsfeo) |
---|
188 | zdelta_fe = ( zfinal + zoutput ) - ( zinit + zinput ) |
---|
189 | |
---|
190 | |
---|
191 | END IF |
---|
192 | |
---|
193 | IF( kt == nitsedend) THEN |
---|
194 | |
---|
195 | IF (lwp) THEN |
---|
196 | WRITE(numsed,*) |
---|
197 | WRITE(numsed,*)'================== General mass balance ================== ' |
---|
198 | WRITE(numsed,*)' ' |
---|
199 | WRITE(numsed,*)' ' |
---|
200 | WRITE(numsed,*)' Initial total solid Masses (mole/dx.dy) ' |
---|
201 | WRITE(numsed,*)' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' |
---|
202 | WRITE(numsed,*)' Opal, Clay, POC, POS, POR, CaCO3, FeOH, FeS' |
---|
203 | WRITE(numsed,'(8x,4(1PE10.3,2X))')zsolcp_inv_i(jsopal),zsolcp_inv_i(jsclay),zsolcp_inv_i(jspoc), & |
---|
204 | & zsolcp_inv_i(jspos),zsolcp_inv_i(jspor),zsolcp_inv_i(jscal),zsolcp_inv_i(jsfeo),zsolcp_inv_i(jsfes) |
---|
205 | WRITE(numsed,*)' ' |
---|
206 | WRITE(numsed,*)' Initial total dissolved Masses (mole/dx.dy) ' |
---|
207 | WRITE(numsed,*)' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' |
---|
208 | WRITE(numsed,*)' Si, O2, DIC, Nit, Phos, Fe2+' |
---|
209 | WRITE(numsed,'(5x,5(1PE10.3,2X))') zpwcp_inv_i(jwsil), zpwcp_inv_i(jwoxy), & |
---|
210 | & zpwcp_inv_i(jwdic), zpwcp_inv_i(jwno3), zpwcp_inv_i(jwpo4), zpwcp_inv_i(jwfe2) |
---|
211 | WRITE(numsed,*)' ' |
---|
212 | WRITE(numsed,*)' Solid inputs : Opale, Clay, POC, CaCO3, Fe' |
---|
213 | WRITE(numsed,'(A4,10X,5(1PE10.3,2X))')'Rain : ',rain_tot(jsopal),rain_tot(jsclay),rain_tot(jspoc) + rain_tot(jspos) + rain_tot(jspor),& |
---|
214 | & rain_tot(jscal), rain_tot(jsfeo) |
---|
215 | WRITE(numsed,'(A12,6x,5(1PE10.3,2X))')' From Sed : ',fromsed_tot(jsopal), fromsed_tot(jsclay), & |
---|
216 | & fromsed_tot(jspoc)+fromsed_tot(jspos)+fromsed_tot(jspor), fromsed_tot(jscal), & |
---|
217 | & fromsed_tot(jsfeo) + fromsed_tot(jsfes) |
---|
218 | WRITE(numsed,*)'Diss. inputs : Si, O2, DIC, Nit, Phos, Fe' |
---|
219 | WRITE(numsed,'(A9,1x,6(1PE10.3,2X))')' From Pisc : ', diss_in_tot(jwsil), & |
---|
220 | & diss_in_tot(jwoxy), diss_in_tot(jwdic), diss_in_tot(jwno3), diss_in_tot(jwpo4), diss_in_tot(jwfe2) |
---|
221 | WRITE(numsed,*)' ' |
---|
222 | WRITE(numsed,*)'Solid output : Opale, Clay, POC, CaCO3, Fe' |
---|
223 | WRITE(numsed,'(A6,8x,5(1PE10.3,2X))')'To sed', tosed_tot(jsopal),tosed_tot(jsclay),tosed_tot(jspoc) & |
---|
224 | & +tosed_tot(jspos)+tosed_tot(jspor),tosed_tot(jscal), tosed_tot(jsfeo)+tosed_tot(jsfes) |
---|
225 | WRITE(numsed,'(A5,9x,5(1PE10.3,2X))')'Perdu', rloss_tot(jsopal),rloss_tot(jsclay),rloss_tot(jspoc) & |
---|
226 | & +rloss_tot(jspos)+rloss_tot(jspor),rloss_tot(jscal),rloss_tot(jsfeo)+rloss_tot(jsfes) |
---|
227 | WRITE(numsed,*)'Diss. output : Si, O2, DIC, Nit, Phos, Fe ' |
---|
228 | WRITE(numsed,'(A7,2x,6(1PE10.3,2X))')'To kbot', diss_out_tot(jwsil), & |
---|
229 | & diss_out_tot(jwoxy), diss_out_tot(jwdic), diss_out_tot(jwno3), diss_out_tot(jwpo4), diss_out_tot(jwfe2) |
---|
230 | WRITE(numsed,*)' ' |
---|
231 | WRITE(numsed,*)'Final solid Masses (mole/dx.dy) ' |
---|
232 | WRITE(numsed,*)' Opale, Clay, POC, CaCO3, Fe' |
---|
233 | WRITE(numsed,'(4x,5(1PE10.3,2X))')zsolcp_inv_f(jsopal),zsolcp_inv_f(jsclay),zsolcp_inv_f(jspoc) & |
---|
234 | & +zsolcp_inv_f(jspos)+zsolcp_inv_f(jspor),zsolcp_inv_f(jscal),zsolcp_inv_f(jsfeo)+zsolcp_inv_f(jsfes) |
---|
235 | WRITE(numsed,*)' ' |
---|
236 | WRITE(numsed,*)'Final dissolved Masses (mole/dx.dy) (k=2-11)' |
---|
237 | WRITE(numsed,*)' Si, O2, DIC, Nit, Phos, Fe' |
---|
238 | WRITE(numsed,'(4x,6(1PE10.3,2X))') zpwcp_inv_f(jwsil), zpwcp_inv_f(jwoxy), & |
---|
239 | & zpwcp_inv_f(jwdic), zpwcp_inv_f(jwno3), zpwcp_inv_f(jwpo4), zpwcp_inv_f(jwfe2) |
---|
240 | WRITE(numsed,*)' ' |
---|
241 | WRITE(numsed,*)'Delta : Opale, Clay, C, Fe, S,' |
---|
242 | WRITE(numsed,'(7x,6(1PE11.3,1X))') zdelta_sil / ( zsolcp_inv_i(jsopal) + zpwcp_inv_i(jwsil) ) , & |
---|
243 | & zdelta_clay / ( zsolcp_inv_i(jsclay) ) , & |
---|
244 | & zdelta_co2 / ( zsolcp_inv_i(jspoc) + zsolcp_inv_i(jspos) + zsolcp_inv_i(jspor) & |
---|
245 | & + zsolcp_inv_i(jscal) + zpwcp_inv_i(jwdic) ), & |
---|
246 | & zdelta_fe / ( zpwcp_inv_i(jwfe2) + zsolcp_inv_i(jsfeo) + zsolcp_inv_i(jsfes) ) , & |
---|
247 | & zdelta_no3 / ( zpwcp_inv_i(jwso4) + zpwcp_inv_i(jwh2s) + zsolcp_inv_i(jsfes) ) |
---|
248 | WRITE(numsed,*)'==========================================================================' |
---|
249 | |
---|
250 | ENDIF |
---|
251 | ENDIF |
---|
252 | |
---|
253 | IF( ln_timing ) CALL timing_stop('sed_mbc') |
---|
254 | |
---|
255 | END SUBROUTINE sed_mbc |
---|
256 | |
---|
257 | END MODULE sedmbc |
---|