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.
p4zsed.F90 in NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90 @ 14276

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

numerous updates to PISCES, PISCES-QUOTA and the sediment module
  • Property svn:keywords set to Id
File size: 30.7 KB
Line 
1MODULE p4zsed
2   !!======================================================================
3   !!                         ***  MODULE p4sed  ***
4   !! TOP :   PISCES Compute loss of biogenic matter in the sediments
5   !!                Compute gain/loss of tracers from dust, rivers and
6   !!                sediments
7   !!                This module is used both by PISCES and PISCES-QUOTA
8   !!======================================================================
9   !! History :   1.0  !  2004-03 (O. Aumont) Original code
10   !!             2.0  !  2007-12 (C. Ethe, G. Madec)  F90
11   !!             3.4  !  2011-06 (C. Ethe) USE of fldread
12   !!             3.5  !  2012-07 (O. Aumont) improvment of river input of nutrients
13   !!-----------------------------------------------------------------------
14   !!   p4z_sed        :  Compute loss of organic matter in the sediments
15   !!                  :  Compute gain/loss of tracers from dust, rivers and
16   !!                     sediments
17   !!-----------------------------------------------------------------------
18   USE oce_trc         !  shared variables between ocean and passive tracers
19   USE trc             !  passive tracers common variables
20   USE sms_pisces      !  PISCES Source Minus Sink variables
21   USE p4zlim          !  Co-limitations of differents nutrients
22   USE p4zsbc          !  External source of nutrients
23   USE p4zint          !  interpolation and computation of various fields
24   USE sed             !  Sediment module
25   USE iom             !  I/O manager
26   USE prtctl_trc      !  print control for debugging
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   p4z_sed 
32   PUBLIC   p4z_sed_alloc
33 
34   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nitrpot    !: Nitrogen fixation
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: sdenit     !: Nitrate reduction in the sediments
36   REAL(wp) :: r1_rday                  !: inverse of rday
37   LOGICAL, SAVE :: lk_sed              !: Explicit sediment module
38
39   !!----------------------------------------------------------------------
40   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
41   !! $Id$
42   !! Software governed by the CeCILL license (see ./LICENSE)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE p4z_sed( kt, knt )
47      !!---------------------------------------------------------------------
48      !!                     ***  ROUTINE p4z_sed  ***
49      !!
50      !! ** Purpose : Compute the loss of biogenic matter in the sediments. This
51      !!              is by no way a real sediment model. The loss is simply
52      !!              computed from metamodels.
53      !!              Loss/gain of tracers are also computed here for
54      !!              dust, rivers, sediments and hydrothermal vents (Fe)
55      !!              N2 fixation is modeled using an implicit approach
56      !!
57      !! ** Method  : - Fluxes with the sediments are mainly modeled using
58      !!                statiscal metamodels.
59      !!---------------------------------------------------------------------
60      !
61      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
62      INTEGER  ::  ji, jj, jk, ikt
63      REAL(wp) ::  zrivalk, zrivsil, zrivno3
64      REAL(wp) ::  zwflux, zlim, zfact, zfactcal
65      REAL(wp) ::  zo2, zno3, zflx, zpdenit, z1pdenit, zolimit
66      REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep
67      REAL(wp) ::  zwstpoc, zwstpon, zwstpop
68      REAL(wp) ::  ztrpo4s, ztrdp, zwdust, zmudia, ztemp
69      REAL(wp) ::  xdiano3, xdianh4
70      !
71      CHARACTER (len=25) :: charout
72      REAL(wp), DIMENSION(jpi,jpj    ) :: zdenit2d, zbureff, zwork
73      REAL(wp), DIMENSION(jpi,jpj    ) :: zwsbio3, zwsbio4
74      REAL(wp), DIMENSION(jpi,jpj    ) :: zsedcal, zsedsi, zsedc
75      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsoufer, zlight, ztrfer
76      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zpdep, zsidep
77      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zironice
78      !!---------------------------------------------------------------------
79      !
80      IF( ln_timing )  CALL timing_start('p4z_sed')
81      !
82      IF( kt == nittrc000 .AND. knt == 1 )   THEN
83          r1_rday  = 1. / rday
84          ! Configuration with an active two-way sediment module
85          IF (ln_sediment .AND. ln_sed_2way) THEN
86             lk_sed = .TRUE.
87          ELSE
88             lk_sed = .FALSE.
89          ENDIF
90      ENDIF
91      !
92      IF( kt == nittrc000 .AND. knt == 1 )   r1_rday  = 1. / rday
93      !
94      ! Allocate temporary workspace
95      ALLOCATE( ztrpo4(jpi,jpj,jpk) )
96      IF( ln_p5z )    ALLOCATE( ztrdop(jpi,jpj,jpk) )
97
98      zdenit2d(:,:) = 0.e0
99      zbureff (:,:) = 0.e0
100      zwork   (:,:) = 0.e0
101      zsedsi  (:,:) = 0.e0
102      zsedcal (:,:) = 0.e0
103      zsedc   (:,:) = 0.e0
104
105      ! Iron input/uptake due to sea ice : Crude parameterization based on
106      ! Lancelot et al. Iron concentration in sea-ice is constant and set
107      ! in the namelist_pisces (icefeinput). ln_ironice is forced to false
108      ! when nn_ice_tr = 1
109      ! ----------------------------------------------------
110      IF( ln_ironice ) THEN 
111         !                                             
112         ALLOCATE( zironice(jpi,jpj) )
113
114         ! Compute the iron flux between sea ice and sea water
115         ! Simple parameterization assuming a fixed constant concentration in
116         ! sea-ice (icefeinput)
117         ! ------------------------------------------------------------------
118         DO jj = 1, jpj
119            DO ji = 1, jpi
120               zdep    = rfact2 / e3t_n(ji,jj,1)
121               zwflux  = fmmflx(ji,jj) / 1000._wp
122               zironice(ji,jj) =  MAX( -0.99 * trb(ji,jj,1,jpfer), -zwflux * icefeinput * zdep )
123            END DO
124         END DO
125         ! Update of the TRA array
126         tra(:,:,1,jpfer) = tra(:,:,1,jpfer) + zironice(:,:) 
127         !
128         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironice" ) )   &
129            &   CALL iom_put( "Ironice", zironice(:,:) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) ) ! iron flux from ice
130         !
131         DEALLOCATE( zironice )
132         !                                             
133      ENDIF
134
135      ! Add the external input of nutrients from dust deposition
136      ! ----------------------------------------------------------
137      IF( ln_dust ) THEN
138         !                                             
139         ALLOCATE( zsidep(jpi,jpj,jpk), zpdep(jpi,jpj,jpk), zirondep(jpi,jpj,jpk) )
140
141         ! Iron, P and Si deposition at the surface
142         ! Iron flux at the surface due to dust deposition. Solubility can be
143         ! be variable if ln_solub is set to true. In that case, solubility
144         ! has to be provided in the specific input file (read in p4zsbc)
145         ! mfrac is the mean iron relative weight content of dust
146         ! ------------------------------------------------------------------
147         IF( ln_solub ) THEN
148            zirondep(:,:,1) = solub(:,:) * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss 
149         ELSE
150            zirondep(:,:,1) = dustsolub  * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss 
151         ENDIF
152
153         ! Si and P flux at the surface due to dust deposition. The content
154         ! and the solubility are hard coded
155         ! ----------------------------------------------------------------
156         ! Si crustal abundance is about 26.9% by mass and 7.5% is soluble
157         ! (see for instance, Tegen and Kohfeld, 2006).
158         zsidep(:,:,1) = 0.269 * 0.075 * dust(:,:) * rfact2 / e3t_n(:,:,1) / 28.1 
159         ! P Crustal abundance is about 0.1% by mass and about 10% of it soluble
160         ! (Paytan and McLaughlin, 2007).
161         zpdep (:,:,1) = 0.1 * 1.e-3 * dust(:,:) * rfact2 / e3t_n(:,:,1) / 31. / po4r 
162
163         ! Iron solubilization of particles in the water column
164         ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/d
165         ! Dust are supposed to sink at wdust sinking speed. 3% of the iron
166         ! in dust is hypothesized to be soluble at a dissolution rate set to
167         ! 1/(250 days). The vertical distribution of iron in dust is computed
168         ! from a steady state assumption. Parameters are very uncertain and
169         ! are estimated from the literature quoted in Raiswell et al. (2011)
170         ! -------------------------------------------------------------------
171         zwdust = 0.03 * rday / ( wdust * 55.85 ) / ( 250. * rday )
172         DO jk = 2, jpkm1
173            zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -gdept_n(:,:,jk) / (250. * wdust) )
174            zpdep   (:,:,jk) = zirondep(:,:,jk) * 1.e-3 / mfrac * 55.85 / 31.
175            zsidep  (:,:,jk) = zirondep(:,:,jk) * 0.269 / mfrac * 55.85 / 28.1
176         END DO
177
178         ! Solubilization of particles in the water column (Si, P, Fe)
179         DO jk = 1, jpkm1
180            tra(:,:,jk,jppo4) = tra(:,:,jk,jppo4) + zpdep   (:,:,jk)
181            tra(:,:,jk,jpfer) = tra(:,:,jk,jpfer) + zirondep(:,:,jk) 
182            tra(:,:,jk,jpsil) = tra(:,:,jk,jpsil) + zsidep  (:,:,jk)
183         ENDDO
184         !
185         IF( lk_iomput ) THEN
186            IF( knt == nrdttrc ) THEN
187                IF( iom_use( "Irondep" ) )   &
188                &  CALL iom_put( "Irondep", zirondep(:,:,1) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron
189                IF( iom_use( "pdust" ) )   &
190                &  CALL iom_put( "pdust"  , dust(:,:) / ( wdust / rday )  * tmask(:,:,1) ) ! dust concentration at surface
191            ENDIF
192         ENDIF
193         DEALLOCATE( zsidep, zpdep, zirondep )
194         !                                             
195      ENDIF
196     
197      ! Add the external input of nutrients from river
198      ! ----------------------------------------------
199      IF( ln_river ) THEN
200         DO jj = 1, jpj
201            DO ji = 1, jpi
202               DO jk = 1, nk_rnf(ji,jj)
203                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) +  rivdip(ji,jj) * rfact2
204                  tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) +  rivdin(ji,jj) * rfact2
205                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) +  rivdic(ji,jj) * 5.e-5 * rfact2
206                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) +  rivdsi(ji,jj) * rfact2
207                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +  rivdic(ji,jj) * rfact2
208                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) +  ( rivalk(ji,jj) - rno3 * rivdin(ji,jj) ) * rfact2
209                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) +  rivdoc(ji,jj) * rfact2
210               ENDDO
211            ENDDO
212         ENDDO
213
214         ! When prognostic ligands are activated, ligands are supplied
215         ! to the ocean by rivers. We assume that the amount of ligands
216         ! is equal to that of iron (iron is completely complexed)
217         ! ------------------------------------------------------------
218         IF (ln_ligand) THEN
219            DO jj = 1, jpj
220               DO ji = 1, jpi
221                  DO jk = 1, nk_rnf(ji,jj)
222                     tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) +  rivdic(ji,jj) * 5.e-5 * rfact2
223                  ENDDO
224               ENDDO
225            ENDDO
226         ENDIF
227         ! PISCES-QUOTA part
228         IF( ln_p5z ) THEN
229            DO jj = 1, jpj
230               DO ji = 1, jpi
231                  DO jk = 1, nk_rnf(ji,jj)
232                     tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + rivdop(ji,jj) * rfact2
233                     tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + rivdon(ji,jj) * rfact2
234                  ENDDO
235               ENDDO
236            ENDDO
237         ENDIF
238      ENDIF
239     
240      ! Add the external input of nutrients from nitrogen deposition
241      ! ----------------------------------------------------------
242      IF( ln_ndepo ) THEN
243         tra(:,:,1,jpno3) = tra(:,:,1,jpno3) + nitdep(:,:) * rfact2
244         tra(:,:,1,jptal) = tra(:,:,1,jptal) - rno3 * nitdep(:,:) * rfact2
245      ENDIF
246
247      ! Add the external input of iron from hydrothermal vents
248      ! Please refer to Tagliabue et al. (2010) for more information
249      ! ------------------------------------------------------------
250      IF( ln_hydrofe ) THEN
251            tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2
252         IF( ln_ligand ) THEN
253            tra(:,:,:,jplgw) = tra(:,:,:,jplgw) + ( hydrofe(:,:,:) * lgw_rath ) * rfact2
254         ENDIF
255         !
256         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "HYDR" ) )   &
257            &   CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input
258      ENDIF
259
260      ! OA: Warning, the following part is necessary to avoid CFL problems
261      ! above the sediments. Vertical sinking speed is limited using the
262      ! typical CFL criterion
263      ! --------------------------------------------------------------------
264      DO jj = 1, jpj
265         DO ji = 1, jpi
266            ikt  = mbkt(ji,jj)
267            zdep = e3t_n(ji,jj,ikt) / xstep
268            zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) )
269            zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) )
270         END DO
271      END DO
272      !
273      ! No sediment module activated
274      IF( .NOT.lk_sed ) THEN
275!
276         ! Add the external input of iron from sediment mobilization
277         ! ------------------------------------------------------
278         IF( ln_ironsed ) THEN
279            tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2
280            !
281            IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironsed" ) )   &
282               &   CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! iron inputs from sediments
283         ENDIF
284
285         ! Computation of the sediment denitrification proportion: The metamodel
286         ! from Middleburg (2006) is used
287         ! Computation of the fraction of organic matter that is permanently
288         ! buried from Dunne's model (2007)
289         ! -------------------------------------------------------
290         DO jj = 1, jpj
291            DO ji = 1, jpi
292              IF( tmask(ji,jj,1) == 1 ) THEN
293                 ikt = mbkt(ji,jj)
294                 zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
295                   &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4
296                 zflx  = LOG10( MAX( 1E-3, zflx ) )
297                 zo2   = LOG10( MAX( 10. , trb(ji,jj,ikt,jpoxy) * 1E6 ) )
298                 zno3  = LOG10( MAX( 1.  , trb(ji,jj,ikt,jpno3) * 1E6 * rno3 ) )
299                 zdep  = LOG10( gdepw_n(ji,jj,ikt+1) )
300                 zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    &
301                   &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2
302                 zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) )
303                   !
304                 zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
305                   &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6
306                 zbureff(ji,jj) = 0.013 + 0.13 * zflx**2 / ( 7.0 + zflx )**2
307              ENDIF
308            END DO
309         END DO 
310         !
311      ENDIF
312
313      ! Fraction of dSi that is dissolved in the sediments. This fraction is 
314      ! set to a constant value in p4zsbc
315      ! --------------------------------------------------------------------
316      IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac
317
318      ! Loss of bSi and CaCO3 to the sediments
319      DO jj = 1, jpj
320         DO ji = 1, jpi
321            ikt  = mbkt(ji,jj)
322            zdep = xstep / e3t_n(ji,jj,ikt) 
323            zwsc = zwsbio4(ji,jj) * zdep
324            zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc
325            zcaloss = trb(ji,jj,ikt,jpcal) * zwsc
326            !
327            tra(ji,jj,ikt,jpgsi) = tra(ji,jj,ikt,jpgsi) - zsiloss
328            tra(ji,jj,ikt,jpcal) = tra(ji,jj,ikt,jpcal) - zcaloss
329         END DO
330      END DO
331      !
332      IF( .NOT.lk_sed ) THEN
333         ! Dissolution of CaCO3 and bSi in the sediments. This is
334         ! instantaneous since here sediments are not explicitly
335         ! modeled. The amount of CaCO3 that dissolves in the sediments
336         ! is computed using a metamodel constructed from Archer (1996)
337         ! A minimum set to sedcalfrac is preserved. This value is defined
338         ! in p4zsbc
339         ! ---------------------------------------------------------------
340         DO jj = 1, jpj
341            DO ji = 1, jpi
342               ikt  = mbkt(ji,jj)
343               zdep = xstep / e3t_n(ji,jj,ikt) 
344               zwsc = zwsbio4(ji,jj) * zdep
345               zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc
346               zcaloss = trb(ji,jj,ikt,jpcal) * zwsc
347               tra(ji,jj,ikt,jpsil) = tra(ji,jj,ikt,jpsil) + zsiloss * zrivsil 
348               !
349               zfactcal = MAX(-0.1, MIN( excess(ji,jj,ikt), 0.2 ) )
350               zfactcal = 0.3 + 0.7 * MIN( 1., (0.1 + zfactcal) / ( 0.5 - zfactcal ) )
351               zrivalk  = sedcalfrac * zfactcal
352               tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0
353               tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk
354               zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss * e3t_n(ji,jj,ikt) 
355               zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss * e3t_n(ji,jj,ikt) 
356            END DO
357         END DO
358      ENDIF
359      !
360      ! Loss of particulate organic carbon and Fe to the sediments
361      DO jj = 1, jpj
362         DO ji = 1, jpi
363            ikt  = mbkt(ji,jj)
364            zdep = xstep / e3t_n(ji,jj,ikt) 
365            zws4 = zwsbio4(ji,jj) * zdep
366            zws3 = zwsbio3(ji,jj) * zdep
367            tra(ji,jj,ikt,jpgoc) = tra(ji,jj,ikt,jpgoc) - trb(ji,jj,ikt,jpgoc) * zws4 
368            tra(ji,jj,ikt,jppoc) = tra(ji,jj,ikt,jppoc) - trb(ji,jj,ikt,jppoc) * zws3
369            tra(ji,jj,ikt,jpbfe) = tra(ji,jj,ikt,jpbfe) - trb(ji,jj,ikt,jpbfe) * zws4
370            tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3
371         END DO
372      END DO
373      !
374      ! Loss of particulate organic N and P to the sediments (p5z)
375      IF( ln_p5z ) THEN
376         DO jj = 1, jpj
377            DO ji = 1, jpi
378               ikt  = mbkt(ji,jj)
379               zdep = xstep / e3t_n(ji,jj,ikt) 
380               zws4 = zwsbio4(ji,jj) * zdep
381               zws3 = zwsbio3(ji,jj) * zdep
382               tra(ji,jj,ikt,jpgon) = tra(ji,jj,ikt,jpgon) - trb(ji,jj,ikt,jpgon) * zws4
383               tra(ji,jj,ikt,jppon) = tra(ji,jj,ikt,jppon) - trb(ji,jj,ikt,jppon) * zws3
384               tra(ji,jj,ikt,jpgop) = tra(ji,jj,ikt,jpgop) - trb(ji,jj,ikt,jpgop) * zws4
385               tra(ji,jj,ikt,jppop) = tra(ji,jj,ikt,jppop) - trb(ji,jj,ikt,jppop) * zws3
386            END DO
387         END DO
388      ENDIF
389
390      IF( .NOT.lk_sed ) THEN
391         ! Degradation of organic matter in the sediments. The metamodel of
392         ! Middleburg (2006) is used here to mimic the diagenetic reactions.
393         ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after
394         ! denitrification in the sediments. Not very clever, but simpliest option.
395         ! ------------------------------------------------------------------------
396         DO jj = 1, jpj
397            DO ji = 1, jpi
398               ikt  = mbkt(ji,jj)
399               zdep = xstep / e3t_n(ji,jj,ikt) 
400               zws4 = zwsbio4(ji,jj) * zdep
401               zws3 = zwsbio3(ji,jj) * zdep
402               ! Fraction that is permanently buried in the sediments
403               zrivno3 = 1. - zbureff(ji,jj)
404               zwstpoc = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3
405               ! Denitrification in the sediments
406               zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 )
407               ! Fraction that is not denitrified
408               z1pdenit = zwstpoc * zrivno3 - zpdenit
409               ! Oxic remineralization of organic matter in the sediments
410               zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) )
411               ! The fraction that cannot be denitrified nor oxidized by O2
412               ! is released back to the water column as DOC
413               tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit
414               ! Update of the tracers concentrations
415               tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit
416               tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit
417               tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) - rdenit * zpdenit
418               tra(ji,jj,ikt,jpoxy) = tra(ji,jj,ikt,jpoxy) - zolimit * o2ut
419               tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * zpdenit )
420               tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit 
421               sdenit(ji,jj) = rdenit * zpdenit * e3t_n(ji,jj,ikt)
422               zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t_n(ji,jj,ikt)
423               ! PISCES-QUOTA (p5z)
424               IF( ln_p5z ) THEN
425                  zwstpop              = trb(ji,jj,ikt,jpgop) * zws4 + trb(ji,jj,ikt,jppop) * zws3
426                  zwstpon              = trb(ji,jj,ikt,jpgon) * zws4 + trb(ji,jj,ikt,jppon) * zws3
427                  tra(ji,jj,ikt,jpdon) = tra(ji,jj,ikt,jpdon) + ( z1pdenit - zolimit ) * zwstpon / (zwstpoc + rtrn)
428                  tra(ji,jj,ikt,jpdop) = tra(ji,jj,ikt,jpdop) + ( z1pdenit - zolimit ) * zwstpop / (zwstpoc + rtrn)
429               ENDIF
430            END DO
431         END DO
432       ENDIF
433
434
435      ! Nitrogen fixation process : light limitation of diazotrophy
436      ! Small source of iron from particulate inorganic iron (photochemistry)
437      ! This is a purely adhoc param.
438      !----------------------------------------------------------------------
439      DO jk = 1, jpkm1
440         zlight (:,:,jk) =  ( 1.- EXP( -etot_ndcy(:,:,jk) / diazolight ) ) * ( 1. - fr_i(:,:) ) 
441         zsoufer(:,:,jk) = zlight(:,:,jk) * 1E-10 / ( 1E-10 + biron(:,:,jk) )
442      ENDDO
443
444      ! Diazotrophy (nitrogen fixation) is modeled according to an empirical
445      ! formulation. This is described in Aumont et al. (2015). Limitation
446      ! by P and Fe is computed. Inhibition by high N concentrations is imposed.
447      ! Diazotrophy sensitivity to temperature is parameterized as in
448      ! Ye et al. (2012) 
449      ! ------------------------------------------------------------------------
450      IF( ln_p4z ) THEN
451         ! PISCES part
452         DO jk = 1, jpkm1
453            DO jj = 1, jpj
454               DO ji = 1, jpi
455                  ! Potential nitrogen fixation dependant on temperature
456                  ztemp = tsn(ji,jj,jk,jp_tem)
457                  zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) / rno3
458                  ! Nitrogen fixation is inhibited when enough NO3 and/or NH4
459                  zlim = ( 1.- xnanonh4(ji,jj,jk) - xnanono3(ji,jj,jk) )
460                  zfact = zlim * rfact2
461                  ! Nitrogen fixation limitation by PO4 and Fe
462                  ztrfer(ji,jj,jk) = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
463                  ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( concnnh4 + trb(ji,jj,jk,jppo4) )
464                  ztrdp = ztrpo4(ji,jj,jk)
465                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer(ji,jj,jk), ztrdp ) * zlight(ji,jj,jk)
466               END DO
467            END DO
468         END DO
469      ELSE       ! p5z
470         ! PISCES-QUOTA part
471         DO jk = 1, jpkm1
472            DO jj = 1, jpj
473               DO ji = 1, jpi
474                  ! Potential nitrogen fixation dependant on temperature
475                  ztemp = tsn(ji,jj,jk,jp_tem)
476                  zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
477                  ! Nitrogen fixation is inhibited when enough NO3 and/or NH4
478                  xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) )
479                  xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
480                  zlim = ( 1.- xdiano3 - xdianh4 )
481                  zfact = zlim * rfact2
482                  ! Nitrogen fixation limitation by PO4/DOP and Fe
483                  ztrfer(ji,jj,jk) = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
484                  ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) )
485                  ztrdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( 1E-6 + trb(ji,jj,jk,jpdop) ) * (1. - ztrpo4(ji,jj,jk))
486                  ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk)
487                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer(ji,jj,jk), ztrdp ) * zlight(ji,jj,jk)
488               END DO
489            END DO
490         END DO
491      ENDIF
492
493      ! Update of the TRA arrays due to nitrogen fixation
494      ! -------------------------------------------------
495      IF( ln_p4z ) THEN
496         ! PISCES part
497         DO jk = 1, jpkm1
498            DO jj = 1, jpj
499               DO ji = 1, jpi
500                  zfact = nitrpot(ji,jj,jk) * nitrfix
501                  ! 1/3 of the diazotrophs growth is supposed to be excreted
502                  ! as NH4. 1/3 as DOC and the rest is routed to POC/GOC as
503                  ! a result of mortality by predation. Completely adhoc param
504                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0
505                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0
506                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zfact * 2.0 / 3.0
507                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zfact * 1.0 / 3.0
508                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zfact * 1.0 / 3.0 * 2.0 / 3.0
509                  tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0
510                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
511                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C at max
512                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0
513                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 2.0 / 3.0
514                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 1.0 / 3.0
515                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.005 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
516                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trb(ji,jj,jk,jppo4) ) &
517                  &                     * 0.001 * trb(ji,jj,jk,jpdoc) * xstep
518              END DO
519            END DO
520         END DO
521      ELSE    ! p5z
522         ! PISCES-QUOTA part
523         DO jk = 1, jpkm1
524            DO jj = 1, jpj
525               DO ji = 1, jpi
526                  zfact = nitrpot(ji,jj,jk) * nitrfix
527                  ! 1/3 of the diazotrophs growth is supposed to be excreted
528                  ! as NH4. 1/3 as DOC and the rest is routed POC and GOC as
529                  ! a result of mortality by predation. Completely adhoc param
530                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0
531                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0
532                  ! N/P ratio of diazotrophs is supposed to be 46
533                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) &
534                  &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
535                  tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zfact * 1.0 / 3.0
536                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zfact * 1.0 / 3.0
537                  tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + 16.0 / 46.0 * zfact / 3.0  &
538                  &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   &
539                  &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
540                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zfact * 1.0 / 3.0 * 2.0 / 3.0
541                  tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + zfact * 1.0 / 3.0 * 2.0 /3.0
542                  tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 2.0 /3.0
543                  tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0
544                  tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) + zfact * 1.0 / 3.0 * 1.0 /3.0
545                  tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0
546                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
547                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C at max
548                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 
549                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 2.0 / 3.0
550                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 1.0 / 3.0
551                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
552              END DO
553            END DO
554         END DO
555         !
556      ENDIF
557
558      IF( lk_iomput ) THEN
559         IF( knt == nrdttrc ) THEN
560            zfact = 1.e+3 * rfact2r !  conversion from molC/l/kt  to molN/m3/s
561            IF( iom_use("Nfix"   ) ) CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation
562            IF( iom_use("INTNFIX") ) THEN   ! nitrogen fixation rate in ocean ( vertically integrated )
563               zwork(:,:) = 0.
564               DO jk = 1, jpkm1
565                 zwork(:,:) = zwork(:,:) + nitrpot(:,:,jk) * nitrfix * rno3 * zfact * e3t_n(:,:,jk) * tmask(:,:,jk)
566               ENDDO
567               CALL iom_put( "INTNFIX" , zwork ) 
568            ENDIF
569            IF( iom_use("SedCal" ) ) CALL iom_put( "SedCal", zsedcal(:,:) * zfact ) ! Permanent burial of CaCO3 in sediments
570            IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * zfact ) ! Permanent burial of bSi in sediments
571            IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * zfact ) ! Permanent burial of OC in sediments
572!            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 ) ! Denitrification in the sediments
573            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", zdenit2d(:,:) ) ! Denitrification in the sediments
574         ENDIF
575      ENDIF
576      !
577      IF(ln_ctl) THEN  ! print mean trends (USEd for debugging)
578         WRITE(charout, fmt="('sed ')")
579         CALL prt_ctl_trc_info(charout)
580         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
581      ENDIF
582      !
583      IF( ln_p5z )    DEALLOCATE( ztrpo4, ztrdop )
584      !
585      IF( ln_timing )  CALL timing_stop('p4z_sed')
586      !
587   END SUBROUTINE p4z_sed
588
589
590   INTEGER FUNCTION p4z_sed_alloc()
591      !!----------------------------------------------------------------------
592      !!                     ***  ROUTINE p4z_sed_alloc  ***
593      !!----------------------------------------------------------------------
594      ALLOCATE( nitrpot(jpi,jpj,jpk), sdenit(jpi,jpj), STAT=p4z_sed_alloc )
595      !
596      IF( p4z_sed_alloc /= 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_alloc: failed to allocate arrays' )
597      !
598   END FUNCTION p4z_sed_alloc
599
600   !!======================================================================
601END MODULE p4zsed
Note: See TracBrowser for help on using the repository browser.