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

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

correction of minor bug in PISCES model ( mpi case ), see ticket:327

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 17.7 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                              !: 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
236      WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
237      zstrn(:,:) = 24. / zstrn(:,:)
238
239!CDIR NOVERRCHK
240      DO jk = 1, jpkm1
241!CDIR NOVERRCHK
242         DO jj = 1, jpj
243!CDIR NOVERRCHK
244            DO ji = 1, jpi
245
246               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
247                  !     Computation of the various production terms for nanophyto.
248                  zetot2 = enano(ji,jj,jk) * zstrn(ji,jj)
249                  zmax = MAX( 0.1, xlimphy(ji,jj,jk) )
250                  zpislopen = zpislopead(ji,jj,jk)          &
251                  &         * trn(ji,jj,jk,jpnch) / ( rtrn + trn(ji,jj,jk,jpphy) * 12.)         &
252                  &         / ( prmax(ji,jj,jk) * rjjss * zmax + rtrn )
253
254                  zprbiochl = prmax(ji,jj,jk) * (  1.- EXP( -zpislopen * zetot2 )  )
255
256                  zprorca(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trn(ji,jj,jk,jpphy) * rfact2
257
258                  zpronew(ji,jj,jk) = zprorca(ji,jj,jk) * xnanono3(ji,jj,jk)    &
259                  &             / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn )
260                  zprod = rjjss * zprorca(ji,jj,jk) * zprbiochl * trn(ji,jj,jk,jpphy) *zmax
261
262                  zprofen(ji,jj,jk) = (fecnm)**2 * zprod / chlcnm            &
263                  &              / (  zpislopead(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpnfe) + rtrn  )
264
265                  zprochln(ji,jj,jk) = chlcnm * 144. * zprod                  &
266                  &              / (  zpislopead(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpnch) + rtrn  )
267               ENDIF
268            END DO
269         END DO
270      END DO
271
272!CDIR NOVERRCHK
273      DO jk = 1, jpkm1
274!CDIR NOVERRCHK
275         DO jj = 1, jpj
276!CDIR NOVERRCHK
277            DO ji = 1, jpi
278               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
279                  !  Computation of the various production terms for diatoms
280                  zetot2 = ediat(ji,jj,jk) * zstrn(ji,jj)
281                  zmax = MAX( 0.1, xlimdia(ji,jj,jk) )
282                  zpislope2n = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)        &
283                  &           / ( rtrn + trn(ji,jj,jk,jpdia) * 12.)        &
284                  &           / ( prmax(ji,jj,jk) * rjjss * zmax + rtrn )
285
286                  zprdiachl = prmax(ji,jj,jk) * (  1.- EXP( -zetot2 * zpislope2n )  )
287
288                  zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trn(ji,jj,jk,jpdia) * rfact2
289
290                  zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk)     &
291                  &              / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn )
292
293                  zprod = rjjss * zprorcad(ji,jj,jk) * zprdiachl * trn(ji,jj,jk,jpdia) * zmax
294
295                  zprofed(ji,jj,jk) = (fecdm)**2 * zprod / chlcdm                   &
296                  &              / ( zpislopead2(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpdfe) + rtrn )
297
298                  zprochld(ji,jj,jk) = chlcdm * 144. * zprod       &
299                  &              / ( zpislopead2(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpdch) + rtrn )
300
301               ENDIF
302            END DO
303         END DO
304      END DO
305      !
306
307      !   Update the arrays TRA which contain the biological sources and sinks
308      DO jk = 1, jpkm1
309         DO jj = 1, jpj
310           DO ji =1 ,jpi
311              zproreg  = zprorca(ji,jj,jk) - zpronew(ji,jj,jk)
312              zproreg2 = zprorcad(ji,jj,jk) - zpronewd(ji,jj,jk)
313              tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk)
314              tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - zpronew(ji,jj,jk) - zpronewd(ji,jj,jk)
315              tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zproreg - zproreg2
316              tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zprorca(ji,jj,jk) * texcret
317              tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) + zprochln(ji,jj,jk) * texcret
318              tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) + zprofen(ji,jj,jk) * texcret
319              tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorcad(ji,jj,jk) * texcret2
320              tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) + zprochld(ji,jj,jk) * texcret2
321              tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcret2
322              tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcret2
323              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + &
324              &                     excret2 * zprorcad(ji,jj,jk) + excret * zprorca(ji,jj,jk)
325              tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) &
326              &                    + ( o2ut + o2nit ) * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) )
327              tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) &
328              &                     - texcret * zprofen(ji,jj,jk) - texcret2 * zprofed(ji,jj,jk)
329              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) &
330              &                     - texcret2 * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
331              tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk)
332              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) &
333              &                    + rno3 * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) )
334          END DO
335        END DO
336     END DO
337
338     ! Total primary production per year
339     DO jk = 1, jpkm1
340        DO jj = 1, jpj
341          DO ji = 1, jpi
342             zvol = cvol(ji,jj,jk)
343#if defined key_off_degrad
344             zvol = zvol * facvol(ji,jj,jk)
345#endif
346             tpp  = tpp + ( zprorca(ji,jj,jk) + zprorcad(ji,jj,jk) ) &
347                          * zvol * tmask(ji,jj,jk) * tmask_i(ji,jj)
348          END DO
349        END DO
350      END DO
351
352
353      IF( MOD( kt, nspyr ) == 0 ) THEN
354        IF( lk_mpp ) CALL mpp_sum( tpp )
355        WRITE(numout,*) 'Total PP :'
356        WRITE(numout,*) '-------------------- : ',tpp * 12. / 1.E12
357        WRITE(numout,*) '(GtC/an)'
358        tpp = 0.
359      ENDIF
360
361#if defined key_trc_dia3d
362      zrfact2 = 1.e3 * rfact2r
363      !   Supplementary diagnostics
364      trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:)
365      trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorca(:,:,:)  * zrfact2
366      trc3d(:,:,:,jp_pcs0_3d + 5)  = zprorcad(:,:,:) * zrfact2
367      trc3d(:,:,:,jp_pcs0_3d + 6)  = zpronew(:,:,:)  * zrfact2
368      trc3d(:,:,:,jp_pcs0_3d + 7)  = zpronewd(:,:,:) * zrfact2
369      trc3d(:,:,:,jp_pcs0_3d + 8)  = zprorcad(:,:,:) * zysopt(:,:,:) * zrfact2
370      trc3d(:,:,:,jp_pcs0_3d + 9) = zprofed(:,:,:) * zrfact2
371#if ! defined key_kriest
372      trc3d(:,:,:,jp_pcs0_3d + 10) = zprofen(:,:,:) * zrfact2
373#endif
374#endif
375
376       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
377         WRITE(charout, FMT="('prod')")
378         CALL prt_ctl_trc_info(charout)
379         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
380       ENDIF
381
382   END SUBROUTINE p4z_prod
383
384   SUBROUTINE p4z_prod_init
385
386      !!----------------------------------------------------------------------
387      !!                  ***  ROUTINE p4z_prod_init  ***
388      !!
389      !! ** Purpose :   Initialization of phytoplankton production parameters
390      !!
391      !! ** Method  :   Read the nampisprod namelist and check the parameters
392      !!      called at the first timestep (nittrc000)
393      !!
394      !! ** input   :   Namelist nampisprod
395      !!
396      !!----------------------------------------------------------------------
397
398      NAMELIST/nampisprod/ pislope, pislope2, excret, excret2, chlcnm, chlcdm,   &
399         &              fecnm, fecdm, grosip
400
401      REWIND( numnat )                     ! read numnat
402      READ  ( numnat, nampisprod )
403
404      IF(lwp) THEN                         ! control print
405         WRITE(numout,*) ' '
406         WRITE(numout,*) ' Namelist parameters for phytoplankton growth, nampisprod'
407         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
408         WRITE(numout,*) '    mean Si/C ratio                           grosip    =', grosip
409         WRITE(numout,*) '    P-I slope                                 pislope   =', pislope
410         WRITE(numout,*) '    excretion ratio of nanophytoplankton      excret    =', excret
411         WRITE(numout,*) '    excretion ratio of diatoms                excret2   =', excret2
412         WRITE(numout,*) '    P-I slope  for diatoms                    pislope2  =', pislope2
413         WRITE(numout,*) '    Minimum Chl/C in nanophytoplankton        chlcnm    =', chlcnm
414         WRITE(numout,*) '    Minimum Chl/C in diatoms                  chlcdm    =', chlcdm
415         WRITE(numout,*) '    Maximum Fe/C in nanophytoplankton         fecnm     =', fecnm
416         WRITE(numout,*) '    Minimum Fe/C in diatoms                   fecdm     =', fecdm
417      ENDIF
418
419      rpis180   = rpi / 180.
420      texcret   = 1.0 - excret
421      texcret2  = 1.0 - excret2
422      tpp       = 0.
423
424   END SUBROUTINE p4z_prod_init
425
426
427
428#else
429   !!======================================================================
430   !!  Dummy module :                                   No PISCES bio-model
431   !!======================================================================
432CONTAINS
433   SUBROUTINE p4z_prod                    ! Empty routine
434   END SUBROUTINE p4z_prod
435#endif 
436
437   !!======================================================================
438END MODULE  p4zprod
Note: See TracBrowser for help on using the repository browser.