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 @ 935

Last change on this file since 935 was 935, checked in by cetlod, 16 years ago

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