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 trunk/NEMO/TOP_SRC/PISCES – NEMO

source: trunk/NEMO/TOP_SRC/PISCES/p4zprod.F90 @ 1284

Last change on this file since 1284 was 1263, checked in by cetlod, 15 years ago

add missing USE lib_mpp in p4zprod routine, see ticket:302
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 18.2 KB
Line 
1MODULE p4zprod
2   !!======================================================================
3   !!                         ***  MODULE p4zprod  ***
4   !! TOP :   PISCES
5   !!======================================================================
6   !! History :   1.0  !  2004     (O. Aumont) Original code
7   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
8   !!----------------------------------------------------------------------
9#if defined key_pisces
10   !!----------------------------------------------------------------------
11   !!   'key_pisces'                                       PISCES bio-model
12   !!----------------------------------------------------------------------
13   !!   p4z_prod       : 
14   !!----------------------------------------------------------------------
15   USE trc
16   USE oce_trc         !
17   USE sms_pisces      !
18   USE prtctl_trc
19   USE p4zopt
20   USE p4zint
21   USE p4zlim
22
23   USE lib_mpp
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   p4z_prod    ! called in p4zbio.F90
29
30   !! * Shared module variables
31   REAL(wp), PUBLIC ::   &
32     pislope   = 3.0_wp          ,  &  !:
33     pislope2  = 3.0_wp          ,  &  !:
34     excret    = 10.e-5_wp       , &   !:
35     excret2   = 0.05_wp         , &   !:
36     chlcnm    = 0.033_wp        , &   !:
37     chlcdm    = 0.05_wp         , &   !:
38     fecnm     = 10.E-6_wp       , &   !:
39     fecdm     = 15.E-6_wp       , &   !:
40     grosip    = 0.151_wp
41
42   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk)  ::        &
43     &                   prmax
44   
45   REAL(wp) ::   &
46      texcret                    ,  &  !: 1 - excret
47      texcret2                   ,  &  !: 1 - excret2       
48      rpis180                    ,  &  !: rpi / 180
49      tpp = 0.                         !: Total primary production
50
51   !!* Substitution
52#  include "domzgr_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
55   !! $Id$
56   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58
59CONTAINS
60
61   SUBROUTINE p4z_prod( kt , jnt )
62      !!---------------------------------------------------------------------
63      !!                     ***  ROUTINE p4z_prod  ***
64      !!
65      !! ** Purpose :   Compute the phytoplankton production depending on
66      !!              light, temperature and nutrient availability
67      !!
68      !! ** Method  : - ???
69      !!---------------------------------------------------------------------
70      INTEGER, INTENT(in) :: kt, jnt
71      INTEGER  ::   ji, jj, jk, nspyr
72      REAL(wp) ::   zsilfac, zfact, zrfact2
73      REAL(wp) ::   zprdiachl, zprbiochl, zsilim, ztn, zadap, zadap2
74      REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zetot2, zmax, zproreg, zproreg2
75      REAL(wp) ::   zmxltst, zmxlday, zlim1
76      REAL(wp) ::   zpislopen  , zpislope2n
77      REAL(wp) ::   zrum, zcodel, zargu, zvol
78      REAL(wp), DIMENSION(jpi,jpj)     ::   zmixnano   , zmixdiat, zstrn
79      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpislopead , zpislopead2
80      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprdia     , zprbio, zysopt
81      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprorca    , zprorcad, zprofed
82      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprofen   , zprochln, zprochld
83      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpronew    , zpronewd
84      CHARACTER (len=25) :: charout
85      !!---------------------------------------------------------------------
86
87
88      IF( ( kt * jnt ) == nittrc000  )   CALL p4z_prod_init      ! Initialization (first time-step only)
89
90
91      zprorca (:,:,:) = 0.0
92      zprorcad(:,:,:) = 0.0
93      zprofed(:,:,:) = 0.0
94      zprofen(:,:,:) = 0.0
95      zprochln(:,:,:) = 0.0
96      zprochld(:,:,:) = 0.0
97      zpronew (:,:,:) = 0.0
98      zpronewd(:,:,:) = 0.0
99      zprdia  (:,:,:) = 0.0
100      zprbio  (:,:,:) = 0.0
101      zysopt  (:,:,:) = 0.0
102
103      nspyr  = INT( raass / rdt )
104
105
106!     Computation of the optimal production
107!     -------------------------------------
108
109# if defined key_off_degrad
110      prmax(:,:,:) = 0.6 / rjjss * tgfunc(:,:,:) * facvol(:,:,:)
111# else
112      prmax(:,:,:) = 0.6 / rjjss * tgfunc(:,:,:)
113# endif
114
115      ! compute the day length depending on latitude and the day
116      !--------------------------------------------------------
117      IF(lwp) write(numout,*)
118      IF(lwp) write(numout,*) 'p4zday : - Julian day ', nday_year
119      IF(lwp) write(numout,*) '~~~~~~'
120
121      IF( nleapy == 1 .AND. MOD( nyear, 4 ) == 0 ) THEN
122         zrum = FLOAT( nday_year - 80 ) / 366.
123      ELSE
124         zrum = FLOAT( nday_year - 80 ) / 365.
125      ENDIF
126      zcodel = ASIN(  SIN( zrum * rpi * 2. ) * SIN( rpis180 * 23.5 )  )
127
128      ! day length in hours
129      zstrn(:,:) = 0.
130      DO jj = 1, jpj
131         DO ji = 1, jpi
132            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rpis180 )
133            zargu = MAX( -1., MIN(  1., zargu ) )
134            zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rpis180 / 15. )
135         END DO
136      END DO
137
138
139!CDIR NOVERRCHK
140      DO jk = 1, jpkm1
141!CDIR NOVERRCHK
142         DO jj = 1, jpj
143!CDIR NOVERRCHK
144            DO ji = 1, jpi
145
146!      Computation of the P-I slope for nanos and diatoms
147!      --------------------------------------------------
148                IF( etot(ji,jj,jk) > 1.E-3 ) THEN
149                   ztn    = MAX( 0., tn(ji,jj,jk) - 15. )
150                   zadap  = 0.+ 1.* ztn / ( 2.+ ztn )
151                   zadap2 = 0.e0
152
153                   zfact  = EXP( -0.21 * emoy(ji,jj,jk) )
154
155                   zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * zfact )
156                   zpislopead2(ji,jj,jk) = pislope2 * ( 1.+ zadap2 * zfact )
157
158                   zpislopen = zpislopead(ji,jj,jk) * trn(ji,jj,jk,jpnch)                 &
159                     &         / ( trn(ji,jj,jk,jpphy) * 12.                   + rtrn )   &
160                     &         / ( prmax(ji,jj,jk) * rjjss * xlimphy(ji,jj,jk) + rtrn )
161
162                   zpislope2n = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)                &
163                     &          / ( trn(ji,jj,jk,jpdia) * 12.                   + rtrn )   &
164                     &          / ( prmax(ji,jj,jk) * rjjss * xlimdia(ji,jj,jk) + rtrn )
165
166!     Computation of production function
167!     ----------------------------------
168
169                   zprbio(ji,jj,jk) = prmax(ji,jj,jk) * &
170                     &                (  1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
171                   zprdia(ji,jj,jk) = prmax(ji,jj,jk) * &
172                     &                (  1.- EXP( -zpislope2n * ediat(ji,jj,jk) )  )
173               ENDIF
174            END DO
175         END DO
176      END DO
177
178
179      DO jk = 1, jpkm1
180         DO jj = 1, jpj
181            DO ji = 1, jpi
182
183                IF( etot(ji,jj,jk) > 1.E-3 ) THEN
184!    Si/C of diatoms
185!    ------------------------
186!    Si/C increases with iron stress and silicate availability
187!    Si/C is arbitrariliy increased for very high Si concentrations
188!    to mimic the very high ratios observed in the Southern Ocean (silpot2)
189
190                  zlim1  = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 )
191                  zlim   = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk)
192
193                  zsilim = MIN( zprdia(ji,jj,jk)    / ( rtrn + prmax(ji,jj,jk) ),                 &
194                  &          trn(ji,jj,jk,jpfer) / ( concdfe(ji,jj,jk) + trn(ji,jj,jk,jpfer) ),   &
195                  &          trn(ji,jj,jk,jppo4) / ( concdnh4 + trn(ji,jj,jk,jppo4) ),            &
196                  &          zlim )
197                  zsilfac = 5.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim1 - 0.5 ) )  ) + 1.e0
198                  zsiborn = MAX( 0.e0, ( trn(ji,jj,jk,jpsil) - 15.e-6 ) )
199                  zsilfac2 = 1.+ 3.* zsiborn / ( zsiborn + xksi2 )
200                  zsilfac = MIN( 6.4,zsilfac * zsilfac2)
201                  zysopt(ji,jj,jk) = grosip * zlim1 * zsilfac
202
203              ENDIF
204            END DO
205         END DO
206      END DO
207
208!    Computation of the limitation term due to
209!    A mixed layer deeper than the euphotic depth
210!    --------------------------------------------
211
212      DO jj = 1, jpj
213         DO ji = 1, jpi
214            zmxltst = MAX( 0.e0, hmld(ji,jj) - heup(ji,jj) )
215            zmxlday = zmxltst**2 / rjjss
216            zmixnano(ji,jj) = 1.- zmxlday / ( 1.+ zmxlday )
217            zmixdiat(ji,jj) = 1.- zmxlday / ( 3.+ zmxlday )
218         END DO
219      END DO
220                                                                               
221      DO jk = 1, jpkm1
222         DO jj = 1, jpj
223            DO ji = 1, jpi
224               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
225
226!     Mixed-layer effect on production
227!     --------------------------------
228                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * zmixnano(ji,jj)
229                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj)
230               ENDIF
231            END DO
232         END DO
233      END DO
234
235!     Computation of the fractionnal day length
236!     -----------------------------------------
237!      zstrn(:,:) = strn(:,:)
238
239      WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
240      zstrn(:,:) = 24. / zstrn(:,:)
241!      DO jj = 1, jpj
242!         DO ji = 1, jpi
243!
244!      Computation of the maximum light intensity
245!      ------------------------------------------
246!            IF( zstrn(ji,jj) < 1.e0 )   zstrn(ji,jj) = 24.
247!            zstrn(ji,jj) = 24. / zstrn(ji,jj)
248!         END DO
249!      END DO
250
251!CDIR NOVERRCHK
252      DO jk = 1, jpkm1
253!CDIR NOVERRCHK
254         DO jj = 1, jpj
255!CDIR NOVERRCHK
256            DO ji = 1, jpi
257
258               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
259!     Computation of the various production terms for nanophyto.
260!     ----------------------------------------------------------
261                  zetot2 = enano(ji,jj,jk) * zstrn(ji,jj)
262                  zmax = MAX( 0.1, xlimphy(ji,jj,jk) )
263                  zpislopen = zpislopead(ji,jj,jk)          &
264                  &         * trn(ji,jj,jk,jpnch) / ( rtrn + trn(ji,jj,jk,jpphy) * 12.)         &
265                  &         / ( prmax(ji,jj,jk) * rjjss * zmax + rtrn )
266
267                  zprbiochl = prmax(ji,jj,jk) * (  1.- EXP( -zpislopen * zetot2 )  )
268
269                  zprorca(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trn(ji,jj,jk,jpphy) * rfact2
270
271                  zpronew(ji,jj,jk) = zprorca(ji,jj,jk) * xnanono3(ji,jj,jk)    &
272                  &             / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn )
273                  zprod = rjjss * zprorca(ji,jj,jk) * zprbiochl * trn(ji,jj,jk,jpphy) *zmax
274
275                  zprofen(ji,jj,jk) = (fecnm)**2 * zprod / chlcnm            &
276                  &              / (  zpislopead(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpnfe) + rtrn  )
277
278                  zprochln(ji,jj,jk) = chlcnm * 144. * zprod                  &
279                  &              / (  zpislopead(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpnch) + rtrn  )
280               ENDIF
281            END DO
282         END DO
283      END DO
284
285!CDIR NOVERRCHK
286      DO jk = 1, jpkm1
287!CDIR NOVERRCHK
288         DO jj = 1, jpj
289!CDIR NOVERRCHK
290            DO ji = 1, jpi
291               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
292!       Computation of the various production terms for diatoms
293!       -------------------------------------------------------
294                  zetot2 = ediat(ji,jj,jk) * zstrn(ji,jj)
295                  zmax = MAX( 0.1, xlimdia(ji,jj,jk) )
296                  zpislope2n = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)        &
297                  &           / ( rtrn + trn(ji,jj,jk,jpdia) * 12.)        &
298                  &           / ( prmax(ji,jj,jk) * rjjss * zmax + rtrn )
299
300                  zprdiachl = prmax(ji,jj,jk) * (  1.- EXP( -zetot2 * zpislope2n )  )
301
302                  zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trn(ji,jj,jk,jpdia) * rfact2
303
304                  zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk)     &
305                  &              / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn )
306
307                  zprod = rjjss * zprorcad(ji,jj,jk) * zprdiachl * trn(ji,jj,jk,jpdia) * zmax
308
309                  zprofed(ji,jj,jk) = (fecdm)**2 * zprod / chlcdm                   &
310                  &              / ( zpislopead2(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpdfe) + rtrn )
311
312                  zprochld(ji,jj,jk) = chlcdm * 144. * zprod       &
313                  &              / ( zpislopead2(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpdch) + rtrn )
314
315               ENDIF
316            END DO
317         END DO
318      END DO
319      !
320
321!
322!   Update the arrays TRA which contain the biological sources and sinks
323!   --------------------------------------------------------------------
324
325      DO jk = 1, jpkm1
326         DO jj = 1, jpj
327           DO ji =1 ,jpi
328              zproreg  = zprorca(ji,jj,jk) - zpronew(ji,jj,jk)
329              zproreg2 = zprorcad(ji,jj,jk) - zpronewd(ji,jj,jk)
330              tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk)
331              tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - zpronew(ji,jj,jk) - zpronewd(ji,jj,jk)
332              tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zproreg - zproreg2
333              tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zprorca(ji,jj,jk) * texcret
334              tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) + zprochln(ji,jj,jk) * texcret
335              tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) + zprofen(ji,jj,jk) * texcret
336              tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorcad(ji,jj,jk) * texcret2
337              tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) + zprochld(ji,jj,jk) * texcret2
338              tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcret2
339              tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcret2
340              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + &
341              &                     excret2 * zprorcad(ji,jj,jk) + excret * zprorca(ji,jj,jk)
342              tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) &
343              &                    + ( o2ut + o2nit ) * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) )
344              tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) &
345              &                     - texcret * zprofen(ji,jj,jk) - texcret2 * zprofed(ji,jj,jk)
346              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) &
347              &                     - texcret2 * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
348              tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk)
349              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) &
350              &                    + rno3 * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) )
351          END DO
352        END DO
353     END DO
354
355     ! Total primary production per year
356     DO jk = 1, jpkm1
357        DO jj = 1, jpj
358          DO ji = 1, jpi
359             zvol = cvol(ji,jj,jk)
360#if defined key_off_degrad
361             zvol = zvol * facvol(ji,jj,jk)
362#endif
363             tpp  = tpp + ( zprorca(ji,jj,jk) + zprorcad(ji,jj,jk) ) * zvol
364          END DO
365        END DO
366      END DO
367
368      IF( lk_mpp ) CALL mpp_sum( tpp )
369
370      IF( MOD( kt, nspyr ) == 0 ) THEN
371        WRITE(numout,*) 'Total PP :'
372        WRITE(numout,*) '-------------------- : ',tpp * 12. / 1.E12
373        WRITE(numout,*) '(GtC/an)'
374        tpp = 0.
375      ENDIF
376
377#if defined key_trc_dia3d
378      zrfact2 = 1.e3 * rfact2r
379!   Supplementary diagnostics
380!   -------------------------
381      trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:)
382      trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorca(:,:,:)  * zrfact2
383      trc3d(:,:,:,jp_pcs0_3d + 5)  = zprorcad(:,:,:) * zrfact2
384      trc3d(:,:,:,jp_pcs0_3d + 6)  = zpronew(:,:,:)  * zrfact2
385      trc3d(:,:,:,jp_pcs0_3d + 7)  = zpronewd(:,:,:) * zrfact2
386      trc3d(:,:,:,jp_pcs0_3d + 8)  = zprorcad(:,:,:) * zysopt(:,:,:) * zrfact2
387      trc3d(:,:,:,jp_pcs0_3d + 9) = zprofed(:,:,:) * zrfact2
388#if ! defined key_kriest
389      trc3d(:,:,:,jp_pcs0_3d + 10) = zprofen(:,:,:) * zrfact2
390#endif
391#endif
392
393       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
394         WRITE(charout, FMT="('prod')")
395         CALL prt_ctl_trc_info(charout)
396         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
397       ENDIF
398
399   END SUBROUTINE p4z_prod
400
401   SUBROUTINE p4z_prod_init
402
403      !!----------------------------------------------------------------------
404      !!                  ***  ROUTINE p4z_prod_init  ***
405      !!
406      !! ** Purpose :   Initialization of phytoplankton production parameters
407      !!
408      !! ** Method  :   Read the nampisprod namelist and check the parameters
409      !!      called at the first timestep (nittrc000)
410      !!
411      !! ** input   :   Namelist nampisprod
412      !!
413      !!----------------------------------------------------------------------
414
415      NAMELIST/nampisprod/ pislope, pislope2, excret, excret2, chlcnm, chlcdm,   &
416         &              fecnm, fecdm, grosip
417
418      REWIND( numnat )                     ! read numnat
419      READ  ( numnat, nampisprod )
420
421      IF(lwp) THEN                         ! control print
422         WRITE(numout,*) ' '
423         WRITE(numout,*) ' Namelist parameters for phytoplankton growth, nampisprod'
424         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
425         WRITE(numout,*) '    mean Si/C ratio                           grosip    =', grosip
426         WRITE(numout,*) '    P-I slope                                 pislope   =', pislope
427         WRITE(numout,*) '    excretion ratio of nanophytoplankton      excret    =', excret
428         WRITE(numout,*) '    excretion ratio of diatoms                excret2   =', excret2
429         WRITE(numout,*) '    P-I slope  for diatoms                    pislope2  =', pislope2
430         WRITE(numout,*) '    Minimum Chl/C in nanophytoplankton        chlcnm    =', chlcnm
431         WRITE(numout,*) '    Minimum Chl/C in diatoms                  chlcdm    =', chlcdm
432         WRITE(numout,*) '    Maximum Fe/C in nanophytoplankton         fecnm     =', fecnm
433         WRITE(numout,*) '    Minimum Fe/C in diatoms                   fecdm     =', fecdm
434      ENDIF
435
436      rpis180   = rpi / 180.
437      texcret   = 1.0 - excret
438      texcret2  = 1.0 - excret2
439
440   END SUBROUTINE p4z_prod_init
441
442
443
444#else
445   !!======================================================================
446   !!  Dummy module :                                   No PISCES bio-model
447   !!======================================================================
448CONTAINS
449   SUBROUTINE p4z_prod                    ! Empty routine
450   END SUBROUTINE p4z_prod
451#endif 
452
453   !!======================================================================
454END MODULE  p4zprod
Note: See TracBrowser for help on using the repository browser.