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.
p4zprod.F90 in NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z – NEMO

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

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

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