source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedmbc.F90 @ 10345

Last change on this file since 10345 was 10345, checked in by smasson, 2 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10344, see #2133

  • Property svn:keywords set to Id
File size: 11.7 KB
Line 
1MODULE 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$
31CONTAINS
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
257END MODULE sedmbc
Note: See TracBrowser for help on using the repository browser.