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

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

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