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

Last change on this file since 13233 was 13233, checked in by aumont, 4 years ago

update of the PISCES comments

  • Property svn:keywords set to Id
File size: 29.6 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) ::  ztrfer, 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
76      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zpdep
77      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zsidep, 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), 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         zsidep(:,:)   = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 28.1 
157         zpdep (:,:,1) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 31. / po4r 
158
159         ! Iron solubilization of particles in the water column
160         ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/d
161         ! Dust are supposed to sink at wdust sinking speed. 3% of the iron
162         ! in dust is hypothesized to be soluble at a dissolution rate set to
163         ! 1/(250 days). The vertical distribution of iron in dust is computed
164         ! from a steady state assumption. Parameters are very uncertain and
165         ! are estimated from the literature quoted in Raiswell et al. (2011)
166         ! -------------------------------------------------------------------
167         zwdust = 0.03 * rday / ( wdust * 55.85 ) / ( 250. * rday )
168         DO jk = 2, jpkm1
169            zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -gdept_n(:,:,jk) / (250. * wdust) )
170            zpdep   (:,:,jk) = zirondep(:,:,jk) * 0.38 / po4r
171         END DO
172
173         ! Solubilization of particles in the water column (Si, P, Fe)
174         tra(:,:,1,jpsil) = tra(:,:,1,jpsil) + zsidep  (:,:)
175         DO jk = 1, jpkm1
176            tra(:,:,jk,jppo4) = tra(:,:,jk,jppo4) + zpdep   (:,:,jk)
177            tra(:,:,jk,jpfer) = tra(:,:,jk,jpfer) + zirondep(:,:,jk) 
178         ENDDO
179         !
180         IF( lk_iomput ) THEN
181            IF( knt == nrdttrc ) THEN
182                IF( iom_use( "Irondep" ) )   &
183                &  CALL iom_put( "Irondep", zirondep(:,:,1) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron
184                IF( iom_use( "pdust" ) )   &
185                &  CALL iom_put( "pdust"  , dust(:,:) / ( wdust / rday )  * tmask(:,:,1) ) ! dust concentration at surface
186            ENDIF
187         ENDIF
188         DEALLOCATE( zsidep, zpdep, zirondep )
189         !                                             
190      ENDIF
191     
192      ! Add the external input of nutrients from river
193      ! ----------------------------------------------
194      IF( ln_river ) THEN
195         DO jj = 1, jpj
196            DO ji = 1, jpi
197               DO jk = 1, nk_rnf(ji,jj)
198                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) +  rivdip(ji,jj) * rfact2
199                  tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) +  rivdin(ji,jj) * rfact2
200                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) +  rivdic(ji,jj) * 5.e-5 * rfact2
201                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) +  rivdsi(ji,jj) * rfact2
202                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +  rivdic(ji,jj) * rfact2
203                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) +  ( rivalk(ji,jj) - rno3 * rivdin(ji,jj) ) * rfact2
204                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) +  rivdoc(ji,jj) * rfact2
205               ENDDO
206            ENDDO
207         ENDDO
208
209         ! When prognostic ligands are activated, ligands are supplied
210         ! to the ocean by rivers. We assume that the amount of ligands
211         ! is equal to that of iron (iron is completely complexed)
212         ! ------------------------------------------------------------
213         IF (ln_ligand) THEN
214            DO jj = 1, jpj
215               DO ji = 1, jpi
216                  DO jk = 1, nk_rnf(ji,jj)
217                     tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) +  rivdic(ji,jj) * 5.e-5 * rfact2
218                  ENDDO
219               ENDDO
220            ENDDO
221         ENDIF
222         ! PISCES-QUOTA part
223         IF( ln_p5z ) THEN
224            DO jj = 1, jpj
225               DO ji = 1, jpi
226                  DO jk = 1, nk_rnf(ji,jj)
227                     tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + rivdop(ji,jj) * rfact2
228                     tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + rivdon(ji,jj) * rfact2
229                  ENDDO
230               ENDDO
231            ENDDO
232         ENDIF
233      ENDIF
234     
235      ! Add the external input of nutrients from nitrogen deposition
236      ! ----------------------------------------------------------
237      IF( ln_ndepo ) THEN
238         tra(:,:,1,jpno3) = tra(:,:,1,jpno3) + nitdep(:,:) * rfact2
239         tra(:,:,1,jptal) = tra(:,:,1,jptal) - rno3 * nitdep(:,:) * rfact2
240      ENDIF
241
242      ! Add the external input of iron from hydrothermal vents
243      ! Please refer to Tagliabue et al. (2010) for more information
244      ! ------------------------------------------------------------
245      IF( ln_hydrofe ) THEN
246            tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2
247         IF( ln_ligand ) THEN
248            tra(:,:,:,jplgw) = tra(:,:,:,jplgw) + ( hydrofe(:,:,:) * lgw_rath ) * rfact2
249         ENDIF
250         !
251         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "HYDR" ) )   &
252            &   CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input
253      ENDIF
254
255      ! OA: Warning, the following part is necessary to avoid CFL problems
256      ! above the sediments. Vertical sinking speed is limited using the
257      ! typical CFL criterion
258      ! --------------------------------------------------------------------
259      DO jj = 1, jpj
260         DO ji = 1, jpi
261            ikt  = mbkt(ji,jj)
262            zdep = e3t_n(ji,jj,ikt) / xstep
263            zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) )
264            zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) )
265         END DO
266      END DO
267      !
268      ! No sediment module activated
269      IF( .NOT.lk_sed ) THEN
270!
271         ! Add the external input of iron from sediment mobilization
272         ! ------------------------------------------------------
273         IF( ln_ironsed ) THEN
274            tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2
275            !
276            IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironsed" ) )   &
277               &   CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! iron inputs from sediments
278         ENDIF
279
280         ! Computation of the sediment denitrification proportion: The metamodel
281         ! from Middleburg (2006) is used
282         ! Computation of the fraction of organic matter that is permanently
283         ! buried from Dunne's model (2007)
284         ! -------------------------------------------------------
285         DO jj = 1, jpj
286            DO ji = 1, jpi
287              IF( tmask(ji,jj,1) == 1 ) THEN
288                 ikt = mbkt(ji,jj)
289                 zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
290                   &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4
291                 zflx  = LOG10( MAX( 1E-3, zflx ) )
292                 zo2   = LOG10( MAX( 10. , trb(ji,jj,ikt,jpoxy) * 1E6 ) )
293                 zno3  = LOG10( MAX( 1.  , trb(ji,jj,ikt,jpno3) * 1E6 * rno3 ) )
294                 zdep  = LOG10( gdepw_n(ji,jj,ikt+1) )
295                 zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    &
296                   &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2
297                 zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) )
298                   !
299                 zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
300                   &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6
301                 zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2
302              ENDIF
303            END DO
304         END DO 
305         !
306      ENDIF
307
308      ! Fraction of dSi that is remineralized in the sediments. This is
309      ! set so that the burial in sediments equals the total input of Si
310      ! by rivers and dust (sedsilfrac)
311      ! ----------------------------------------------------------------
312      IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac
313
314      ! Loss of bSi and CaCO3 to the sediments
315      DO jj = 1, jpj
316         DO ji = 1, jpi
317            ikt  = mbkt(ji,jj)
318            zdep = xstep / e3t_n(ji,jj,ikt) 
319            zwsc = zwsbio4(ji,jj) * zdep
320            zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc
321            zcaloss = trb(ji,jj,ikt,jpcal) * zwsc
322            !
323            tra(ji,jj,ikt,jpgsi) = tra(ji,jj,ikt,jpgsi) - zsiloss
324            tra(ji,jj,ikt,jpcal) = tra(ji,jj,ikt,jpcal) - zcaloss
325         END DO
326      END DO
327      !
328      IF( .NOT.lk_sed ) THEN
329         ! Dissolution of CaCO3 and bSi in the sediments. This is
330         ! instantaneous since here sediments are not explicitly
331         ! modeled. The amount of CaCO3 that dissolves in the sediments
332         ! is computed using a metamodel constructed from Archer (1996)
333         ! ------------------------------------------------------------
334         DO jj = 1, jpj
335            DO ji = 1, jpi
336               ikt  = mbkt(ji,jj)
337               zdep = xstep / e3t_n(ji,jj,ikt) 
338               zwsc = zwsbio4(ji,jj) * zdep
339               zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc
340               zcaloss = trb(ji,jj,ikt,jpcal) * zwsc
341               tra(ji,jj,ikt,jpsil) = tra(ji,jj,ikt,jpsil) + zsiloss * zrivsil 
342               !
343               zfactcal = MIN( excess(ji,jj,ikt), 0.2 )
344               zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) )
345               zrivalk  = sedcalfrac * zfactcal
346               tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0
347               tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk
348               zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss * e3t_n(ji,jj,ikt) 
349               zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss * e3t_n(ji,jj,ikt) 
350            END DO
351         END DO
352      ENDIF
353      !
354      ! Loss of particulate organic carbon and Fe to the sediments
355      DO jj = 1, jpj
356         DO ji = 1, jpi
357            ikt  = mbkt(ji,jj)
358            zdep = xstep / e3t_n(ji,jj,ikt) 
359            zws4 = zwsbio4(ji,jj) * zdep
360            zws3 = zwsbio3(ji,jj) * zdep
361            tra(ji,jj,ikt,jpgoc) = tra(ji,jj,ikt,jpgoc) - trb(ji,jj,ikt,jpgoc) * zws4 
362            tra(ji,jj,ikt,jppoc) = tra(ji,jj,ikt,jppoc) - trb(ji,jj,ikt,jppoc) * zws3
363            tra(ji,jj,ikt,jpbfe) = tra(ji,jj,ikt,jpbfe) - trb(ji,jj,ikt,jpbfe) * zws4
364            tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3
365         END DO
366      END DO
367      !
368      ! Loss of particulate organic N and P to the sediments (p5z)
369      IF( ln_p5z ) THEN
370         DO jj = 1, jpj
371            DO ji = 1, jpi
372               ikt  = mbkt(ji,jj)
373               zdep = xstep / e3t_n(ji,jj,ikt) 
374               zws4 = zwsbio4(ji,jj) * zdep
375               zws3 = zwsbio3(ji,jj) * zdep
376               tra(ji,jj,ikt,jpgon) = tra(ji,jj,ikt,jpgon) - trb(ji,jj,ikt,jpgon) * zws4
377               tra(ji,jj,ikt,jppon) = tra(ji,jj,ikt,jppon) - trb(ji,jj,ikt,jppon) * zws3
378               tra(ji,jj,ikt,jpgop) = tra(ji,jj,ikt,jpgop) - trb(ji,jj,ikt,jpgop) * zws4
379               tra(ji,jj,ikt,jppop) = tra(ji,jj,ikt,jppop) - trb(ji,jj,ikt,jppop) * zws3
380            END DO
381         END DO
382      ENDIF
383
384      IF( .NOT.lk_sed ) THEN
385         ! Degradation of organic matter in the sediments. The metamodel of
386         ! Middleburg (2006) is used here to mimic the diagenetic reactions.
387         ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after
388         ! denitrification in the sediments. Not very clever, but simpliest option.
389         ! ------------------------------------------------------------------------
390         DO jj = 1, jpj
391            DO ji = 1, jpi
392               ikt  = mbkt(ji,jj)
393               zdep = xstep / e3t_n(ji,jj,ikt) 
394               zws4 = zwsbio4(ji,jj) * zdep
395               zws3 = zwsbio3(ji,jj) * zdep
396               zrivno3 = 1. - zbureff(ji,jj)
397               zwstpoc = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3
398               zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 )
399               z1pdenit = zwstpoc * zrivno3 - zpdenit
400               zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) )
401               tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit
402               tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit
403               tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit
404               tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) - rdenit * zpdenit
405               tra(ji,jj,ikt,jpoxy) = tra(ji,jj,ikt,jpoxy) - zolimit * o2ut
406               tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * zpdenit )
407               tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit 
408               sdenit(ji,jj) = rdenit * zpdenit * e3t_n(ji,jj,ikt)
409               zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t_n(ji,jj,ikt)
410               IF( ln_p5z ) THEN
411                  zwstpop              = trb(ji,jj,ikt,jpgop) * zws4 + trb(ji,jj,ikt,jppop) * zws3
412                  zwstpon              = trb(ji,jj,ikt,jpgon) * zws4 + trb(ji,jj,ikt,jppon) * zws3
413                  tra(ji,jj,ikt,jpdon) = tra(ji,jj,ikt,jpdon) + ( z1pdenit - zolimit ) * zwstpon / (zwstpoc + rtrn)
414                  tra(ji,jj,ikt,jpdop) = tra(ji,jj,ikt,jpdop) + ( z1pdenit - zolimit ) * zwstpop / (zwstpoc + rtrn)
415               ENDIF
416            END DO
417         END DO
418       ENDIF
419
420
421      ! Nitrogen fixation process : light limitation of diazotrophy
422      ! Small source of iron from particulate inorganic iron (photochemistry)
423      !----------------------------------------------------------------------
424      DO jk = 1, jpkm1
425         zlight (:,:,jk) =  ( 1.- EXP( -etot_ndcy(:,:,jk) / diazolight ) ) * ( 1. - fr_i(:,:) ) 
426         zsoufer(:,:,jk) = zlight(:,:,jk) * 2E-11 / ( 2E-11 + biron(:,:,jk) )
427      ENDDO
428
429      ! Diazotrophy (nitrogen fixation) is modeled according to an empirical
430      ! formulation. This is described in Aumont et al. (2015). Limitation
431      ! by P and Fe is computed. Inhibition by high N concentrations is imposed.
432      ! Diazotrophy sensitivity to temperature is parameterized as in
433      ! Ye et al. (2012) 
434      ! ------------------------------------------------------------------------
435      IF( ln_p4z ) THEN
436         ! PISCES part
437         DO jk = 1, jpkm1
438            DO jj = 1, jpj
439               DO ji = 1, jpi
440                  ! Potential nitrogen fixation dependant on temperature and iron
441                  ztemp = tsn(ji,jj,jk,jp_tem)
442                  zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
443                  xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) )
444                  xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
445                  zlim = ( 1.- xdiano3 - xdianh4 )
446                  ! Nitrogen fixation is almost fully halted when the N
447                  ! limitation term (xdiano3+xdianh4) is > 0.9
448                  IF( zlim <= 0.1 )   zlim = 0.01
449                  zfact = zlim * rfact2
450                  ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
451                  ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) )
452                  ztrdp = ztrpo4(ji,jj,jk)
453                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk)
454               END DO
455            END DO
456         END DO
457      ELSE       ! p5z
458         ! PISCES-QUOTA part
459         DO jk = 1, jpkm1
460            DO jj = 1, jpj
461               DO ji = 1, jpi
462                  ! Potential nitrogen fixation dependant on temperature and iron
463                  ztemp = tsn(ji,jj,jk,jp_tem)
464                  zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
465                  xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) )
466                  xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
467                  zlim = ( 1.- xdiano3 - xdianh4 )
468
469                  ! Nitrogen fixation is almost fully halted when the N
470                  ! limitation term (xdiano3+xdianh4) is > 0.9
471                  IF( zlim <= 0.1 )   zlim = 0.01
472                  zfact = zlim * rfact2
473                  ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
474                  ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) )
475                  ztrdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( 1E-6 + trb(ji,jj,jk,jpdop) ) * (1. - ztrpo4(ji,jj,jk))
476                  ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk)
477                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk)
478               END DO
479            END DO
480         END DO
481      ENDIF
482
483      ! Update of the TRA arrays due to nitrogen fixation
484      ! -------------------------------------------------
485      IF( ln_p4z ) THEN
486         ! PISCES part
487         DO jk = 1, jpkm1
488            DO jj = 1, jpj
489               DO ji = 1, jpi
490                  zfact = nitrpot(ji,jj,jk) * nitrfix
491                  ! 1/3 of the diazotrophs growth is supposed to be excreted
492                  ! as NH4. 1/3 as DOC and the rest is routed POC and GOC as
493                  ! a result of mortality by predation. Completely adhoc param
494                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0
495                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0
496                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zfact * 2.0 / 3.0
497                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zfact * 1.0 / 3.0
498                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zfact * 1.0 / 3.0 * 2.0 / 3.0
499                  tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0
500                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
501                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C
502                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0
503                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
504                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
505                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
506              END DO
507            END DO
508         END DO
509      ELSE    ! p5z
510         ! PISCES-QUOTA part
511         DO jk = 1, jpkm1
512            DO jj = 1, jpj
513               DO ji = 1, jpi
514                  zfact = nitrpot(ji,jj,jk) * nitrfix
515                  ! 1/3 of the diazotrophs growth is supposed to be excreted
516                  ! as NH4. 1/3 as DOC and the rest is routed POC and GOC as
517                  ! a result of mortality by predation. Completely adhoc param
518                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0
519                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0
520                  ! N/P ratio of diazotrophs is supposed to be 46
521                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) &
522                  &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
523                  tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zfact * 1.0 / 3.0
524                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zfact * 1.0 / 3.0
525                  tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + 16.0 / 46.0 * zfact / 3.0  &
526                  &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   &
527                  &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
528                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zfact * 1.0 / 3.0 * 2.0 / 3.0
529                  tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + zfact * 1.0 / 3.0 * 2.0 /3.0
530                  tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 2.0 /3.0
531                  tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0
532                  tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) + zfact * 1.0 / 3.0 * 1.0 /3.0
533                  tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0
534                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
535                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C
536                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0 
537                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
538                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
539                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
540              END DO
541            END DO
542         END DO
543         !
544      ENDIF
545
546      IF( lk_iomput ) THEN
547         IF( knt == nrdttrc ) THEN
548            zfact = 1.e+3 * rfact2r !  conversion from molC/l/kt  to molN/m3/s
549            IF( iom_use("Nfix"   ) ) CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation
550            IF( iom_use("INTNFIX") ) THEN   ! nitrogen fixation rate in ocean ( vertically integrated )
551               zwork(:,:) = 0.
552               DO jk = 1, jpkm1
553                 zwork(:,:) = zwork(:,:) + nitrpot(:,:,jk) * nitrfix * rno3 * zfact * e3t_n(:,:,jk) * tmask(:,:,jk)
554               ENDDO
555               CALL iom_put( "INTNFIX" , zwork ) 
556            ENDIF
557            IF( iom_use("SedCal" ) ) CALL iom_put( "SedCal", zsedcal(:,:) * zfact ) ! Permanent burial of CaCO3 in sediments
558            IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * zfact ) ! Permanent burial of bSi in sediments
559            IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * zfact ) ! Permanent burial of OC in sediments
560            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 ) ! Denitrification in the sediments
561         ENDIF
562      ENDIF
563      !
564      IF(ln_ctl) THEN  ! print mean trends (USEd for debugging)
565         WRITE(charout, fmt="('sed ')")
566         CALL prt_ctl_trc_info(charout)
567         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
568      ENDIF
569      !
570      IF( ln_p5z )    DEALLOCATE( ztrpo4, ztrdop )
571      !
572      IF( ln_timing )  CALL timing_stop('p4z_sed')
573      !
574   END SUBROUTINE p4z_sed
575
576
577   INTEGER FUNCTION p4z_sed_alloc()
578      !!----------------------------------------------------------------------
579      !!                     ***  ROUTINE p4z_sed_alloc  ***
580      !!----------------------------------------------------------------------
581      ALLOCATE( nitrpot(jpi,jpj,jpk), sdenit(jpi,jpj), STAT=p4z_sed_alloc )
582      !
583      IF( p4z_sed_alloc /= 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_alloc: failed to allocate arrays' )
584      !
585   END FUNCTION p4z_sed_alloc
586
587   !!======================================================================
588END MODULE p4zsed
Note: See TracBrowser for help on using the repository browser.