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

source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90 @ 2567

Last change on this file since 2567 was 2528, checked in by rblod, 13 years ago

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