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 @ 15127

Last change on this file since 15127 was 15127, checked in by cetlod, 2 years ago

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