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

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

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