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.
p5zprod.F90 in branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z – NEMO

source: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zprod.F90 @ 7190

Last change on this file since 7190 was 7190, checked in by aumont, 7 years ago

various important bug fixes
  • Property svn:executable set to *
File size: 35.7 KB
Line 
1MODULE p5zprod
2   !!======================================================================
3   !!                         ***  MODULE p5zprod  ***
4   !! TOP :  Growth Rate of the two phytoplanktons groups
5   !!======================================================================
6   !! History :   1.0  !  2004     (O. Aumont) Original code
7   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
8   !!             3.4  !  2011-05  (O. Aumont, C. Ethe) New parameterization of light limitation
9   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
10   !!----------------------------------------------------------------------
11#if defined key_pisces_quota
12   !!----------------------------------------------------------------------
13   !!   'key_pisces_quota'     PISCES bio-model with variable stoichiometry
14   !!----------------------------------------------------------------------
15   !!   p5z_prod       :   Compute the growth Rate of the two phytoplanktons groups
16   !!   p5z_prod_init  :   Initialization of the parameters for growth
17   !!   p5z_prod_alloc :   Allocate variables for growth
18   !!----------------------------------------------------------------------
19   USE oce_trc         !  shared variables between ocean and passive tracers
20   USE trc             !  passive tracers common variables
21   USE sms_pisces      !  PISCES Source Minus Sink variables
22   USE p4zopt          !  optical model
23   USE p5zlim          !  Co-limitations of differents nutrients
24   USE prtctl_trc      !  print control for debugging
25   USE iom             !  I/O manager
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   p5z_prod         ! called in p5zbio.F90
31   PUBLIC   p5z_prod_init    ! called in trcsms_pisces.F90
32   PUBLIC   p5z_prod_alloc
33
34   !! * Shared module variables
35   REAL(wp), PUBLIC ::  pislopen        !:
36   REAL(wp), PUBLIC ::  pislopep        !:
37   REAL(wp), PUBLIC ::  pisloped        !:
38   REAL(wp), PUBLIC ::  xadap           !:
39   REAL(wp), PUBLIC ::  excretn         !:
40   REAL(wp), PUBLIC ::  excretp         !:
41   REAL(wp), PUBLIC ::  excretd         !:
42   REAL(wp), PUBLIC ::  bresp           !:
43   REAL(wp), PUBLIC ::  thetanpm        !:
44   REAL(wp), PUBLIC ::  thetannm        !:
45   REAL(wp), PUBLIC ::  thetandm        !:
46   REAL(wp), PUBLIC ::  chlcmin         !:
47   REAL(wp), PUBLIC ::  grosip          !:
48
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prmaxn   !: optimal production = f(temperature)
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prmaxp   !: optimal production = f(temperature)
51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prmaxd   !: optimal production = f(temperature)
52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   zdaylen
53   
54   REAL(wp) :: r1_rday                !: 1 / rday
55   REAL(wp) :: texcretn               !: 1 - excret
56   REAL(wp) :: texcretp               !: 1 - excretp
57   REAL(wp) :: texcretd               !: 1 - excret2       
58
59
60   !!* Substitution
61#  include "top_substitute.h90"
62   !!----------------------------------------------------------------------
63   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
64   !! $Id: p4zprod.F90 3160 2011-11-20 14:27:18Z cetlod $
65   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE p5z_prod( kt , knt )
70      !!---------------------------------------------------------------------
71      !!                     ***  ROUTINE p5z_prod  ***
72      !!
73      !! ** Purpose :   Compute the phytoplankton production depending on
74      !!              light, temperature and nutrient availability
75      !!
76      !! ** Method  : - ???
77      !!---------------------------------------------------------------------
78      !
79      INTEGER, INTENT(in) :: kt, knt
80      !
81      INTEGER  ::   ji, jj, jk
82      REAL(wp) ::   zsilfac, znanotot, zpicotot, zdiattot, zconctemp, zconctemp2
83      REAL(wp) ::   zration, zratiop, zratiof, zmax, zmax2, zsilim, ztn, zadap
84      REAL(wp) ::   zpronmax, zpropmax, zprofmax, zrat
85      REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zprontot, zproptot, zprodtot
86      REAL(wp) ::   zprnutmax, zdocprod, zprochln, zprochld, zprochlp
87      REAL(wp) ::   zpislopen, zpislopep, zpisloped, thetannm_n, thetandm_n, thetanpm_n
88      REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
89      REAL(wp) ::   zrfact2
90      CHARACTER (len=25) :: charout
91      REAL(wp), POINTER, DIMENSION(:,:  ) :: zmixnano, zmixpico, zmixdiat, zstrn
92      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpislopeadn, zpislopeadp, zpislopeadd
93      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprbio, zprpic, zprdia, zysopt
94      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprchln, zprchlp, zprchld
95      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprorcan, zprorcap, zprorcad 
96      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprofed, zprofep, zprofen
97      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpronewn, zpronewp, zpronewd
98      REAL(wp), POINTER, DIMENSION(:,:,:) :: zproregn, zproregp, zproregd
99      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpropo4n, zpropo4p, zpropo4d
100      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprodopn, zprodopp, zprodopd
101      REAL(wp), POINTER, DIMENSION(:,:,:) :: zrespn, zrespp, zrespd, zprnut
102      REAL(wp), POINTER, DIMENSION(:,:,:) :: zcroissn, zcroissp, zcroissd
103      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmxl_fac, zmxl_chl
104      !!---------------------------------------------------------------------
105      !
106      IF( nn_timing == 1 )  CALL timing_start('p5z_prod')
107      !
108      !  Allocate temporary workspace
109      CALL wrk_alloc( jpi, jpj,      zmixnano, zmixpico, zmixdiat, zstrn )
110      CALL wrk_alloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
111      CALL wrk_alloc( jpi, jpj, jpk, zpislopeadn, zpislopeadp, zpislopeadd, zysopt ) 
112      CALL wrk_alloc( jpi, jpj, jpk, zprdia, zprpic, zprbio, zprorcan, zprorcap, zprorcad )
113      CALL wrk_alloc( jpi, jpj, jpk, zprofed, zprofep, zprofen )
114      CALL wrk_alloc( jpi, jpj, jpk, zpronewn, zpronewp, zpronewd, zproregn, zproregp, zproregd )
115      CALL wrk_alloc( jpi, jpj, jpk, zpropo4n, zpropo4p, zpropo4d, zrespn, zrespp, zrespd, zprnut )
116      CALL wrk_alloc( jpi, jpj, jpk, zprchln, zprchlp, zprchld, zprodopn, zprodopp, zprodopd )
117      CALL wrk_alloc( jpi, jpj, jpk, zcroissp, zcroissn, zcroissd )
118      !
119      zprorcan(:,:,:) = 0._wp ; zprorcap(:,:,:) = 0._wp ; zprorcad(:,:,:) = 0._wp
120      zprofed (:,:,:) = 0._wp ; zprofep (:,:,:) = 0._wp ; zprofen (:,:,:) = 0._wp
121      zpronewn(:,:,:) = 0._wp ; zpronewp(:,:,:) = 0._wp ; zpronewd(:,:,:) = 0._wp
122      zproregn(:,:,:) = 0._wp ; zproregp(:,:,:) = 0._wp ; zproregd(:,:,:) = 0._wp 
123      zpropo4n(:,:,:) = 0._wp ; zpropo4p(:,:,:) = 0._wp ; zpropo4d(:,:,:) = 0._wp
124      zprdia  (:,:,:) = 0._wp ; zprpic  (:,:,:) = 0._wp ; zprbio  (:,:,:) = 0._wp
125      zysopt  (:,:,:) = 0._wp
126      zrespn  (:,:,:) = 0._wp ; zrespp  (:,:,:) = 0._wp ; zrespd  (:,:,:) = 0._wp 
127
128      ! Computation of the optimal production
129      prmaxn(:,:,:) = ( 0.65_wp * (1. + zpsino3 * qnpmax ) ) * r1_rday * tgfunc(:,:,:)
130      prmaxp(:,:,:) = 0.5 / 0.65 * prmaxn(:,:,:) 
131      prmaxd(:,:,:) = prmaxn(:,:,:) 
132      zprnut(:,:,:) = 0.65_wp * r1_rday * tgfunc(:,:,:)
133
134      IF( lk_degrad )  THEN
135         prmaxn(:,:,:) = prmaxn(:,:,:) * facvol(:,:,:) 
136         prmaxp(:,:,:) = prmaxp(:,:,:) * facvol(:,:,:)
137         prmaxd(:,:,:) = prmaxd(:,:,:) * facvol(:,:,:)
138      ENDIF
139
140      ! compute the day length depending on latitude and the day
141      zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
142      zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
143
144      ! day length in hours
145      zstrn(:,:) = 0.
146      DO jj = 1, jpj
147         DO ji = 1, jpi
148            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
149            zargu = MAX( -1., MIN(  1., zargu ) )
150            zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
151         END DO
152      END DO
153
154         ! Impact of the day duration on phytoplankton growth
155      DO jk = 1, jpkm1
156         DO jj = 1 ,jpj
157            DO ji = 1, jpi
158               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
159                  zval = MAX( 1., zstrn(ji,jj) )
160                  IF( fsdept(ji,jj,jk) <= hmld(ji,jj) ) THEN
161                     zval = MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
162                  ENDIF
163                  zmxl_chl(ji,jj,jk) = zval / 24.
164                  zmxl_fac(ji,jj,jk) = 1.5 * zval / ( 12. + zval )
165               ENDIF
166            END DO
167         END DO
168      END DO
169
170      zprbio(:,:,:) = prmaxn(:,:,:) * zmxl_fac(:,:,:)
171      zprdia(:,:,:) = prmaxd(:,:,:) * zmxl_fac(:,:,:)
172      zprpic(:,:,:) = prmaxp(:,:,:) * zmxl_fac(:,:,:)
173
174
175      ! Maximum light intensity
176      zdaylen(:,:) = MAX(1., zstrn(:,:)) / 24.
177      WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
178
179      DO jk = 1, jpkm1
180!CDIR NOVERRCHK
181         DO jj = 1, jpj
182!CDIR NOVERRCHK
183            DO ji = 1, jpi
184               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
185                  ! Computation of the P-I slope for nanos and diatoms
186                  ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
187                  zadap       = xadap * ztn / ( 2.+ ztn )
188                  !
189                  zpislopeadn(ji,jj,jk) = pislopen * trb(ji,jj,jk,jpnch)    &
190                  &                       /( trb(ji,jj,jk,jpphy) * 12. + rtrn)
191                  zpislopeadp(ji,jj,jk) = pislopep * ( 1. + zadap * EXP( -0.25 * epico(ji,jj,jk) ) )   &
192                  &                       * trb(ji,jj,jk,jppch) /( trb(ji,jj,jk,jppic) * 12. + rtrn)
193                  zpislopeadd(ji,jj,jk) = pisloped * trb(ji,jj,jk,jpdch)    &
194                     &                    /( trb(ji,jj,jk,jpdia) * 12. + rtrn)
195                  !
196                  zpislopen = zpislopeadn(ji,jj,jk) / ( zprbio(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
197                  zpislopep = zpislopeadp(ji,jj,jk) / ( zprpic(ji,jj,jk) * rday * xlimpic(ji,jj,jk) + rtrn )
198                  zpisloped = zpislopeadd(ji,jj,jk) / ( zprdia(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
199
200                  ! Computation of production function for Carbon
201                  !  ---------------------------------------------
202                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
203                  zprpic(ji,jj,jk) = zprpic(ji,jj,jk) * ( 1.- EXP( -zpislopep * epico(ji,jj,jk) )  )
204                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
205
206                  ! Computation of production function for Chlorophyll
207                  !  -------------------------------------------------
208                  zpislopen = zpislopen * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
209                  zpisloped = zpisloped * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
210                  zpislopep = zpislopep * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
211                  zprchln(ji,jj,jk) = prmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
212                  zprchlp(ji,jj,jk) = prmaxp(ji,jj,jk) * ( 1.- EXP( -zpislopep * epico(ji,jj,jk) )  )
213                  zprchld(ji,jj,jk) = prmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
214               ENDIF
215            END DO
216         END DO
217      END DO
218
219      DO jk = 1, jpkm1
220         DO jj = 1, jpj
221            DO ji = 1, jpi
222
223                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
224                  !    Si/C of diatoms
225                  !    ------------------------
226                  !    Si/C increases with iron stress and silicate availability
227                  !    Si/C is arbitrariliy increased for very high Si concentrations
228                  !    to mimic the very high ratios observed in the Southern Ocean (silpot2)
229                  zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
230                  zsilim = MIN( zprdia(ji,jj,jk) / ( prmaxd(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) )
231                  zsilfac = 3.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim - 0.5 ) )  ) + 1.e0
232                  zsiborn = trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil)
233                  IF (gphit(ji,jj) < -30 ) THEN
234                    zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
235                  ELSE
236                    zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 )
237                  ENDIF
238                  zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2
239              ENDIF
240            END DO
241         END DO
242      END DO
243
244      !  Sea-ice effect on production                                                                               
245      DO jk = 1, jpkm1
246         DO jj = 1, jpj
247            DO ji = 1, jpi
248               zprbio(ji,jj,jk)  = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
249               zprpic(ji,jj,jk)  = zprpic(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
250               zprdia(ji,jj,jk)  = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
251               zprnut(ji,jj,jk)  = zprnut(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
252            END DO
253         END DO
254      END DO
255
256      ! Computation of the various production terms of nanophytoplankton
257      DO jk = 1, jpkm1
258         DO jj = 1, jpj
259            DO ji = 1, jpi
260               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
261                  !  production terms for nanophyto.
262                  zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trb(ji,jj,jk,jpphy) * rfact2
263                  !
264                  zration = trb(ji,jj,jk,jpnph) / ( trb(ji,jj,jk,jpphy) + rtrn )
265                  zratiop = trb(ji,jj,jk,jppph) / ( trb(ji,jj,jk,jpphy) + rtrn )
266                  zratiof = trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) + rtrn )
267                  zprnutmax = zprnut(ji,jj,jk) * fvnuptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jpphy) * rfact2
268                  ! Uptake of nitrogen
269                  zrat = MIN( 1., zration / (xqnnmax(ji,jj,jk) + rtrn) ) 
270                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
271                  zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqpnmin(ji,jj,jk) )   &
272                  &          / ( xqpnmax(ji,jj,jk) - xqpnmin(ji,jj,jk) ), xlimnfe(ji,jj,jk) ) )
273                  zpronewn(ji,jj,jk) = zpronmax * zdaylen(ji,jj) * xnanono3(ji,jj,jk)
274                  zproregn(ji,jj,jk) = zpronmax * xnanonh4(ji,jj,jk)
275                  ! Uptake of phosphorus
276                  zrat = MIN( 1., zratiop / (xqpnmax(ji,jj,jk) + rtrn) )
277                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
278                  zpropmax = zprnutmax * zmax * xlimnfe(ji,jj,jk)
279                  zpropo4n(ji,jj,jk) = zpropmax * xnanopo4(ji,jj,jk)
280                  zprodopn(ji,jj,jk) = zpropmax * xnanodop(ji,jj,jk)
281                  ! Uptake of iron
282                  zrat = MIN( 1., zratiof / qfnmax )
283                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
284                  zprofmax = zprnutmax * qfnmax * zmax
285                  zprofen(ji,jj,jk) = zprofmax * xnanofer(ji,jj,jk) * ( 3. - 2.4 * xlimnfe(ji,jj,jk)    &
286                  &          / ( xlimnfe(ji,jj,jk) + 0.2 ) ) * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn  &
287                  &          + xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) ) * (1. - xnanofer(ji,jj,jk) ) )
288               ENDIF
289            END DO
290         END DO
291      END DO
292
293      ! Computation of the various production terms of picophytoplankton
294      DO jk = 1, jpkm1
295         DO jj = 1, jpj
296            DO ji = 1, jpi
297               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
298                  !  production terms for picophyto.
299                  zprorcap(ji,jj,jk) = zprpic(ji,jj,jk)  * xlimpic(ji,jj,jk) * trb(ji,jj,jk,jppic) * rfact2
300                  !
301                  zration = trb(ji,jj,jk,jpnpi) / ( trb(ji,jj,jk,jppic) + rtrn )
302                  zratiop = trb(ji,jj,jk,jpppi) / ( trb(ji,jj,jk,jppic) + rtrn )
303                  zratiof = trb(ji,jj,jk,jppfe) / ( trb(ji,jj,jk,jppic) + rtrn )
304                  zprnutmax = zprnut(ji,jj,jk) * fvpuptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jppic) * rfact2
305                  ! Uptake of nitrogen
306                  zrat = MIN( 1., zration / (xqnpmax(ji,jj,jk) + rtrn) )
307                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
308                  zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqppmin(ji,jj,jk) )   &
309                  &          / ( xqppmax(ji,jj,jk) - xqppmin(ji,jj,jk) ), xlimpfe(ji,jj,jk) ) )
310                  zpronewp(ji,jj,jk) = zpronmax * zdaylen(ji,jj) * xpicono3(ji,jj,jk) 
311                  zproregp(ji,jj,jk) = zpronmax * xpiconh4(ji,jj,jk)
312                  ! Uptake of phosphorus
313                  zrat = MIN( 1., zratiop / (xqppmax(ji,jj,jk) + rtrn) )
314                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
315                  zpropmax = zprnutmax * zmax * xlimpfe(ji,jj,jk)
316                  zpropo4p(ji,jj,jk) = zpropmax * xpicopo4(ji,jj,jk)
317                  zprodopp(ji,jj,jk) = zpropmax * xpicodop(ji,jj,jk)
318                  ! Uptake of iron
319                  zrat = MIN( 1., zratiof / qfpmax )
320                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
321                  zprofmax = zprnutmax * qfpmax * zmax
322                  zprofep(ji,jj,jk) = zprofmax * xpicofer(ji,jj,jk) * ( 3. - 2.4 * xlimpfe(ji,jj,jk)   &
323                  &          / ( xlimpfe(ji,jj,jk) + 0.2 ) ) * (1. + 0.8 * xpicono3(ji,jj,jk) / ( rtrn   &
324                  &          + xpicono3(ji,jj,jk) + xpiconh4(ji,jj,jk) ) * (1. - xpicofer(ji,jj,jk) ) )
325               ENDIF
326            END DO
327         END DO
328      END DO
329
330      ! Computation of the various production terms of diatoms
331      DO jk = 1, jpkm1
332         DO jj = 1, jpj
333            DO ji = 1, jpi
334               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
335                  !  production terms for diatomees
336                  zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trb(ji,jj,jk,jpdia) * rfact2
337                  ! Computation of the respiration term according to pahlow
338                  ! & oschlies (2013)
339                  !
340                  zration = trb(ji,jj,jk,jpndi) / ( trb(ji,jj,jk,jpdia) + rtrn )
341                  zratiop = trb(ji,jj,jk,jppdi) / ( trb(ji,jj,jk,jpdia) + rtrn )
342                  zratiof = trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) + rtrn )
343                  zprnutmax = zprnut(ji,jj,jk) * fvduptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jpdia) * rfact2
344                  ! Uptake of nitrogen
345                  zrat = MIN( 1., zration / (xqndmax(ji,jj,jk) + rtrn) )
346                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05)) 
347                  zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqpdmin(ji,jj,jk) )   &
348                  &          / ( xqpdmax(ji,jj,jk) - xqpdmin(ji,jj,jk) ), xlimdfe(ji,jj,jk) ) )
349                  zpronewd(ji,jj,jk) = zpronmax * zdaylen(ji,jj) * xdiatno3(ji,jj,jk)
350                  zproregd(ji,jj,jk) = zpronmax * xdiatnh4(ji,jj,jk)
351                  ! Uptake of phosphorus
352                  zrat = MIN( 1., zratiop / (xqpdmax(ji,jj,jk) + rtrn) )
353                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05)) 
354                  zpropmax = zprnutmax * zmax * xlimdfe(ji,jj,jk)
355                  zpropo4d(ji,jj,jk) = zpropmax * xdiatpo4(ji,jj,jk)
356                  zprodopd(ji,jj,jk) = zpropmax * xdiatdop(ji,jj,jk)
357                  ! Uptake of iron
358                  zrat = MIN( 1., zratiof / qfdmax )
359                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
360                  zprofmax = zprnutmax * qfdmax * zmax
361                  zprofed(ji,jj,jk) = zprofmax * xdiatfer(ji,jj,jk) * ( 3. - 2.4 * xlimdfe(ji,jj,jk)     &
362                  &          / ( xlimdfe(ji,jj,jk) + 0.2 ) ) * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn   &
363                  &          + xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) ) * (1. - xdiatfer(ji,jj,jk) ) )
364               ENDIF
365            END DO
366         END DO
367      END DO
368
369      DO jk = 1, jpkm1
370         DO jj = 1, jpj
371            DO ji = 1, jpi
372               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
373                     !  production terms for nanophyto. ( chlorophyll )
374                  znanotot = enano(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
375                  zprod = rday * (zpronewn(ji,jj,jk) + zproregn(ji,jj,jk)) * zprchln(ji,jj,jk) * xlimphy(ji,jj,jk)
376                  thetannm_n   = MIN ( thetannm, ( thetannm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
377                  &               * (1. - 1.14 / 43.4 * 20.))
378                  zprochln = thetannm_n * zprod / ( zpislopeadn(ji,jj,jk) * znanotot + rtrn )
379                  zprochln = MAX(zprochln, chlcmin * 12. * zprorcan (ji,jj,jk) )
380                     !  production terms for picophyto. ( chlorophyll )
381                  zpicotot = epico(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
382                  zprod = rday * (zpronewp(ji,jj,jk) + zproregp(ji,jj,jk)) * zprchlp(ji,jj,jk) * xlimpic(ji,jj,jk)
383                  thetanpm_n   = MIN ( thetanpm, ( thetanpm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
384                  &               * (1. - 1.14 / 43.4 * 20.))
385                  zprochlp = thetanpm_n * zprod / ( zpislopeadp(ji,jj,jk) * zpicotot + rtrn )
386                  zprochlp = MAX(zprochlp, chlcmin * 12. * zprorcap(ji,jj,jk) )
387                  !  production terms for diatomees ( chlorophyll )
388                  zdiattot = ediat(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
389                  zprod = rday * (zpronewd(ji,jj,jk) + zproregd(ji,jj,jk)) * zprchld(ji,jj,jk) * xlimdia(ji,jj,jk)
390                  thetandm_n   = MIN ( thetandm, ( thetandm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
391                  &               * (1. - 1.14 / 43.4 * 20.))
392                  zprochld = thetandm_n * zprod / ( zpislopeadd(ji,jj,jk) * zdiattot + rtrn )
393                  zprochld = MAX(zprochld, chlcmin * 12. * zprorcad(ji,jj,jk) )
394                  !   Update the arrays TRA which contain the Chla sources and sinks
395                  tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) + zprochln * texcretn
396                  tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) + zprochld * texcretd
397                  tra(ji,jj,jk,jppch) = tra(ji,jj,jk,jppch) + zprochlp * texcretp
398               ENDIF
399            END DO
400         END DO
401      END DO
402
403      !   Update the arrays TRA which contain the biological sources and sinks
404      DO jk = 1, jpkm1
405         DO jj = 1, jpj
406           DO ji =1 ,jpi
407              zprontot = zpronewn(ji,jj,jk) + zproregn(ji,jj,jk)
408              zproptot = zpronewp(ji,jj,jk) + zproregp(ji,jj,jk)
409              zprodtot = zpronewd(ji,jj,jk) + zproregd(ji,jj,jk)
410              zdocprod = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)  &
411              &          + excretp * zprorcap(ji,jj,jk)
412              tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zpropo4n(ji,jj,jk) - zpropo4d(ji,jj,jk)  &
413              &                     - zpropo4p(ji,jj,jk)
414              tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - zpronewn(ji,jj,jk) - zpronewd(ji,jj,jk)  &
415              &                     - zpronewp(ji,jj,jk)
416              tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zproregn(ji,jj,jk) - zproregd(ji,jj,jk)  &
417              &                     - zproregp(ji,jj,jk)
418              tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zprorcan(ji,jj,jk) * texcretn    &
419                 &                  - zpsino3 * zpronewn(ji,jj,jk) - zpsinh4 * zproregn(ji,jj,jk)   &
420                 &                  - zrespn(ji,jj,jk) 
421              zcroissn(ji,jj,jk) = tra(ji,jj,jk,jpphy) / rfact2/ (trb(ji,jj,jk,jpphy) + rtrn)
422              tra(ji,jj,jk,jpnph) = tra(ji,jj,jk,jpnph) + zprontot * texcretn
423              tra(ji,jj,jk,jppph) = tra(ji,jj,jk,jppph) + zpropo4n(ji,jj,jk) * texcretn   &
424              &                     + zprodopn(ji,jj,jk) * texcretn
425              tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) + zprofen(ji,jj,jk) * texcretn
426              tra(ji,jj,jk,jppic) = tra(ji,jj,jk,jppic) + zprorcap(ji,jj,jk) * texcretp     &
427                 &                  - zpsino3 * zpronewp(ji,jj,jk) - zpsinh4 * zproregp(ji,jj,jk)   &
428                 &                  - zrespp(ji,jj,jk) 
429              zcroissp(ji,jj,jk) = tra(ji,jj,jk,jppic) / rfact2/ (trb(ji,jj,jk,jppic) + rtrn)
430              tra(ji,jj,jk,jpnpi) = tra(ji,jj,jk,jpnpi) + zproptot * texcretp
431              tra(ji,jj,jk,jpppi) = tra(ji,jj,jk,jpppi) + zpropo4p(ji,jj,jk) * texcretp   &
432              &                     + zprodopp(ji,jj,jk) * texcretp
433              tra(ji,jj,jk,jppfe) = tra(ji,jj,jk,jppfe) + zprofep(ji,jj,jk) * texcretp
434              tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorcad(ji,jj,jk) * texcretd   &
435                 &                  - zpsino3 * zpronewd(ji,jj,jk) - zpsinh4 * zproregd(ji,jj,jk)   &
436                 &                  - zrespd(ji,jj,jk) 
437              zcroissd(ji,jj,jk) = tra(ji,jj,jk,jpdia) / rfact2 / (trb(ji,jj,jk,jpdia) + rtrn)
438              tra(ji,jj,jk,jpndi) = tra(ji,jj,jk,jpndi) + zprodtot * texcretd
439              tra(ji,jj,jk,jppdi) = tra(ji,jj,jk,jppdi) + zpropo4d(ji,jj,jk) * texcretd   &
440              &                     + zprodopd(ji,jj,jk) * texcretd
441              tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcretd
442              tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcretd
443              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)  &
444              &                     + excretp * zprorcap(ji,jj,jk)
445              tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + excretd * zprodtot + excretn * zprontot   &
446              &                     + excretp * zproptot
447              tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + excretd * zpropo4d(ji,jj,jk) + excretn * zpropo4n(ji,jj,jk)   &
448              &    - texcretn * zprodopn(ji,jj,jk) - texcretd * zprodopd(ji,jj,jk) + excretp * zpropo4p(ji,jj,jk)     &
449              &    - texcretp * zprodopp(ji,jj,jk)
450              tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproregn(ji,jj,jk) + zproregd(ji,jj,jk)   &
451                 &                + zproregp(ji,jj,jk) ) + ( o2ut + o2nit ) * ( zpronewn(ji,jj,jk)           &
452                 &                + zpronewd(ji,jj,jk) + zpronewp(ji,jj,jk) )   &
453                 &                - o2ut * ( zrespn(ji,jj,jk) + zrespp(ji,jj,jk) + zrespd(ji,jj,jk) )
454              zfeup    = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)  &
455              &          + texcretp * zprofep(ji,jj,jk)
456              tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zfeup
457              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - texcretd * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
458              tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk) - zprorcap(ji,jj,jk)  &
459              &                     + zpsino3 * zpronewn(ji,jj,jk) + zpsinh4 * zproregn(ji,jj,jk)   &
460              &                     + zpsino3 * zpronewp(ji,jj,jk) + zpsinh4 * zproregp(ji,jj,jk)   &
461              &                     + zpsino3 * zpronewd(ji,jj,jk) + zpsinh4 * zproregd(ji,jj,jk)  &
462              &                     + zrespn(ji,jj,jk) + zrespd(ji,jj,jk) + zrespp(ji,jj,jk) 
463              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk)  &
464              &                     + zpronewp(ji,jj,jk) ) - rno3 * ( zproregn(ji,jj,jk) + zproregd(ji,jj,jk)     &
465              &                     + zproregp(ji,jj,jk) ) 
466#if defined key_ligand
467              tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp - zfeup * plig(ji,jj,jk) * lthet
468#endif
469          END DO
470        END DO
471     END DO
472
473     ! Total primary production per year
474     tpp = tpp + glob_sum( ( zprorcan(:,:,:) + zprorcad(:,:,:) + zprorcap(:,:,:) ) * cvol(:,:,:) )
475
476     IF( ln_diatrc ) THEN
477         !
478         zrfact2 = 1.e3 * rfact2r  ! conversion from mol/L/timestep into mol/m3/s
479         IF( lk_iomput ) THEN
480           IF( knt == nrdttrc ) THEN
481              CALL iom_put( "PPPHY"  , zprorcan(:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by nanophyto
482              CALL iom_put( "PPPHYP" , zprorcap(:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by nanophyto
483              CALL iom_put( "PPPHY2" , zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by diatom
484              CALL iom_put( "PPNEWN" , zpronewn(:,:,:) * zrfact2 * tmask(:,:,:) )  ! new primary production by nanophyto
485              CALL iom_put( "PPNEWP" , zpronewp(:,:,:) * zrfact2 * tmask(:,:,:) )  ! new primary production by nanophyto
486              CALL iom_put( "PPNEWD" , zpronewd(:,:,:) * zrfact2 * tmask(:,:,:) )  ! new primary production by diatom
487              CALL iom_put( "PBSi"   , zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:) ) ! biogenic silica production
488              CALL iom_put( "PFeD"   , zprofed (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by diatom
489              CALL iom_put( "PFeP"   , zprofep (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by diatom
490              CALL iom_put( "PFeN"   , zprofen (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by nanophyto
491              CALL iom_put( "Mumax"  , prmaxn(:,:,:) * tmask(:,:,:) )  ! Maximum growth rate
492              CALL iom_put( "MuN"    , zprbio(:,:,:) * xlimphy(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for nanophyto
493              CALL iom_put( "MuP"    , zprpic(:,:,:) * xlimpic(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for nanophyto
494              CALL iom_put( "MuD"    , zprdia(:,:,:) * xlimdia(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for diatoms
495              CALL iom_put( "MunetN"    , zcroissn(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for nanophyto
496              CALL iom_put( "MunetP"    , zcroissp(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for nanophyto
497              CALL iom_put( "MunetD"    , zcroissd(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for diatoms
498              CALL iom_put( "LNlight", zprbio (:,:,:) / (prmaxn(:,:,:) + rtrn) * tmask(:,:,:) )  ! light limitation term
499              CALL iom_put( "LPlight", zprpic (:,:,:) / (prmaxp(:,:,:) + rtrn) * tmask(:,:,:) )  ! light limitation term
500              CALL iom_put( "LDlight", zprdia (:,:,:) / (prmaxd(:,:,:) + rtrn) * tmask(:,:,:) )  ! light limitation term
501           ENDIF
502         ELSE
503              trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorcan(:,:,:) * zrfact2 * tmask(:,:,:)
504              trc3d(:,:,:,jp_pcs0_3d + 5)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:)
505              trc3d(:,:,:,jp_pcs0_3d + 6)  = zpronewn(:,:,:) * zrfact2 * tmask(:,:,:)
506              trc3d(:,:,:,jp_pcs0_3d + 7)  = zpronewd(:,:,:) * zrfact2 * tmask(:,:,:)
507              trc3d(:,:,:,jp_pcs0_3d + 8)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:)
508              trc3d(:,:,:,jp_pcs0_3d + 9)  = zprofed (:,:,:) * zrfact2 * tmask(:,:,:)
509#  if ! defined key_kriest
510              trc3d(:,:,:,jp_pcs0_3d + 10) = zprofen (:,:,:) * zrfact2 * tmask(:,:,:)
511#  endif
512         ENDIF
513         !
514      ENDIF
515
516      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
517         WRITE(charout, FMT="('prod')")
518         CALL prt_ctl_trc_info(charout)
519         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
520      ENDIF
521      !
522      CALL wrk_dealloc( jpi, jpj,      zmixnano, zmixpico, zmixdiat, zstrn )
523      CALL wrk_dealloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
524      CALL wrk_dealloc( jpi, jpj, jpk, zpislopeadn, zpislopeadp, zpislopeadd, zysopt )                           
525      CALL wrk_dealloc( jpi, jpj, jpk, zprdia, zprpic, zprbio, zprorcan, zprorcap, zprorcad )
526      CALL wrk_dealloc( jpi, jpj, jpk, zprofed, zprofep, zprofen ) 
527      CALL wrk_dealloc( jpi, jpj, jpk, zpronewn, zpronewp, zpronewd, zproregn, zproregp, zproregd )
528      CALL wrk_dealloc( jpi, jpj, jpk, zpropo4n, zpropo4p, zpropo4d, zrespn, zrespp, zrespd, zprnut )
529      CALL wrk_dealloc( jpi, jpj, jpk, zprchln, zprchlp, zprchld, zprodopn, zprodopp, zprodopd )
530      CALL wrk_dealloc( jpi, jpj, jpk, zcroissp, zcroissn, zcroissd )
531      !
532      IF( nn_timing == 1 )  CALL timing_stop('p5z_prod')
533      !
534   END SUBROUTINE p5z_prod
535
536
537   SUBROUTINE p5z_prod_init
538      !!----------------------------------------------------------------------
539      !!                  ***  ROUTINE p5z_prod_init  ***
540      !!
541      !! ** Purpose :   Initialization of phytoplankton production parameters
542      !!
543      !! ** Method  :   Read the nampisprod namelist and check the parameters
544      !!      called at the first timestep (nittrc000)
545      !!
546      !! ** input   :   Namelist nampisprod
547      !!----------------------------------------------------------------------
548      !
549      NAMELIST/nampisprod/ pislopen, pislopep, pisloped, xadap, bresp, excretn, excretp, excretd,  &
550         &                 thetannm, thetanpm, thetandm, chlcmin, grosip
551
552      INTEGER :: ios                 ! Local integer output status for namelist read
553      !!----------------------------------------------------------------------
554
555      REWIND( numnatp_ref )              ! Namelist nampisprod in reference namelist : Pisces phytoplankton production
556      READ  ( numnatp_ref, nampisprod, IOSTAT = ios, ERR = 901)
557901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisprod in reference namelist', lwp )
558
559      REWIND( numnatp_cfg )              ! Namelist nampisprod in configuration namelist : Pisces phytoplankton production
560      READ  ( numnatp_cfg, nampisprod, IOSTAT = ios, ERR = 902 )
561902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisprod in configuration namelist', lwp )
562      IF(lwm) WRITE ( numonp, nampisprod )
563
564      IF(lwp) THEN                         ! control print
565         WRITE(numout,*) ' '
566         WRITE(numout,*) ' Namelist parameters for phytoplankton growth, nampisprod'
567         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
568         WRITE(numout,*) '    mean Si/C ratio                           grosip       =', grosip
569         WRITE(numout,*) '    P-I slope                                 pislopen     =', pislopen
570         WRITE(numout,*) '    Acclimation factor to low light           xadap        =', xadap
571         WRITE(numout,*) '    excretion ratio of nanophytoplankton      excretn      =', excretn
572         WRITE(numout,*) '    excretion ratio of picophytoplankton      excretp      =', excretp
573         WRITE(numout,*) '    excretion ratio of diatoms                excretd      =', excretd
574         WRITE(numout,*) '    basal respiration in phytoplankton        bresp        =', bresp
575         WRITE(numout,*) '    Maximum Chl/C in phytoplankton            chlcmin      =', chlcmin
576         WRITE(numout,*) '    P-I slope  for diatoms                    pisloped     =', pisloped
577         WRITE(numout,*) '    P-I slope  for picophytoplankton          pislopep     =', pislopep
578         WRITE(numout,*) '    Minimum Chl/N in nanophytoplankton        thetannm     =', thetannm
579         WRITE(numout,*) '    Minimum Chl/N in picophytoplankton        thetanpm     =', thetanpm
580         WRITE(numout,*) '    Minimum Chl/N in diatoms                  thetandm     =', thetandm
581      ENDIF
582      !
583      r1_rday   = 1._wp / rday 
584      texcretn  = 1._wp - excretn
585      texcretp  = 1._wp - excretp
586      texcretd  = 1._wp - excretd
587      tpp       = 0._wp
588      !
589   END SUBROUTINE p5z_prod_init
590
591
592   INTEGER FUNCTION p5z_prod_alloc()
593      !!----------------------------------------------------------------------
594      !!                     ***  ROUTINE p5z_prod_alloc  ***
595      !!----------------------------------------------------------------------
596      ALLOCATE( prmaxn(jpi,jpj,jpk), prmaxp(jpi,jpj,jpk), prmaxd(jpi,jpj,jpk),  &
597      &          zdaylen(jpi,jpj), STAT = p5z_prod_alloc )
598      !
599      IF( p5z_prod_alloc /= 0 ) CALL ctl_warn('p5z_prod_alloc : failed to allocate arrays.')
600      !
601   END FUNCTION p5z_prod_alloc
602
603#else
604   !!======================================================================
605   !!  Dummy module :                                   No PISCES bio-model
606   !!======================================================================
607CONTAINS
608   SUBROUTINE p5z_prod                    ! Empty routine
609   END SUBROUTINE p5z_prod
610#endif 
611
612   !!======================================================================
613END MODULE  p5zprod
Note: See TracBrowser for help on using the repository browser.