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.
seddsr.F90 in NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddsr.F90

Last change on this file was 15297, checked in by aumont, 3 years ago

major update of the sediment module
  • Property svn:keywords set to Id
File size: 11.3 KB
Line 
1MODULE seddsr
2   !!======================================================================
3   !!              ***  MODULE  seddsr  ***
4   !!    Sediment : dissolution and reaction in pore water related
5   !!    related to organic matter
6   !!=====================================================================
7   !! * Modules used
8   USE sed     ! sediment global variable
9   USE sed_oce
10   USE sedini
11   USE lib_mpp         ! distribued memory computing library
12   USE lib_fortran
13
14   IMPLICIT NONE
15   PRIVATE
16
17   PUBLIC sed_dsr
18
19   !! * Module variables
20
21   REAL(wp), DIMENSION(jpsol), PUBLIC      :: dens_mol_wgt  ! molecular density
22
23   !! $Id$
24CONTAINS
25   
26   SUBROUTINE sed_dsr( accmask )
27      !!----------------------------------------------------------------------
28      !!                   ***  ROUTINE sed_dsr  ***
29      !!
30      !!  ** Purpose :  computes pore water dissolution and reaction
31      !!
32      !!  ** Methode :  Computation of the redox reactions in sediment.
33      !!                The main redox reactions are solved in sed_dsr whereas
34      !!                the secondary reactions are solved in sed_dsr_redoxb.
35      !!                A strand spliting approach is being used here (see
36      !!                sed_dsr_redoxb for more information).
37      !!
38      !!   History :
39      !!        !  98-08 (E. Maier-Reimer, Christoph Heinze )  Original code
40      !!        !  04-10 (N. Emprin, M. Gehlen ) f90
41      !!        !  06-04 (C. Ethe)  Re-organization
42      !!        !  19-08 (O. Aumont) Debugging and improvement of the model.
43      !!                             The original method is replaced by a
44      !!                              Strand splitting method which deals
45      !!                              well with stiff reactions.
46      !!----------------------------------------------------------------------
47      !! Arguments
48      ! --- local variables
49      INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask
50      INTEGER :: ji, jk, js, jw, jn   ! dummy looop indices
51
52      REAL(wp), DIMENSION(jpoce,jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo    ! undersaturation ; indice jpwatp1 is for calcite   
53      REAL(wp)  ::  zreasat
54      !!
55      !!----------------------------------------------------------------------
56
57      IF( ln_timing )  CALL timing_start('sed_dsr')
58!
59     ! Initializations
60     !----------------------
61     
62      zlimo2 (:,:)    = 0.    ;   zlimno3(:,:)  = 0.
63      zlimso4(:,:)    = 0.
64 
65      !----------------------------------------------------------
66      ! 5.  Beginning of solid reaction
67      !---------------------------------------------------------
68     
69      ! Computes the SMS of the species which do not affect the remin
70      ! processes (Alkalinity, PO4, NH4, and the release of Fe from
71      ! biogenic Fe
72      ! In the following, all loops start from jpk = 2 as reactions
73      ! only occur in the sediment
74      ! --------------------------------------------------------------
75      DO jk = 2, jpksed
76         DO ji = 1, jpoce
77            IF (accmask(ji) == 0 ) THEN
78               zreasat = rearatpom(ji,jk)
79               ! For alkalinity
80               pwcpa(ji,jk,jwalk)  = pwcpa(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r )
81               ! For Phosphate (in mol/l)
82               pwcpa(ji,jk,jwpo4)  = pwcpa(ji,jk,jwpo4) + zreasat * spo4r
83               ! For Ammonium
84               pwcpa(ji,jk,jwnh4)  = pwcpa(ji,jk,jwnh4) + zreasat * srno3 * radssol(jk,jwnh4)
85               ! For iron (in mol/l)
86               pwcpa(ji,jk,jwfe2)  = pwcpa(ji,jk,jwfe2) + fecratio(ji) * zreasat * radssol(jk,jwfe2)
87            ENDIF
88         END DO
89      ENDDO
90
91      ! Computes SMS of oxygen
92      DO jk = 2, jpksed
93         DO ji = 1, jpoce
94            IF (accmask(ji) == 0 ) THEN
95               zlimo2(ji,jk) = pwcp(ji,jk,jwoxy) / ( pwcp(ji,jk,jwoxy) + xksedo2 )
96               zlimo2(ji,jk) = MAX( 0., zlimo2(ji,jk) )
97               zreasat = rearatpom(ji,jk) * zlimo2(ji,jk)
98               ! Acid Silicic
99               pwcpa(ji,jk,jwoxy)  = pwcpa(ji,jk,jwoxy) - so2ut * zreasat
100            ENDIF
101         END DO
102      ENDDO
103
104      !--------------------------------------------------------------------
105      ! Denitrification
106      !--------------------------------------------------------------------
107      DO jk = 2, jpksed
108         DO ji = 1, jpoce
109            IF (accmask(ji) == 0 ) THEN
110               zlimno3(ji,jk) = ( 1.0 - zlimo2(ji,jk) ) * pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 ) 
111               zlimno3(ji,jk) = MAX(0., zlimno3(ji,jk) )
112               zreasat = rearatpom(ji,jk) * zlimno3(ji,jk) * srDnit
113               ! For nitrates
114               pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) - zreasat
115               ! For alkalinity
116               pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat
117            ENDIF
118         END DO
119      ENDDO
120
121
122      !--------------------------------------------------------------------
123      ! Begining POC iron reduction
124      !--------------------------------------------------------------------
125
126      DO jk = 2, jpksed
127         DO ji = 1, jpoce
128            IF (accmask(ji) == 0 ) THEN
129               zlimfeo(ji,jk) = ( 1.0 - zlimno3(ji,jk) - zlimo2(ji,jk) ) * solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo )
130               zlimfeo(ji,jk) = MAX(0., zlimfeo(ji,jk) )
131               zreasat = 4.0 * rearatpom(ji,jk) * zlimfeo(ji,jk)
132               ! For FEOH
133               solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - zreasat / volc(ji,jk,jsfeo)
134               ! For PO4
135               pwcpa(ji,jk,jwpo4)  = pwcpa(ji,jk,jwpo4) + zreasat * redfep
136               ! For alkalinity
137               pwcpa(ji,jk,jwalk)  = pwcpa(ji,jk,jwalk) + 2.0 * zreasat
138               ! Iron
139               pwcpa(ji,jk,jwfe2)  = pwcpa(ji,jk,jwfe2) + zreasat * radssol(jk,jwfe2)
140            ENDIF
141         END DO
142      ENDDO
143
144      !--------------------------------------------------------------------
145      ! Begining POC denitrification and NO3- diffusion
146      !--------------------------------------------------------------------
147
148      DO jk = 2, jpksed
149         DO ji = 1, jpoce
150            IF (accmask(ji) == 0 ) THEN
151               zlimso4(ji,jk) = ( 1.0 - zlimno3(ji,jk) - zlimo2(ji,jk) - zlimfeo(ji,jk) ) * pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 )
152               zlimso4(ji,jk) = MAX(0., zlimso4(ji,jk) )
153               zreasat = 0.5 * rearatpom(ji,jk) * zlimso4(ji,jk)
154               ! For sulfur
155               pwcpa(ji,jk,jwso4)  =  pwcpa(ji,jk,jwso4) - zreasat 
156               pwcpa(ji,jk,jwh2s)  = pwcpa(ji,jk,jwh2s) + zreasat
157               ! For alkalinity
158               pwcpa(ji,jk,jwalk)  = pwcpa(ji,jk,jwalk) + 2.0 * zreasat
159            ENDIF
160         END DO
161      ENDDO
162
163      ! Secondary redox reactions
164      ! -------------------------
165
166      call sed_dsr_redoxb( accmask )
167
168      IF( ln_timing )  CALL timing_stop('sed_dsr')
169!     
170   END SUBROUTINE sed_dsr
171
172   SUBROUTINE sed_dsr_redoxb( accmask )
173      !!----------------------------------------------------------------------
174      !!                   ***  ROUTINE sed_dsr_redox  ***
175      !!
176      !!  ** Purpose :  computes secondary redox reactions
177      !!
178      !!   History :
179      !!        !  18-08 (O. Aumont)  Original code
180      !!----------------------------------------------------------------------
181      !! Arguments
182      ! --- local variables
183      INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask
184      INTEGER   ::  ji, jni, jnj, jk   ! dummy looop indices
185
186      REAL(wp)  ::  zalpha, zexcess, zh2seq, zsedfer, zreasat
187      !!
188      !!----------------------------------------------------------------------
189
190      IF( ln_timing )  CALL timing_start('sed_dsr_redoxb')
191
192      DO jk = 2, jpksed
193         DO ji = 1, jpoce
194            IF (accmask(ji) == 0 ) THEN
195               zalpha = ( pwcp(ji,jk,jwfe2) - pwcp(ji,jk,jwlgw) ) * 1E9
196               zsedfer = ( zalpha + SQRT( zalpha**2 + 1.E-10 ) ) /2.0 * 1E-9
197
198               ! First pass of the scheme. At the end, it is 1st order
199               ! -----------------------------------------------------
200               ! Fe (both adsorbed and solute) + O2
201               zreasat = reac_fe2 * zsedfer / radssol(jk,jwfe2) * pwcp(ji,jk,jwoxy)
202               pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2)
203               pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 0.25 * zreasat 
204               pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) - redfep * zreasat 
205               pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat
206               solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) + zreasat / volc(ji,jk,jsfeo)
207
208               ! H2S + O2
209               zreasat = reac_h2s * pwcp(ji,jk,jwoxy) * pwcp(ji,jk,jwh2s)
210               pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat
211               pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat
212               pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat
213               pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat
214
215               ! NH4 + O2
216               zreasat = reac_nh4 * pwcp(ji,jk,jwnh4) / radssol(jk,jwnh4) * pwcp(ji,jk,jwoxy)
217               pwcpa(ji,jk,jwnh4) = pwcpa(ji,jk,jwnh4) - zreasat * radssol(jk,jwnh4) 
218               pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 
219               pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) + zreasat
220               pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat
221
222               ! FeS - O2
223               zreasat = reac_feso * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) * pwcp(ji,jk,jwoxy) 
224               solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) - zreasat / volc(ji,jk,jsfes)
225               pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat
226               pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * radssol(jk,jwfe2)
227               pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat
228
229               ! FEOH + H2S
230               zreasat = reac_feh2s * solcp(ji,jk,jsfeo) * volc(ji,jk,jsfeo) * pwcp(ji,jk,jwh2s) 
231               solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - 8.0 * zreasat / volc(ji,jk,jsfeo)
232               pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat
233               pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + 8.0 * zreasat * radssol(jk,jwfe2) 
234               pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 14.0 * zreasat 
235               pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat
236               pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + 8.0 * redfep * zreasat
237
238               ! Fe + H2S
239               zh2seq     = MAX(rtrn, 2.1E-3 * hipor(ji,jk))
240               zexcess = pwcp(ji,jk,jwh2s) * zsedfer / zh2seq - 1.0
241               IF ( zexcess >= 0.0 ) THEN
242                  zreasat = reac_fesp * zexcess
243                  pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2)
244                  solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes)
245                  pwcpa(ji,jk,jwh2s)  = pwcpa(ji,jk,jwh2s) - zreasat
246               ELSE
247                  zreasat = reac_fesd * zexcess * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes)
248                  pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2)
249                  solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes)
250                  pwcpa(ji,jk,jwh2s)  = pwcpa(ji,jk,jwh2s) - zreasat
251               ENDIF
252               ! For alkalinity
253               pwcpa(ji,jk,jwalk)  = pwcpa(ji,jk,jwalk) - zreasat * 2.0
254            ENDIF
255         END DO
256     END DO
257
258      IF( ln_timing )  CALL timing_stop('sed_dsr_redoxb')
259
260  END SUBROUTINE sed_dsr_redoxb
261
262END MODULE seddsr
Note: See TracBrowser for help on using the repository browser.