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/trunk/src/TOP/PISCES/P4Z – NEMO

source: NEMO/trunk/src/TOP/PISCES/P4Z/p4zsed.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 19.6 KB
Line 
1MODULE p4zsed
2   !!======================================================================
3   !!                         ***  MODULE p4sed  ***
4   !! TOP :   PISCES Compute loss of organic matter in the sediments
5   !!======================================================================
6   !! History :   1.0  !  2004-03 (O. Aumont) Original code
7   !!             2.0  !  2007-12 (C. Ethe, G. Madec)  F90
8   !!             3.4  !  2011-06 (C. Ethe) USE of fldread
9   !!             3.5  !  2012-07 (O. Aumont) improvment of river input of nutrients
10   !!----------------------------------------------------------------------
11   !!   p4z_sed        :  Compute loss of organic matter in the sediments
12   !!----------------------------------------------------------------------
13   USE oce_trc         !  shared variables between ocean and passive tracers
14   USE trc             !  passive tracers common variables
15   USE sms_pisces      !  PISCES Source Minus Sink variables
16   USE p4zlim          !  Co-limitations of differents nutrients
17   USE p4zint          !  interpolation and computation of various fields
18   USE sed             !  Sediment module
19   USE iom             !  I/O manager
20   USE prtctl_trc      !  print control for debugging
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   p4z_sed 
26   PUBLIC   p4z_sed_init
27   PUBLIC   p4z_sed_alloc
28 
29   REAL(wp), PUBLIC ::   nitrfix      !: Nitrogen fixation rate
30   REAL(wp), PUBLIC ::   diazolight   !: Nitrogen fixation sensitivty to light
31   REAL(wp), PUBLIC ::   concfediaz   !: Fe half-saturation Cste for diazotrophs
32
33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nitrpot    !: Nitrogen fixation
34   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: sdenit     !: Nitrate reduction in the sediments
35   !
36   REAL(wp), SAVE :: r1_rday         
37   REAL(wp), SAVE :: sedsilfrac, sedcalfrac
38
39   !! * Substitutions
40#  include "do_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
43   !! $Id$
44   !! Software governed by the CeCILL license (see ./LICENSE)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE p4z_sed( kt, knt, Kbb, Kmm, Krhs )
49      !!---------------------------------------------------------------------
50      !!                     ***  ROUTINE p4z_sed  ***
51      !!
52      !! ** Purpose :   Compute loss of organic matter in the sediments. This
53      !!              is by no way a sediment model. The loss is simply
54      !!              computed to balance the inout from rivers and dust
55      !!
56      !! ** Method  : - ???
57      !!---------------------------------------------------------------------
58      !
59      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
60      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices
61      INTEGER  ::  ji, jj, jk, ikt
62      REAL(wp) ::  zrivalk, zrivsil, zrivno3
63      REAL(wp) ::  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      !!---------------------------------------------------------------------
77      !
78      IF( ln_timing )  CALL timing_start('p4z_sed')
79      !
80
81      ! Allocate temporary workspace
82      ALLOCATE( ztrpo4(jpi,jpj,jpk) )
83      IF( ln_p5z )    ALLOCATE( ztrdop(jpi,jpj,jpk) )
84
85      zdenit2d(:,:) = 0.e0
86      zbureff (:,:) = 0.e0
87      zwork   (:,:) = 0.e0
88      zsedsi  (:,:) = 0.e0
89      zsedcal (:,:) = 0.e0
90      zsedc   (:,:) = 0.e0
91
92      IF( .NOT.lk_sed ) THEN
93         ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments
94         ! --------------------------------------------------------------------
95         DO_2D_11_11
96            ikt  = mbkt(ji,jj)
97            zdep = e3t(ji,jj,ikt,Kmm) / xstep
98            zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) )
99            zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) )
100         END_2D
101
102         ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used
103         ! Computation of the fraction of organic matter that is permanently buried from Dunne's model
104         ! -------------------------------------------------------
105         DO_2D_11_11
106           IF( tmask(ji,jj,1) == 1 ) THEN
107              ikt = mbkt(ji,jj)
108              zflx = (  tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj)   &
109                &     + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4
110              zflx  = LOG10( MAX( 1E-3, zflx ) )
111              zo2   = LOG10( MAX( 10. , tr(ji,jj,ikt,jpoxy,Kbb) * 1E6 ) )
112              zno3  = LOG10( MAX( 1.  , tr(ji,jj,ikt,jpno3,Kbb) * 1E6 * rno3 ) )
113              zdep  = LOG10( gdepw(ji,jj,ikt+1,Kmm) )
114              zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    &
115                &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2
116              zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) )
117                !
118              zflx = (  tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj)   &
119                &     + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) ) * 1E6
120              zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2
121           ENDIF
122         END_2D
123         !
124      ENDIF
125
126      ! This loss is scaled at each bottom grid cell for equilibrating the total budget of silica in the ocean.
127      ! Thus, the amount of silica lost in the sediments equal the supply at the surface (dust+rivers)
128      ! ------------------------------------------------------
129      IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac
130
131      DO_2D_11_11
132         ikt  = mbkt(ji,jj)
133         zdep = xstep / e3t(ji,jj,ikt,Kmm) 
134         zwsc = zwsbio4(ji,jj) * zdep
135         zsiloss = tr(ji,jj,ikt,jpgsi,Kbb) * zwsc
136         zcaloss = tr(ji,jj,ikt,jpcal,Kbb) * zwsc
137         !
138         tr(ji,jj,ikt,jpgsi,Krhs) = tr(ji,jj,ikt,jpgsi,Krhs) - zsiloss
139         tr(ji,jj,ikt,jpcal,Krhs) = tr(ji,jj,ikt,jpcal,Krhs) - zcaloss
140      END_2D
141      !
142      IF( .NOT.lk_sed ) THEN
143         DO_2D_11_11
144            ikt  = mbkt(ji,jj)
145            zdep = xstep / e3t(ji,jj,ikt,Kmm) 
146            zwsc = zwsbio4(ji,jj) * zdep
147            zsiloss = tr(ji,jj,ikt,jpgsi,Kbb) * zwsc
148            zcaloss = tr(ji,jj,ikt,jpcal,Kbb) * zwsc
149            tr(ji,jj,ikt,jpsil,Krhs) = tr(ji,jj,ikt,jpsil,Krhs) + zsiloss * zrivsil 
150            !
151            zfactcal = MIN( excess(ji,jj,ikt), 0.2 )
152            zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) )
153            zrivalk  = sedcalfrac * zfactcal
154            tr(ji,jj,ikt,jptal,Krhs) =  tr(ji,jj,ikt,jptal,Krhs) + zcaloss * zrivalk * 2.0
155            tr(ji,jj,ikt,jpdic,Krhs) =  tr(ji,jj,ikt,jpdic,Krhs) + zcaloss * zrivalk
156            zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss * e3t(ji,jj,ikt,Kmm) 
157            zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss * e3t(ji,jj,ikt,Kmm) 
158         END_2D
159      ENDIF
160      !
161      DO_2D_11_11
162         ikt  = mbkt(ji,jj)
163         zdep = xstep / e3t(ji,jj,ikt,Kmm) 
164         zws4 = zwsbio4(ji,jj) * zdep
165         zws3 = zwsbio3(ji,jj) * zdep
166         tr(ji,jj,ikt,jpgoc,Krhs) = tr(ji,jj,ikt,jpgoc,Krhs) - tr(ji,jj,ikt,jpgoc,Kbb) * zws4 
167         tr(ji,jj,ikt,jppoc,Krhs) = tr(ji,jj,ikt,jppoc,Krhs) - tr(ji,jj,ikt,jppoc,Kbb) * zws3
168         tr(ji,jj,ikt,jpbfe,Krhs) = tr(ji,jj,ikt,jpbfe,Krhs) - tr(ji,jj,ikt,jpbfe,Kbb) * zws4
169         tr(ji,jj,ikt,jpsfe,Krhs) = tr(ji,jj,ikt,jpsfe,Krhs) - tr(ji,jj,ikt,jpsfe,Kbb) * zws3
170      END_2D
171      !
172      IF( ln_p5z ) THEN
173         DO_2D_11_11
174            ikt  = mbkt(ji,jj)
175            zdep = xstep / e3t(ji,jj,ikt,Kmm) 
176            zws4 = zwsbio4(ji,jj) * zdep
177            zws3 = zwsbio3(ji,jj) * zdep
178            tr(ji,jj,ikt,jpgon,Krhs) = tr(ji,jj,ikt,jpgon,Krhs) - tr(ji,jj,ikt,jpgon,Kbb) * zws4
179            tr(ji,jj,ikt,jppon,Krhs) = tr(ji,jj,ikt,jppon,Krhs) - tr(ji,jj,ikt,jppon,Kbb) * zws3
180            tr(ji,jj,ikt,jpgop,Krhs) = tr(ji,jj,ikt,jpgop,Krhs) - tr(ji,jj,ikt,jpgop,Kbb) * zws4
181            tr(ji,jj,ikt,jppop,Krhs) = tr(ji,jj,ikt,jppop,Krhs) - tr(ji,jj,ikt,jppop,Kbb) * zws3
182         END_2D
183      ENDIF
184
185      IF( .NOT.lk_sed ) THEN
186         ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after
187         ! denitrification in the sediments. Not very clever, but simpliest option.
188         DO_2D_11_11
189            ikt  = mbkt(ji,jj)
190            zdep = xstep / e3t(ji,jj,ikt,Kmm) 
191            zws4 = zwsbio4(ji,jj) * zdep
192            zws3 = zwsbio3(ji,jj) * zdep
193            zrivno3 = 1. - zbureff(ji,jj)
194            zwstpoc = tr(ji,jj,ikt,jpgoc,Kbb) * zws4 + tr(ji,jj,ikt,jppoc,Kbb) * zws3
195            zpdenit  = MIN( 0.5 * ( tr(ji,jj,ikt,jpno3,Kbb) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 )
196            z1pdenit = zwstpoc * zrivno3 - zpdenit
197            zolimit = MIN( ( tr(ji,jj,ikt,jpoxy,Kbb) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) )
198            tr(ji,jj,ikt,jpdoc,Krhs) = tr(ji,jj,ikt,jpdoc,Krhs) + z1pdenit - zolimit
199            tr(ji,jj,ikt,jppo4,Krhs) = tr(ji,jj,ikt,jppo4,Krhs) + zpdenit + zolimit
200            tr(ji,jj,ikt,jpnh4,Krhs) = tr(ji,jj,ikt,jpnh4,Krhs) + zpdenit + zolimit
201            tr(ji,jj,ikt,jpno3,Krhs) = tr(ji,jj,ikt,jpno3,Krhs) - rdenit * zpdenit
202            tr(ji,jj,ikt,jpoxy,Krhs) = tr(ji,jj,ikt,jpoxy,Krhs) - zolimit * o2ut
203            tr(ji,jj,ikt,jptal,Krhs) = tr(ji,jj,ikt,jptal,Krhs) + rno3 * (zolimit + (1.+rdenit) * zpdenit )
204            tr(ji,jj,ikt,jpdic,Krhs) = tr(ji,jj,ikt,jpdic,Krhs) + zpdenit + zolimit 
205            sdenit(ji,jj) = rdenit * zpdenit * e3t(ji,jj,ikt,Kmm)
206            zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t(ji,jj,ikt,Kmm)
207            IF( ln_p5z ) THEN
208               zwstpop              = tr(ji,jj,ikt,jpgop,Kbb) * zws4 + tr(ji,jj,ikt,jppop,Kbb) * zws3
209               zwstpon              = tr(ji,jj,ikt,jpgon,Kbb) * zws4 + tr(ji,jj,ikt,jppon,Kbb) * zws3
210               tr(ji,jj,ikt,jpdon,Krhs) = tr(ji,jj,ikt,jpdon,Krhs) + ( z1pdenit - zolimit ) * zwstpon / (zwstpoc + rtrn)
211               tr(ji,jj,ikt,jpdop,Krhs) = tr(ji,jj,ikt,jpdop,Krhs) + ( z1pdenit - zolimit ) * zwstpop / (zwstpoc + rtrn)
212            ENDIF
213         END_2D
214       ENDIF
215
216
217      ! Nitrogen fixation process
218      ! Small source iron from particulate inorganic iron
219      !-----------------------------------
220      DO jk = 1, jpkm1
221         zlight (:,:,jk) =  ( 1.- EXP( -etot_ndcy(:,:,jk) / diazolight ) ) * ( 1. - fr_i(:,:) ) 
222         zsoufer(:,:,jk) = zlight(:,:,jk) * 2E-11 / ( 2E-11 + biron(:,:,jk) )
223      ENDDO
224      IF( ln_p4z ) THEN
225         DO_3D_11_11( 1, jpkm1 )
226            !                      ! Potential nitrogen fixation dependant on temperature and iron
227            ztemp = ts(ji,jj,jk,jp_tem,Kmm)
228            zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
229            !       Potential nitrogen fixation dependant on temperature and iron
230            xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) )
231            xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4)
232            zlim = ( 1.- xdiano3 - xdianh4 )
233            IF( zlim <= 0.1 )   zlim = 0.01
234            zfact = zlim * rfact2
235            ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
236            ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) )
237            ztrdp = ztrpo4(ji,jj,jk)
238            nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk)
239         END_3D
240      ELSE       ! p5z
241         DO_3D_11_11( 1, jpkm1 )
242            !                      ! Potential nitrogen fixation dependant on temperature and iron
243            ztemp = ts(ji,jj,jk,jp_tem,Kmm)
244            zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
245            !       Potential nitrogen fixation dependant on temperature and iron
246            xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) )
247            xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4)
248            zlim = ( 1.- xdiano3 - xdianh4 )
249            IF( zlim <= 0.1 )   zlim = 0.01
250            zfact = zlim * rfact2
251            ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
252            ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) )
253            ztrdop(ji,jj,jk) = tr(ji,jj,jk,jpdop,Kbb) / ( 1E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (1. - ztrpo4(ji,jj,jk))
254            ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk)
255            nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk)
256         END_3D
257      ENDIF
258
259      ! Nitrogen change due to nitrogen fixation
260      ! ----------------------------------------
261      IF( ln_p4z ) THEN
262         DO_3D_11_11( 1, jpkm1 )
263            zfact = nitrpot(ji,jj,jk) * nitrfix
264            tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0
265            tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0
266            tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - zfact * 2.0 / 3.0
267            tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zfact * 1.0 / 3.0
268            tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zfact * 1.0 / 3.0 * 2.0 / 3.0
269            tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfact * 1.0 / 3.0 * 1.0 / 3.0
270            tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
271            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0
272            tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
273            tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
274            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
275            tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + concdnh4 / ( concdnh4 + tr(ji,jj,jk,jppo4,Kbb) ) &
276            &                     * 0.001 * tr(ji,jj,jk,jpdoc,Kbb) * xstep
277         END_3D
278      ELSE    ! p5z
279         DO_3D_11_11( 1, jpkm1 )
280            zfact = nitrpot(ji,jj,jk) * nitrfix
281            tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0
282            tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0
283            tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) &
284            &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
285            tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zfact * 1.0 / 3.0
286            tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zfact * 1.0 / 3.0
287            tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + 16.0 / 46.0 * zfact / 3.0  &
288            &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   &
289            &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
290            tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zfact * 1.0 / 3.0 * 2.0 / 3.0
291            tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) + zfact * 1.0 / 3.0 * 2.0 /3.0
292            tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 2.0 /3.0
293            tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfact * 1.0 / 3.0 * 1.0 / 3.0
294            tr(ji,jj,jk,jpgon,Krhs) = tr(ji,jj,jk,jpgon,Krhs) + zfact * 1.0 / 3.0 * 1.0 /3.0
295            tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0
296            tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
297            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0 
298            tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
299            tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
300            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
301         END_3D
302         !
303      ENDIF
304
305      IF( lk_iomput .AND. knt == nrdttrc ) THEN
306         zfact = 1.e+3 * rfact2r !  conversion from molC/l/kt  to molN/m3/s
307         CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation
308         CALL iom_put( "SedCal", zsedcal(:,:) * zfact )
309         CALL iom_put( "SedSi" , zsedsi (:,:) * zfact )
310         CALL iom_put( "SedC"  , zsedc  (:,:) * zfact )
311         CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 )
312      ENDIF
313      !
314      IF(sn_cfctl%l_prttrc) THEN  ! print mean trends (USEd for debugging)
315         WRITE(charout, fmt="('sed ')")
316         CALL prt_ctl_trc_info(charout)
317         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
318      ENDIF
319      !
320      IF( ln_p5z )    DEALLOCATE( ztrpo4, ztrdop )
321      !
322      IF( ln_timing )  CALL timing_stop('p4z_sed')
323      !
324   END SUBROUTINE p4z_sed
325
326   SUBROUTINE p4z_sed_init
327      !!----------------------------------------------------------------------
328      !!                  ***  routine p4z_sed_init  ***
329      !!
330      !! ** purpose :   initialization of some parameters
331      !!
332      !!----------------------------------------------------------------------
333      !!----------------------------------------------------------------------
334      INTEGER  :: ji, jj, jk, jm
335      INTEGER  :: ios                 ! Local integer output status for namelist read
336      !
337      !!
338      NAMELIST/nampissed/ nitrfix, diazolight, concfediaz
339      !!----------------------------------------------------------------------
340      !
341      IF(lwp) THEN
342         WRITE(numout,*)
343         WRITE(numout,*) 'p4z_sed_init : initialization of sediment mobilisation '
344         WRITE(numout,*) '~~~~~~~~~~~~ '
345      ENDIF
346      !                            !* set file information
347      READ  ( numnatp_ref, nampissed, IOSTAT = ios, ERR = 901)
348901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampissed in reference namelist' )
349      READ  ( numnatp_cfg, nampissed, IOSTAT = ios, ERR = 902 )
350902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampissed in configuration namelist' )
351      IF(lwm) WRITE ( numonp, nampissed )
352
353      IF(lwp) THEN
354         WRITE(numout,*) '   Namelist : nampissed '
355         WRITE(numout,*) '      nitrogen fixation rate                       nitrfix = ', nitrfix
356         WRITE(numout,*) '      nitrogen fixation sensitivty to light    diazolight  = ', diazolight
357         WRITE(numout,*) '      Fe half-saturation cste for diazotrophs  concfediaz  = ', concfediaz
358      ENDIF
359      !
360      r1_rday  = 1. / rday
361      !
362      sedsilfrac = 0.03     ! percentage of silica loss in the sediments
363      sedcalfrac = 0.6      ! percentage of calcite loss in the sediments
364      !
365      lk_sed = ln_sediment .AND. ln_sed_2way 
366      !
367   END SUBROUTINE p4z_sed_init
368
369   INTEGER FUNCTION p4z_sed_alloc()
370      !!----------------------------------------------------------------------
371      !!                     ***  ROUTINE p4z_sed_alloc  ***
372      !!----------------------------------------------------------------------
373      ALLOCATE( nitrpot(jpi,jpj,jpk), sdenit(jpi,jpj), STAT=p4z_sed_alloc )
374      !
375      IF( p4z_sed_alloc /= 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_alloc: failed to allocate arrays' )
376      !
377   END FUNCTION p4z_sed_alloc
378
379   !!======================================================================
380END MODULE p4zsed
Note: See TracBrowser for help on using the repository browser.