source: NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90 @ 12759

Last change on this file since 12759 was 12759, checked in by aumont, 9 months ago

make parameterizations in PISCES-operationnal more similar to thos of PISCES-QUOTA (prey switching, optimal allocation, size, …)

  • Property svn:keywords set to Id
File size: 31.5 KB
Line 
1MODULE p4zprod
2   !!======================================================================
3   !!                         ***  MODULE p4zprod  ***
4   !! TOP :  Growth Rate of the two phytoplankton groups of PISCES
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   !!----------------------------------------------------------------------
10   !!   p4z_prod       : Compute the growth Rate of the two phytoplanktons groups
11   !!   p4z_prod_init  : Initialization of the parameters for growth
12   !!   p4z_prod_alloc : Allocate variables for growth
13   !!----------------------------------------------------------------------
14   USE oce_trc         ! shared variables between ocean and passive tracers
15   USE trc             ! passive tracers common variables
16   USE sms_pisces      ! PISCES Source Minus Sink variables
17   USE p4zlim          ! Co-limitations of differents nutrients
18   USE prtctl_trc      ! print control for debugging
19   USE iom             ! I/O manager
20
21   IMPLICIT NONE
22   PRIVATE
23
24   PUBLIC   p4z_prod         ! called in p4zbio.F90
25   PUBLIC   p4z_prod_init    ! called in trcsms_pisces.F90
26   PUBLIC   p4z_prod_alloc   ! called in trcini_pisces.F90
27
28   REAL(wp), PUBLIC ::   pislopen     !:
29   REAL(wp), PUBLIC ::   pisloped     !:
30   REAL(wp), PUBLIC ::   xadap        !:
31   REAL(wp), PUBLIC ::   excretn      !:
32   REAL(wp), PUBLIC ::   excretd      !:
33   REAL(wp), PUBLIC ::   bresp        !:
34   REAL(wp), PUBLIC ::   chlcnm       !:
35   REAL(wp), PUBLIC ::   chlcdm       !:
36   REAL(wp), PUBLIC ::   chlcmin      !:
37   REAL(wp), PUBLIC ::   fecnm        !:
38   REAL(wp), PUBLIC ::   fecdm        !:
39   REAL(wp), PUBLIC ::   grosip       !:
40
41   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotan   !: proxy of N quota in Nanophyto
42   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotad   !: proxy of N quota in diatoms
43   
44   REAL(wp) ::   r1_rday    ! 1 / rday
45   REAL(wp) ::   texcretn   ! 1 - excretn
46   REAL(wp) ::   texcretd   ! 1 - excretd       
47
48   !!----------------------------------------------------------------------
49   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
50   !! $Id$
51   !! Software governed by the CeCILL license (see ./LICENSE)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE p4z_prod( kt , knt )
56      !!---------------------------------------------------------------------
57      !!                     ***  ROUTINE p4z_prod  ***
58      !!
59      !! ** Purpose :   Computes phytoplankton production depending on
60      !!                light, temperature and nutrient availability
61      !!                Computes also the uptake of Iron and Si as well
62      !!                as the chlorophyll content of the cells
63      !!---------------------------------------------------------------------
64      INTEGER, INTENT(in) ::   kt, knt   !
65      !
66      INTEGER  ::   ji, jj, jk
67      REAL(wp) ::   zsilfac, znanotot, zdiattot, zconctemp, zconctemp2
68      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap, zlim, zsilfac2, zsiborn
69      REAL(wp) ::   zprod, zproreg, zproreg2, zprochln, zprochld
70      REAL(wp) ::   zdocprod, zpislopen, zpisloped, zfact
71      REAL(wp) ::   zratiosi, zmaxsi, zlimfac, zsizetmp
72      REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup, chlcnm_n, chlcdm_n
73      CHARACTER (len=25) :: charout
74      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zw2d
75      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
76      REAL(wp), DIMENSION(jpi,jpj    ) :: zstrn
77      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprmaxn,zprmaxd
78      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadd, zysopt 
79      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprdia, zprbio, zprchld, zprchln   
80      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprorcan, zprorcad, zprofed, zprofen
81      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpronewn, zpronewd
82      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmxl_fac, zmxl_chl
83      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpligprod1, zpligprod2
84      !!---------------------------------------------------------------------
85      !
86      IF( ln_timing )   CALL timing_start('p4z_prod')
87      !
88      !  Allocate temporary workspace
89      !
90      zprorcan(:,:,:) = 0._wp ; zprorcad(:,:,:) = 0._wp ; zprofed (:,:,:) = 0._wp
91      zprofen (:,:,:) = 0._wp ; zysopt  (:,:,:) = 0._wp
92      zpronewn(:,:,:) = 0._wp ; zpronewd(:,:,:) = 0._wp ; zprdia  (:,:,:) = 0._wp
93      zprbio  (:,:,:) = 0._wp ; zprchld (:,:,:) = 0._wp ; zprchln (:,:,:) = 0._wp 
94      zmxl_fac(:,:,:) = 0._wp ; zmxl_chl(:,:,:) = 0._wp 
95
96      ! Computation of the maximimum production
97      ! Parameters are taken from Bissinger et al. (2008)
98      zprmaxn(:,:,:) = 0.8_wp * r1_rday * tgfunc(:,:,:)
99      zprmaxd(:,:,:) = zprmaxn(:,:,:)
100
101      ! compute the day length depending on latitude and the day
102      zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
103      zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
104
105      ! day length in hours
106      zstrn(:,:) = 0.
107      DO jj = 1, jpj
108         DO ji = 1, jpi
109            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
110            zargu = MAX( -1., MIN(  1., zargu ) )
111            zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
112         END DO
113      END DO
114
115      ! Impact of the day duration and light intermittency on phytoplankton growth
116      ! Intermittency is supposed to have a similar effect on production as
117      ! day length. The correcting factor is zmxl_fac. zmxl_chl is the fractional
118      ! day length and is used to compute the mean PAR during daytime.
119      ! Formulation for the impact of day length on PP is from Thompson (1999)
120      ! -------------------------------------------------------------------------
121      DO jk = 1, jpkm1
122         DO jj = 1 ,jpj
123            DO ji = 1, jpi
124               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
125                  zval = MAX( 1., zstrn(ji,jj) )
126                  IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
127                     zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
128                  ENDIF
129                  zmxl_chl(ji,jj,jk) = zval / 24.
130                  zmxl_fac(ji,jj,jk) = 1.0 - exp( -0.26 * zval )
131               ENDIF
132            END DO
133         END DO
134      END DO
135
136      zprbio(:,:,:) = zprmaxn(:,:,:) * zmxl_fac(:,:,:)
137      zprdia(:,:,:) = zprmaxd(:,:,:) * zmxl_fac(:,:,:)
138
139      WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
140
141      ! Computation of the P-I slope for nanos and diatoms
142      ! The formulation proposed by Geider et al. (1997) has been modified
143      ! to exclude the effect of nutrient limitation and temperature in the PI
144      ! curve following Vichi et al. (2007)
145      ! -----------------------------------------------------------------------
146      DO jk = 1, jpkm1
147         DO jj = 1, jpj
148            DO ji = 1, jpi
149               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
150                  ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
151                  zadap       = xadap * ztn / ( 2.+ ztn )
152                  zconctemp   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
153                  zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp
154
155                  ! The initial slope of the PI curve can be increased for nano
156                  ! to account for photadaptation, for instance in the DCM
157                  zpislopeadn(ji,jj,jk) = pislopen * ( 1.+ zadap  * EXP( -0.25 * enano(ji,jj,jk) ) )  &
158                  &                   * trb(ji,jj,jk,jpnch) /( trb(ji,jj,jk,jpphy) * 12. + rtrn)
159                  !
160                  zpislopeadd(ji,jj,jk) = (pislopen * zconctemp2 + pisloped * zconctemp) / ( trb(ji,jj,jk,jpdia) + rtrn )   &
161                  &                   * trb(ji,jj,jk,jpdch) /( trb(ji,jj,jk,jpdia) * 12. + rtrn)
162               ENDIF
163            END DO
164         END DO
165      END DO
166
167      DO jk = 1, jpkm1
168         DO jj = 1, jpj
169            DO ji = 1, jpi
170               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
171                   ! Computation of production function for Carbon
172                   !  ---------------------------------------------
173                   zpislopen = zpislopeadn(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
174                   &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
175                   zpisloped = zpislopeadd(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
176                   &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
177                   zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
178                   zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
179                   !  Computation of production function for Chlorophyll
180                   !--------------------------------------------------
181                   zpislopen = zpislopeadn(ji,jj,jk) / ( zprmaxn(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
182                   zpisloped = zpislopeadd(ji,jj,jk) / ( zprmaxd(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
183                   zprchln(ji,jj,jk) = zprmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) ) )
184                   zprchld(ji,jj,jk) = zprmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) ) )
185               ENDIF
186            END DO
187         END DO
188      END DO
189
190      !  Computation of a proxy of the N/C ratio
191      !  Steady state is assumed
192      !  ---------------------------------------
193      DO jk = 1, jpkm1
194         DO jj = 1, jpj
195            DO ji = 1, jpi
196                zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) )   &
197                &      * zprmaxn(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn )
198                quotan(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
199                zval = MIN( xdiatpo4(ji,jj,jk), ( xdiatnh4(ji,jj,jk) + xdiatno3(ji,jj,jk) ) )   &
200                &      * zprmaxd(ji,jj,jk) / ( zprdia(ji,jj,jk) + rtrn )
201                quotad(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
202            END DO
203         END DO
204      END DO
205
206
207      DO jk = 1, jpkm1
208         DO jj = 1, jpj
209            DO ji = 1, jpi
210
211               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
212                   ! Si/C of diatoms
213                   ! ------------------------
214                   ! Si/C increases with iron stress and silicate availability (zsilfac)
215                   ! Si/C is arbitrariliy increased for very high Si concentrations
216                   ! to mimic the very high ratios observed in the Southern Ocean (zsilfac2)
217                   ! -----------------------------------------------------------------------
218                 zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
219                 zsilim = xlimdia(ji,jj,jk) * zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn )
220                 zsiborn = trb(ji,jj,1,jpsil) * trb(ji,jj,1,jpsil) * trb(ji,jj,1,jpsil)
221                 IF (gphit(ji,jj) < -30 ) THEN
222                   zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
223                 ELSE
224                   zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 )
225                 ENDIF
226                 zratiosi = 1.0 - trb(ji,jj,jk,jpdsi) / ( trb(ji,jj,jk,jpdia) + rtrn ) / ( zsilfac2 * grosip * 3.0 + rtrn )
227                 zratiosi = MAX(0., MIN(1.0, zratiosi) )
228                 zmaxsi  = (1.0 + 0.1**4) * zratiosi**4 / ( zratiosi**4 + 0.1**4 )
229                 IF ( xlimsi(ji,jj,jk) /= xlimdia(ji,jj,jk) ) THEN
230                    zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zmaxsi
231                 ELSE
232                    zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zsilim**0.7 * zmaxsi
233                 ENDIF
234              ENDIF
235            END DO
236         END DO
237      END DO
238
239      ! Sea-ice effect on production
240      ! No production is assumed below sea ice
241      ! --------------------------------------
242      DO jk = 1, jpkm1
243         DO jj = 1, jpj
244            DO ji = 1, jpi
245               zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
246               zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
247            END DO
248         END DO
249      END DO
250
251      ! Computation of the various production  and nutrient uptake terms
252      ! ---------------------------------------------------------------
253      DO jk = 1, jpkm1
254         DO jj = 1, jpj
255            DO ji = 1, jpi
256               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
257                  !  production term of nanophyto. (C)
258                  zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trb(ji,jj,jk,jpphy) * rfact2
259
260                  !  New production (uptake of NO3)
261                  zpronewn(ji,jj,jk)  = zprorcan(ji,jj,jk)* xnanono3(ji,jj,jk) / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn )
262
263                  ! Size computation
264                  ! Size is made a function of the limitation of of phytoplankton growth
265                  ! Strongly limited cells are supposed to be smaller. sizena is the
266                  ! size at time step t+1 and is thus updated at the end of the
267                  ! current time step
268                  ! --------------------------------------------------------------------
269                  zlimfac = xlimphy(ji,jj,jk) * zprchln(ji,jj,jk) / ( zprmaxn(ji,jj,jk) + rtrn )
270                  zsizetmp = 1.0 + 1.3 * ( xsizern - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3)
271                  sizena(ji,jj,jk) = min(xsizern, max( sizena(ji,jj,jk), zsizetmp ) )
272
273                  !  Iron uptake rates of nanophytoplankton. Upregulation 
274                  !  is parameterized at low iron concentrations. Typical
275                  !  formulation used in quota formulations. Uptake is downregulated
276                  !  when the quota is close to the maximum quota
277                  zratio = 1.0 - MIN(1.0,trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * fecnm + rtrn ) )
278                  zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) ) 
279                  zprofen(ji,jj,jk) = fecnm * zprmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
280                  &          * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn + xnanono3(ji,jj,jk)  &
281                  &          + xnanonh4(ji,jj,jk) ) * (1. - xnanofer(ji,jj,jk) ) )   &
282                  &          * xnanofer(ji,jj,jk) * zmax * trb(ji,jj,jk,jpphy) * rfact2
283                  !  production terms of diatoms (C)
284                  zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trb(ji,jj,jk,jpdia) * rfact2
285
286                  ! New production (uptake of NO3)
287                  zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk) / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn )
288
289                  ! Size computation
290                  ! Size is made a function of the limitation of of phytoplankton growth
291                  ! Strongly limited cells are supposed to be smaller. sizeda is
292                  ! size at time step t+1 and is thus updated at the end of the
293                  ! current time step.
294                  ! --------------------------------------------------------------------
295                  zlimfac = zprchld(ji,jj,jk) * xlimdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn )
296                  zsizetmp = 1.0 + 1.3 * ( xsizerd - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3)
297                  sizeda(ji,jj,jk) = min(xsizerd, max( sizeda(ji,jj,jk), zsizetmp ) )
298
299                  !  Iron uptake rates of nanophytoplankton. Upregulation 
300                  !  is parameterized at low iron concentrations. Typical
301                  !  formulation used in quota formulations. Uptake is downregulated
302                  !  when the quota is close to the maximum quota
303                  zratio = 1.0 - MIN(1.0, trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * fecdm + rtrn ) )
304                  zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) ) 
305                  zprofed(ji,jj,jk) = fecdm * zprmaxd(ji,jj,jk) * (1.0 - fr_i(ji,jj) )  &
306                  &          * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn + xdiatno3(ji,jj,jk)  &
307                  &          + xdiatnh4(ji,jj,jk) ) * (1. - xdiatfer(ji,jj,jk) ) )   &
308                  &          * xdiatfer(ji,jj,jk) * zmax * trb(ji,jj,jk,jpdia) * rfact2
309               ENDIF
310            END DO
311         END DO
312      END DO
313
314      ! Computation of the chlorophyll production terms
315      ! The parameterization is taken from Geider et al. (1997)
316      ! -------------------------------------------------------
317      DO jk = 1, jpkm1
318         DO jj = 1, jpj
319            DO ji = 1, jpi
320               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
321                  !  production term for nanophyto. ( chlorophyll )
322                  znanotot = enanom(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
323                  zprod    = rday * zprorcan(ji,jj,jk) * zprchln(ji,jj,jk) * xlimphy(ji,jj,jk)
324                  zprochln = chlcmin * 12. * zprorcan (ji,jj,jk)
325                  ! The maximum reachable Chl quota is modulated by temperature
326                  ! following Geider (1987)
327                  chlcnm_n   = MIN ( chlcnm, ( chlcnm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem))) * (1. - 1.14 / 43.4 * 20.))
328                  zprochln = zprochln + (chlcnm_n-chlcmin) * 12. * zprod / &
329                                        & (  zpislopeadn(ji,jj,jk) * znanotot +rtrn)
330                  !  production terms of diatoms ( chlorophyll )
331                  zdiattot = ediatm(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
332                  zprod    = rday * zprorcad(ji,jj,jk) * zprchld(ji,jj,jk) * xlimdia(ji,jj,jk)
333                  zprochld = chlcmin * 12. * zprorcad(ji,jj,jk)
334                  ! The maximum reachable Chl quota is modulated by temperature
335                  ! following Geider (1987)
336                  chlcdm_n   = MIN ( chlcdm, ( chlcdm / (1. - 1.14 / 43.4 * tsn(ji,jj,jk,jp_tem))) * (1. - 1.14 / 43.4 * 20.))
337                  zprochld = zprochld + (chlcdm_n-chlcmin) * 12. * zprod / &
338                                        & ( zpislopeadd(ji,jj,jk) * zdiattot +rtrn )
339                  !   Update the arrays TRA which contain the Chla sources and sinks
340                  tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) + zprochln * texcretn
341                  tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) + zprochld * texcretd
342               ENDIF
343            END DO
344         END DO
345      END DO
346
347      !   Update the arrays TRA which contain the biological sources and sinks
348      DO jk = 1, jpkm1
349         DO jj = 1, jpj
350           DO ji =1 ,jpi
351              IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
352                 zproreg  = zprorcan(ji,jj,jk) - zpronewn(ji,jj,jk)
353                 zproreg2 = zprorcad(ji,jj,jk) - zpronewd(ji,jj,jk)
354                 zdocprod = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)
355                 tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk)
356                 tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - zpronewn(ji,jj,jk) - zpronewd(ji,jj,jk)
357                 tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zproreg - zproreg2
358                 tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zprorcan(ji,jj,jk) * texcretn
359                 tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) + zprofen(ji,jj,jk) * texcretn
360                 tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorcad(ji,jj,jk) * texcretd
361                 tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcretd
362                 tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk)   &
363                 &                   * rfact2 * trb(ji,jj,jk,jpdia)
364                 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zdocprod
365                 tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) &
366                 &                   + ( o2ut + o2nit ) * ( zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk) )
367                 !
368                 zfeup = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)
369                 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zfeup
370                 tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk)   &
371                 &                   * rfact2 * trb(ji,jj,jk,jpdia)
372                 tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk)
373                 tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk) ) &
374                 &                                         - rno3 * ( zproreg + zproreg2 )
375              ENDIF
376           END DO
377        END DO
378     END DO
379     !
380     ! Production and uptake of ligands by phytoplankton. This part is activated
381     ! when ln_ligand is set to .true. in the namelist. Ligand uptake is small
382     ! and based on the FeL model by Morel et al. (2008) and on the study of
383     ! Shaked and Lis (2012)
384     ! -------------------------------------------------------------------------
385     IF( ln_ligand ) THEN
386         zpligprod1(:,:,:) = 0._wp    ;    zpligprod2(:,:,:) = 0._wp
387         DO jk = 1, jpkm1
388            DO jj = 1, jpj
389              DO ji =1 ,jpi
390                 IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
391                    zdocprod = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)
392                    zfeup    = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)
393                    tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp    &
394                    &       - zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
395                    zpligprod1(ji,jj,jk) = zdocprod * ldocp
396                    zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) &
397                    &                      + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
398                 ENDIF
399              END DO
400           END DO
401        END DO
402     ENDIF
403
404
405    ! Total primary production per year
406    IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
407         & tpp = glob_sum( 'p4zprod', ( zprorcan(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
408
409    IF( lk_iomput ) THEN
410       IF( knt == nrdttrc ) THEN
411          ALLOCATE( zw2d(jpi,jpj), zw3d(jpi,jpj,jpk) )
412          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
413          !
414          IF( iom_use( "PPPHYN" ) .OR. iom_use( "PPPHYD" ) )  THEN
415              zw3d(:,:,:) = zprorcan(:,:,:) * zfact * tmask(:,:,:)  ! primary production by nanophyto
416              CALL iom_put( "PPPHYN"  , zw3d )
417              !
418              zw3d(:,:,:) = zprorcad(:,:,:) * zfact * tmask(:,:,:)  ! primary production by diatomes
419              CALL iom_put( "PPPHYD"  , zw3d )
420          ENDIF
421          IF( iom_use( "PPNEWN" ) .OR. iom_use( "PPNEWD" ) )  THEN
422              zw3d(:,:,:) = zpronewn(:,:,:) * zfact * tmask(:,:,:)  ! new primary production by nanophyto
423              CALL iom_put( "PPNEWN"  , zw3d )
424              !
425              zw3d(:,:,:) = zpronewd(:,:,:) * zfact * tmask(:,:,:)  ! new primary production by diatomes
426              CALL iom_put( "PPNEWD"  , zw3d )
427          ENDIF
428          IF( iom_use( "PBSi" ) )  THEN
429              zw3d(:,:,:) = zprmaxd(:,:,:) * 1.E3 * tmask(:,:,:) * zysopt(:,:,:) * trb(:,:,:,jpdia) ! biogenic silica production
430              CALL iom_put( "PBSi"  , zw3d )
431          ENDIF
432          IF( iom_use( "PFeN" ) .OR. iom_use( "PFeD" ) )  THEN
433              zw3d(:,:,:) = zprofen(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron production by nanophyto
434              CALL iom_put( "PFeN"  , zw3d )
435              !
436              zw3d(:,:,:) = zprofed(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron production by  diatomes
437              CALL iom_put( "PFeD"  , zw3d )
438          ENDIF
439          IF( iom_use( "LPRODP" ) )  THEN
440              zw3d(:,:,:) = zpligprod1(:,:,:) * 1e9 * zfact * tmask(:,:,:)  ! Ligand production by phytoplankton
441              CALL iom_put( "LPRODP"  , zw3d )
442          ENDIF
443          IF( iom_use( "LDETP" ) )  THEN
444              zw3d(:,:,:) = zpligprod2(:,:,:) * 1e9 * zfact * tmask(:,:,:)  ! Uptake of ligands by phytoplankton
445              CALL iom_put( "LDETP"  , zw3d )
446          ENDIF
447          IF( iom_use( "Mumax" ) )  THEN
448              zw3d(:,:,:) = zprmaxn(:,:,:) * tmask(:,:,:)   ! Maximum growth rate
449              CALL iom_put( "Mumax"  , zw3d )
450          ENDIF
451          IF( iom_use( "MuN" ) .OR. iom_use( "MuD" ) )  THEN
452              zw3d(:,:,:) = zprbio(:,:,:) * xlimphy(:,:,:) * tmask(:,:,:)  ! Realized growth rate for nanophyto
453              CALL iom_put( "MuN"  , zw3d )
454              !
455              zw3d(:,:,:) =  zprdia(:,:,:) * xlimdia(:,:,:) * tmask(:,:,:)  ! Realized growth rate for diatoms
456              CALL iom_put( "MuD"  , zw3d )
457          ENDIF
458          IF( iom_use( "LNlight" ) .OR. iom_use( "LDlight" ) )  THEN
459              zw3d(:,:,:) = zprbio (:,:,:) / (zprmaxn(:,:,:) + rtrn) * tmask(:,:,:) ! light limitation term of nanophytoplankton
460              CALL iom_put( "LNlight"  , zw3d )
461              !
462              zw3d(:,:,:) = zprdia (:,:,:) / (zprmaxd(:,:,:) + rtrn) * tmask(:,:,:)  ! light limitation term of diatoms
463              CALL iom_put( "LDlight"  , zw3d )
464          ENDIF
465          IF( iom_use( "TPP" ) )  THEN
466              zw3d(:,:,:) = ( zprorcan(:,:,:) + zprorcad(:,:,:) ) * zfact * tmask(:,:,:)  ! total primary production
467              CALL iom_put( "TPP"  , zw3d )
468          ENDIF
469          IF( iom_use( "TPNEW" ) )  THEN
470              zw3d(:,:,:) = ( zpronewn(:,:,:) + zpronewd(:,:,:) ) * zfact * tmask(:,:,:)  ! total new production
471              CALL iom_put( "TPNEW"  , zw3d )
472          ENDIF
473          IF( iom_use( "TPBFE" ) )  THEN
474              zw3d(:,:,:) = ( zprofen(:,:,:) + zprofed(:,:,:) ) * zfact * tmask(:,:,:)  ! total biogenic iron production
475              CALL iom_put( "TPBFE"  , zw3d )
476          ENDIF
477          IF( iom_use( "INTPPPHYN" ) .OR. iom_use( "INTPPPHYD" ) ) THEN 
478             zw2d(:,:) = 0.
479             DO jk = 1, jpkm1
480               zw2d(:,:) = zw2d(:,:) + zprorcan(:,:,jk) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated  primary produc. by nano
481             ENDDO
482             CALL iom_put( "INTPPPHYN" , zw2d )
483             !
484             zw2d(:,:) = 0.
485             DO jk = 1, jpkm1
486                zw2d(:,:) = zw2d(:,:) + zprorcad(:,:,jk) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk) ! vert. integrated  primary produc. by diatom
487             ENDDO
488             CALL iom_put( "INTPPPHYD" , zw2d )
489          ENDIF
490          IF( iom_use( "INTPP" ) ) THEN   
491             zw2d(:,:) = 0.
492             DO jk = 1, jpkm1
493                zw2d(:,:) = zw2d(:,:) + ( zprorcan(:,:,jk) + zprorcad(:,:,jk) ) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk) ! vert. integrated pp
494             ENDDO
495             CALL iom_put( "INTPP" , zw2d )
496          ENDIF
497          IF( iom_use( "INTPNEW" ) ) THEN   
498             zw2d(:,:) = 0.
499             DO jk = 1, jpkm1
500                zw2d(:,:) = zw2d(:,:) + ( zpronewn(:,:,jk) + zpronewd(:,:,jk) ) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated new prod
501             ENDDO
502             CALL iom_put( "INTPNEW" , zw2d )
503          ENDIF
504          IF( iom_use( "INTPBFE" ) ) THEN           !   total biogenic iron production  ( vertically integrated )
505             zw2d(:,:) = 0.
506             DO jk = 1, jpkm1
507                zw2d(:,:) = zw2d(:,:) + ( zprofen(:,:,jk) + zprofed(:,:,jk) ) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk) ! vert integr. bfe prod
508             ENDDO
509            CALL iom_put( "INTPBFE" , zw2d )
510          ENDIF
511          IF( iom_use( "INTPBSI" ) ) THEN           !   total biogenic silica production  ( vertically integrated )
512             zw2d(:,:) = 0.
513             DO jk = 1, jpkm1
514                zw2d(:,:) = zw2d(:,:) + zprorcad(:,:,jk) * zysopt(:,:,jk) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk)  ! vert integr. bsi prod
515             ENDDO
516             CALL iom_put( "INTPBSI" , zw2d )
517          ENDIF
518          IF( iom_use( "tintpp" ) )  CALL iom_put( "tintpp" , tpp * zfact )  !  global total integrated primary production molC/s
519          !
520          DEALLOCATE( zw2d, zw3d )
521       ENDIF
522     ENDIF
523
524     IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
525         WRITE(charout, FMT="('prod')")
526         CALL prt_ctl_trc_info(charout)
527         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
528     ENDIF
529      !
530      IF( ln_timing )  CALL timing_stop('p4z_prod')
531      !
532   END SUBROUTINE p4z_prod
533
534
535   SUBROUTINE p4z_prod_init
536      !!----------------------------------------------------------------------
537      !!                  ***  ROUTINE p4z_prod_init  ***
538      !!
539      !! ** Purpose :   Initialization of phytoplankton production parameters
540      !!
541      !! ** Method  :   Read the namp4zprod namelist and check the parameters
542      !!      called at the first timestep (nittrc000)
543      !!
544      !! ** input   :   Namelist namp4zprod
545      !!----------------------------------------------------------------------
546      INTEGER ::   ios   ! Local integer
547      !
548      NAMELIST/namp4zprod/ pislopen, pisloped, xadap, bresp, excretn, excretd,  &
549         &                 chlcnm, chlcdm, chlcmin, fecnm, fecdm, grosip
550      !!----------------------------------------------------------------------
551      !
552      IF(lwp) THEN                         ! control print
553         WRITE(numout,*)
554         WRITE(numout,*) 'p4z_prod_init : phytoplankton growth'
555         WRITE(numout,*) '~~~~~~~~~~~~~'
556      ENDIF
557      !
558      REWIND( numnatp_ref )              ! Namelist namp4zprod in reference namelist : Pisces phytoplankton production
559      READ  ( numnatp_ref, namp4zprod, IOSTAT = ios, ERR = 901)
560901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zprod in reference namelist' )
561      REWIND( numnatp_cfg )              ! Namelist namp4zprod in configuration namelist : Pisces phytoplankton production
562      READ  ( numnatp_cfg, namp4zprod, IOSTAT = ios, ERR = 902 )
563902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zprod in configuration namelist' )
564      IF(lwm) WRITE( numonp, namp4zprod )
565
566      IF(lwp) THEN                         ! control print
567         WRITE(numout,*) '   Namelist : namp4zprod'
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 diatoms                excretd      =', excretd
573         WRITE(numout,*) '      basal respiration in phytoplankton        bresp        =', bresp
574         WRITE(numout,*) '      Maximum Chl/C in phytoplankton            chlcmin      =', chlcmin
575         WRITE(numout,*) '      P-I slope  for diatoms                    pisloped     =', pisloped
576         WRITE(numout,*) '      Minimum Chl/C in nanophytoplankton        chlcnm       =', chlcnm
577         WRITE(numout,*) '      Minimum Chl/C in diatoms                  chlcdm       =', chlcdm
578         WRITE(numout,*) '      Maximum Fe/C in nanophytoplankton         fecnm        =', fecnm
579         WRITE(numout,*) '      Minimum Fe/C in diatoms                   fecdm        =', fecdm
580      ENDIF
581      !
582      r1_rday   = 1._wp / rday 
583      texcretn  = 1._wp - excretn
584      texcretd  = 1._wp - excretd
585      tpp       = 0._wp
586      !
587   END SUBROUTINE p4z_prod_init
588
589
590   INTEGER FUNCTION p4z_prod_alloc()
591      !!----------------------------------------------------------------------
592      !!                     ***  ROUTINE p4z_prod_alloc  ***
593      !!----------------------------------------------------------------------
594      ALLOCATE( quotan(jpi,jpj,jpk), quotad(jpi,jpj,jpk), STAT = p4z_prod_alloc )
595      !
596      IF( p4z_prod_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_prod_alloc : failed to allocate arrays.' )
597      !
598   END FUNCTION p4z_prod_alloc
599
600   !!======================================================================
601END MODULE p4zprod
Note: See TracBrowser for help on using the repository browser.