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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zsed.F90 @ 12340

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

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • 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.