Index: branches/2011/dev_r2787_PISCES_improvment/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90
===================================================================
--- branches/2011/dev_r2787_PISCES_improvment/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 2822)
+++ branches/2011/dev_r2787_PISCES_improvment/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 2823)
@@ -6,4 +6,5 @@
    !! History :   1.0  !  2004     (O. Aumont) Original code
    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!             3.4  !  2011-05  (O. Aumont, C. Ethe) New parameterization of light limitation
    !!----------------------------------------------------------------------
 #if defined key_pisces
@@ -11,14 +12,16 @@
    !!   'key_pisces'                                       PISCES bio-model
    !!----------------------------------------------------------------------
-   !!   p4z_prod       :  
-   !!----------------------------------------------------------------------
-   USE trc
-   USE oce_trc         !
-   USE sms_pisces      ! 
-   USE prtctl_trc
-   USE p4zopt
-   USE p4zint
-   USE p4zlim
-   USE iom
+   !!----------------------------------------------------------------------
+   !!   p4z_prod       :   Compute the growth Rate of the two phytoplanktons groups
+   !!   p4z_prod_init  :   Initialization of the parameters for growth
+   !!   p4z_prod_alloc :   Allocate variables for growth
+   !!----------------------------------------------------------------------
+   USE oce_trc         !  shared variables between ocean and passive tracers
+   USE trc             !  passive tracers common variables 
+   USE sms_pisces      !  PISCES Source Minus Sink variables
+   USE p4zopt          !  optical model
+   USE p4zlim          !  Co-limitations of differents nutrients
+   USE prtctl_trc      !  print control for debugging
+   USE iom             !  I/O manager
 
    IMPLICIT NONE
@@ -29,22 +32,27 @@
    PUBLIC   p4z_prod_alloc
 
-   REAL(wp), PUBLIC ::   &
-     pislope   = 3.0_wp          ,  &  !:
-     pislope2  = 3.0_wp          ,  &  !:
-     excret    = 10.e-5_wp       , &   !:
-     excret2   = 0.05_wp         , &   !:
-     chlcnm    = 0.033_wp        , &   !:
-     chlcdm    = 0.05_wp         , &   !:
-     fecnm     = 10.E-6_wp       , &   !:
-     fecdm     = 15.E-6_wp       , &   !:
-     grosip    = 0.151_wp
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prmax   !:
+   !! * Shared module variables
+   REAL(wp), PUBLIC ::  pislope   = 3.0_wp            !:
+   REAL(wp), PUBLIC ::  pislope2  = 3.0_wp            !:
+   REAL(wp), PUBLIC ::  excret    = 10.e-5_wp         !:
+   REAL(wp), PUBLIC ::  excret2   = 0.05_wp           !:
+   REAL(wp), PUBLIC ::  bresp     = 0.00333_wp        !:
+   REAL(wp), PUBLIC ::  chlcnm    = 0.033_wp          !:
+   REAL(wp), PUBLIC ::  chlcdm    = 0.05_wp           !:
+   REAL(wp), PUBLIC ::  chlcmin   = 0.00333_wp        !:
+   REAL(wp), PUBLIC ::  fecnm     = 10.E-6_wp         !:
+   REAL(wp), PUBLIC ::  fecdm     = 15.E-6_wp         !:
+   REAL(wp), PUBLIC ::  grosip    = 0.151_wp          !:
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prmax     !: optimal prduction = f(temperature)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotan   !: proxy of N quota in Nanophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotad   !: proxy of N quota in diatomee
    
-   REAL(wp) ::   &
-      rday1                      ,  &  !: 0.6 / rday
-      texcret                    ,  &  !: 1 - excret 
-      texcret2                   ,  &  !: 1 - excret2        
-      tpp                              !: Total primary production
+   REAL(wp) :: r1_rday                !: 1 / rday
+   REAL(wp) :: r1_bresp               !: 1 / bresp 
+   REAL(wp) :: texcret                !: 1 - excret 
+   REAL(wp) :: texcret2               !: 1 - excret2        
+   REAL(wp) :: tpp                    !: Total primary production
+
 
    !!* Substitution
@@ -67,31 +75,33 @@
       !!---------------------------------------------------------------------
       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zmixnano   => wrk_2d_1  , zmixdiat    => wrk_2d_2  , zstrn  => wrk_2d_3
-      USE wrk_nemo, ONLY:   zpislopead => wrk_3d_2  , zpislopead2 => wrk_3d_3
-      USE wrk_nemo, ONLY:   zprdia     => wrk_3d_4  , zprbio      => wrk_3d_5  , zysopt => wrk_3d_6
-      USE wrk_nemo, ONLY:   zprorca    => wrk_3d_7  , zprorcad    => wrk_3d_8
-      USE wrk_nemo, ONLY:   zprofed    => wrk_3d_9  , zprofen     => wrk_3d_10
-      USE wrk_nemo, ONLY:   zprochln   => wrk_3d_11 , zprochld    => wrk_3d_12
-      USE wrk_nemo, ONLY:   zpronew    => wrk_3d_13 , zpronewd    => wrk_3d_14
+      USE wrk_nemo, ONLY:   zmixnano   => wrk_2d_1 , zmixdiat    => wrk_2d_2, zstrn => wrk_2d_3
+      USE wrk_nemo, ONLY:   zpislopead => wrk_3d_2 , zpislopead2 => wrk_3d_3
+      USE wrk_nemo, ONLY:   zprdia     => wrk_3d_4 , zprbio      => wrk_3d_5 
+      USE wrk_nemo, ONLY:   zprdch     => wrk_3d_6 , zprnch      => wrk_3d_7
+      USE wrk_nemo, ONLY:   zprorca    => wrk_3d_8 , zprorcad    => wrk_3d_9
+      USE wrk_nemo, ONLY:   zprofed    => wrk_3d_10, zprofen     => wrk_3d_11
+      USE wrk_nemo, ONLY:   zprochln   => wrk_3d_12, zprochld    => wrk_3d_13
+      USE wrk_nemo, ONLY:   zpronew    => wrk_3d_14, zpronewd    => wrk_3d_15
       !
       INTEGER, INTENT(in) :: kt, jnt
       !
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zsilfac, zfact
-      REAL(wp) ::   zprdiachl, zprbiochl, zsilim, ztn, zadap, zadap2
-      REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zetot2, zmax, zproreg, zproreg2
-      REAL(wp) ::   zmxltst, zmxlday, zlim1
+      REAL(wp) ::   zsilfac, zfact, znanotot, zdiattot, zconctemp, zconctemp2
+      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap
+      REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zetot2, zproreg, zproreg2
+      REAL(wp) ::   zmxltst, zmxlday
       REAL(wp) ::   zpislopen  , zpislope2n
-      REAL(wp) ::   zrum, zcodel, zargu, zval, zvol
-#if defined key_diatrc
+      REAL(wp) ::   zrum, zcodel, zargu, zval
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zysopt 
       REAL(wp) ::   zrfact2
-#endif
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
 
       IF( wrk_in_use(2, 1,2,3)                             .OR.  &
-          wrk_in_use(3, 2,3,4,5,6,7,8,9,10,11,12,13,14)  ) THEN
+          wrk_in_use(3, 2,3,4,5,6,7,8,9,10,11,12,13,14,15)  ) THEN
           CALL ctl_stop('p4z_prod: requested workspace arrays unavailable')   ;   RETURN
       ENDIF
+   
+      ALLOCATE( zysopt(jpi,jpj,jpk) )
 
       zprorca (:,:,:) = 0._wp
@@ -105,12 +115,13 @@
       zprdia  (:,:,:) = 0._wp
       zprbio  (:,:,:) = 0._wp
+      zprdch  (:,:,:) = 0._wp
+      zprnch  (:,:,:) = 0._wp
       zysopt  (:,:,:) = 0._wp
 
       ! Computation of the optimal production
-# if defined key_degrad
-      prmax(:,:,:) = rday1 * tgfunc(:,:,:) * facvol(:,:,:)
-# else
-      prmax(:,:,:) = rday1 * tgfunc(:,:,:)
-# endif
+      prmax(:,:,:) = 0.6_wp * r1_rday * tgfunc(:,:,:) 
+      IF( lk_degrad ) THEN
+        prmax(:,:,:) = prmax(:,:,:) * facvol(:,:,:)
+      ENDIF
 
       ! compute the day length depending on latitude and the day
@@ -119,15 +130,27 @@
 
       ! day length in hours
-      zstrn(:,:) = 0._wp
+      zstrn(:,:) = 0.
       DO jj = 1, jpj
          DO ji = 1, jpi
             zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
             zargu = MAX( -1., MIN(  1., zargu ) )
-            zval  = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
-            IF( zval < 1.e0 )   zval = 24.
-            zstrn(ji,jj) = 24. / zval
-         END DO
-      END DO
-
+            zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
+         END DO
+      END DO
+
+     ! Impact of the day duration on phytoplankton growth
+      DO jk = 1, jpkm1
+         DO jj = 1 ,jpj
+            DO ji = 1, jpi
+               zval = MAX( 1., zstrn(ji,jj) )
+               zval = 1.5 * zval / ( 12. + zval )
+               zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval
+               zprdia(ji,jj,jk) = zprbio(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
+      zstrn(:,:) = 24. / zstrn(:,:)
 
 !CDIR NOVERRCHK
@@ -141,26 +164,50 @@
                IF( etot(ji,jj,jk) > 1.E-3 ) THEN
                    ztn    = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
-                   zadap  = 0.+ 1.* ztn / ( 2.+ ztn )
-                   zadap2 = 0.e0
-
-                   zfact  = EXP( -0.21 * emoy(ji,jj,jk) )
-
-                   zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * zfact )
-                   zpislopead2(ji,jj,jk) = pislope2 * ( 1.+ zadap2 * zfact )
-
-                   zpislopen = zpislopead(ji,jj,jk) * trn(ji,jj,jk,jpnch)                 &
-                     &         / ( trn(ji,jj,jk,jpphy) * 12.                   + rtrn )   &
-                     &         / ( prmax(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
-
-                   zpislope2n = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)                &
-                     &          / ( trn(ji,jj,jk,jpdia) * 12.                   + rtrn )   &
-                     &          / ( prmax(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
-
-                   ! Computation of production function
-                   zprbio(ji,jj,jk) = prmax(ji,jj,jk) * &
-                     &                (  1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
-                   zprdia(ji,jj,jk) = prmax(ji,jj,jk) * &
-                     &                (  1.- EXP( -zpislope2n * ediat(ji,jj,jk) )  )
+                   zadap  = ztn / ( 2.+ ztn )
+
+                   zconctemp   = MAX( 0.e0 , trn(ji,jj,jk,jpdia) - 5e-7 )
+                   zconctemp2  = trn(ji,jj,jk,jpdia) - zconctemp
+
+                   znanotot = enano(ji,jj,jk) * zstrn(ji,jj)
+                   zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj)
+                   zfact  = EXP( -0.21 * znanotot )
+
+                   zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * zfact )  &
+                   &                       * trn(ji,jj,jk,jpnch) /( trn(ji,jj,jk,jpphy) * 12. + rtrn)
+
+                   zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)   &
+                     &                     / ( trn(ji,jj,jk,jpdia) + rtrn )   &
+                     &                     * trn(ji,jj,jk,jpdch) /( trn(ji,jj,jk,jpdia) * 12. + rtrn)
+
+                  !  Computation of production function for Chlorophyll
+                  !--------------------------------------------------
+                   zpislopen  = zpislopead (ji,jj,jk) / ( prmax(ji,jj,jk) * rday + rtrn )
+                   zpislope2n = zpislopead2(ji,jj,jk) / ( prmax(ji,jj,jk) * rday + rtrn )
+                   zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen * znanotot ) )
+                   zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) )
+
+                   ! Computation of production function for Carbon
+                   !  ---------------------------------------------
+                   zpislopen  = zpislopead (ji,jj,jk) / ( ( r1_rday + r1_bresp / chlcnm ) * rday + rtrn)
+                   zpislope2n = zpislopead2(ji,jj,jk) / ( ( r1_rday + r1_bresp / chlcdm ) * rday + rtrn)
+                   zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot )  )
+                   zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot )  )
                ENDIF
+            END DO
+         END DO
+      END DO
+
+      !  Computation of a proxy of the N/C ratio
+      !  ---------------------------------------
+!CDIR NOVERRCHK
+      DO jk = 1, jpkm1
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+                zval = ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) * prmax(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn )
+                quotan(ji,jj,jk) = MIN( 1., 0.5 + 0.5 * zval )
+                zval = ( xdiatnh4(ji,jj,jk) + xdiatno3(ji,jj,jk) ) * prmax(ji,jj,jk) / ( zprdia(ji,jj,jk) + rtrn )
+                quotad(ji,jj,jk) = MIN( 1., 0.5 + 0.5 * zval )
             END DO
          END DO
@@ -178,17 +225,11 @@
                    !    Si/C is arbitrariliy increased for very high Si concentrations
                    !    to mimic the very high ratios observed in the Southern Ocean (silpot2)
-
-                  zlim1  = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 )
-                  zlim   = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk)
-
-                  zsilim = MIN( zprdia(ji,jj,jk)    / ( rtrn + prmax(ji,jj,jk) ),                 &
-                  &          trn(ji,jj,jk,jpfer) / ( concdfe(ji,jj,jk) + trn(ji,jj,jk,jpfer) ),   &
-                  &          trn(ji,jj,jk,jppo4) / ( concdnh4 + trn(ji,jj,jk,jppo4) ),            &
-                  &          zlim )
-                  zsilfac = 5.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim1 - 0.5 ) )  ) + 1.e0
+                  zlim  = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 )
+                  zsilim = MIN( zprdia(ji,jj,jk) / ( prmax(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) )
+                  zsilfac = 4.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim - 0.5 ) )  ) + 1.e0
                   zsiborn = MAX( 0.e0, ( trn(ji,jj,jk,jpsil) - 15.e-6 ) )
-                  zsilfac2 = 1.+ 3.* zsiborn / ( zsiborn + xksi2 )
-                  zsilfac = MIN( 6.4,zsilfac * zsilfac2)
-                  zysopt(ji,jj,jk) = grosip * zlim1 * zsilfac
+                  zsilfac2 = 1.+ 2.* zsiborn / ( zsiborn + xksi2 )
+                  zsilfac = MIN( 5.4, zsilfac * zsilfac2)
+                  zysopt(ji,jj,jk) = grosip * zlim * zsilfac
               ENDIF
             END DO
@@ -201,7 +242,7 @@
          DO ji = 1, jpi
             zmxltst = MAX( 0.e0, hmld(ji,jj) - heup(ji,jj) )
-            zmxlday = zmxltst**2 / rday
-            zmixnano(ji,jj) = 1.- zmxlday / ( 1.+ zmxlday )
-            zmixdiat(ji,jj) = 1.- zmxlday / ( 3.+ zmxlday )
+            zmxlday = zmxltst * zmxltst * r1_rday
+            zmixnano(ji,jj) = 1. - zmxlday / ( 3. + zmxlday )
+            zmixdiat(ji,jj) = 1. - zmxlday / ( 4. + zmxlday )
          END DO
       END DO
@@ -214,9 +255,10 @@
                   zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * zmixnano(ji,jj)
                   zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj)
+                  zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj)
+                  zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj)
                ENDIF
             END DO
          END DO
       END DO
-
 
 !CDIR NOVERRCHK
@@ -230,22 +272,18 @@
                   !     Computation of the various production terms for nanophyto.
                   zetot2 = enano(ji,jj,jk) * zstrn(ji,jj)
-                  zmax = MAX( 0.1, xlimphy(ji,jj,jk) )
-                  zpislopen = zpislopead(ji,jj,jk)          &
-                  &         * trn(ji,jj,jk,jpnch) / ( rtrn + trn(ji,jj,jk,jpphy) * 12.)         &
-                  &         / ( prmax(ji,jj,jk) * rday * zmax + rtrn )
-
-                  zprbiochl = prmax(ji,jj,jk) * (  1.- EXP( -zpislopen * zetot2 )  )
-
                   zprorca(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trn(ji,jj,jk,jpphy) * rfact2
-
-                  zpronew(ji,jj,jk) = zprorca(ji,jj,jk) * xnanono3(ji,jj,jk)    &
-                  &             / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn )
-                  zprod = rday * zprorca(ji,jj,jk) * zprbiochl * trn(ji,jj,jk,jpphy) *zmax
-
-                  zprofen(ji,jj,jk) = (fecnm)**2 * zprod / chlcnm            &
-                  &              / (  zpislopead(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpnfe) + rtrn  )
-
-                  zprochln(ji,jj,jk) = chlcnm * 144. * zprod                  &
-                  &              / (  zpislopead(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpnch) + rtrn  )
+                  zpronew(ji,jj,jk) = zprorca(ji,jj,jk) * xnanono3(ji,jj,jk) / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn )
+                  !
+                  zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk)
+                  zprochln(ji,jj,jk) = chlcmin * 12. * zprorca(ji,jj,jk)
+                  zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + chlcnm * 12. * zprod / (  zpislopead(ji,jj,jk) * zetot2 +rtrn)
+                  !
+                  zratio = trn(ji,jj,jk,jpnfe) / ( trn(ji,jj,jk,jpphy) + rtrn )
+                  zratio = zratio / fecnm 
+                  zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
+                  zprofen(ji,jj,jk) = fecnm * prmax(ji,jj,jk)  &
+                  &             * ( 4. - 4.5 * xlimnfe(ji,jj,jk) / ( xlimnfe(ji,jj,jk) + 0.5 ) )    &
+                  &             * trn(ji,jj,jk,jpfer) / ( trn(ji,jj,jk,jpfer) + concnfe(ji,jj,jk) )  &
+                  &             * zmax * trn(ji,jj,jk,jpphy) * rfact2
                ENDIF
             END DO
@@ -262,29 +300,22 @@
                   !  Computation of the various production terms for diatoms
                   zetot2 = ediat(ji,jj,jk) * zstrn(ji,jj)
-                  zmax = MAX( 0.1, xlimdia(ji,jj,jk) )
-                  zpislope2n = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)        &
-                  &           / ( rtrn + trn(ji,jj,jk,jpdia) * 12.)        &
-                  &           / ( prmax(ji,jj,jk) * rday * zmax + rtrn )
-
-                  zprdiachl = prmax(ji,jj,jk) * (  1.- EXP( -zetot2 * zpislope2n )  )
-
                   zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trn(ji,jj,jk,jpdia) * rfact2
-
-                  zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk)     &
-                  &              / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn )
-
-                  zprod = rday * zprorcad(ji,jj,jk) * zprdiachl * trn(ji,jj,jk,jpdia) * zmax
-
-                  zprofed(ji,jj,jk) = (fecdm)**2 * zprod / chlcdm                   &
-                  &              / ( zpislopead2(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpdfe) + rtrn )
-
-                  zprochld(ji,jj,jk) = chlcdm * 144. * zprod       &
-                  &              / ( zpislopead2(ji,jj,jk) * zetot2 * trn(ji,jj,jk,jpdch) + rtrn )
-
+                  zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk) / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn )
+                  !
+                  zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk)
+                  zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk)
+                  zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + chlcdm * 12. * zprod / ( zpislopead2(ji,jj,jk) * zetot2 +rtrn )
+                  !
+                  zratio = trn(ji,jj,jk,jpdfe) / ( trn(ji,jj,jk,jpdia) + rtrn )
+                  zratio = zratio / fecdm 
+                  zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
+                  zprofed(ji,jj,jk) = fecdm * prmax(ji,jj,jk)  &
+                  &             * ( 4. - 4.5 * xlimdfe(ji,jj,jk) / ( xlimdfe(ji,jj,jk) + 0.5 ) )    &
+                  &             * trn(ji,jj,jk,jpfer) / ( trn(ji,jj,jk,jpfer) + concdfe(ji,jj,jk) )  &
+                  &             * zmax * trn(ji,jj,jk,jpdia) * rfact2
                ENDIF
             END DO
          END DO
       END DO
-      !
 
       !   Update the arrays TRA which contain the biological sources and sinks
@@ -304,15 +335,12 @@
               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcret2
               tra(ji,jj,jk,jpbsi) = tra(ji,jj,jk,jpbsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcret2
-              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + &
-              &                     excret2 * zprorcad(ji,jj,jk) + excret * zprorca(ji,jj,jk)
+              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + excret2 * zprorcad(ji,jj,jk) + excret * zprorca(ji,jj,jk)
               tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) &
-              &                    + ( o2ut + o2nit ) * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) )
-              tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) &
-              &                     - texcret * zprofen(ji,jj,jk) - texcret2 * zprofed(ji,jj,jk)
-              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) &
-              &                     - texcret2 * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
+                 &                + ( o2ut + o2nit ) * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) )
+              tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - texcret * zprofen(ji,jj,jk) - texcret2 * zprofed(ji,jj,jk)
+              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - texcret2 * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk)
-              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) &
-              &                    + rno3 * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) )
+              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) ) &
+                 &                                      - rno3 * ( zproreg + zproreg2 )
           END DO
         END DO
@@ -320,10 +348,9 @@
 
      ! Total primary production per year
-
-#if defined key_degrad
-     tpp = tpp + glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) * facvol(:,:,:) )
-#else
-     tpp = tpp + glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
-#endif
+     IF( lk_degrad )  THEN
+        tpp = tpp + glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) * facvol(:,:,:) )
+     ELSE
+        tpp = tpp + glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
+     ENDIF
 
      IF( kt == nitend .AND. jnt == nrdttrc .AND. lwp ) THEN
@@ -333,21 +360,8 @@
       ENDIF
 
-#if defined key_diatrc && ! defined key_iomput
-      !   Supplementary diagnostics
+#if defined key_diatrc 
       zrfact2 = 1.e3 * rfact2r
-      trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorca (:,:,:) * zrfact2 * tmask(:,:,:)
-      trc3d(:,:,:,jp_pcs0_3d + 5)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:)
-      trc3d(:,:,:,jp_pcs0_3d + 6)  = zpronew (:,:,:) * zrfact2 * tmask(:,:,:)
-      trc3d(:,:,:,jp_pcs0_3d + 7)  = zpronewd(:,:,:) * zrfact2 * tmask(:,:,:)
-      trc3d(:,:,:,jp_pcs0_3d + 8)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:)
-      trc3d(:,:,:,jp_pcs0_3d + 9)  = zprofed (:,:,:) * zrfact2 * tmask(:,:,:)
-#  if ! defined key_kriest
-      trc3d(:,:,:,jp_pcs0_3d + 10) = zprofen (:,:,:) * zrfact2 * tmask(:,:,:)
-#  endif
-#endif
-
-#if defined key_diatrc && defined key_iomput
-      zrfact2 = 1.e3 * rfact2r
-      IF ( jnt == nrdttrc ) then
+#if defined key_iomput
+      IF( jnt == nrdttrc ) THEN
          CALL iom_put( "PPPHY" , zprorca (:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by nanophyto
          CALL iom_put( "PPPHY2", zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by diatom
@@ -357,5 +371,15 @@
          CALL iom_put( "PFeD"  , zprofed (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by diatom
          CALL iom_put( "PFeN"  , zprofen (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by nanophyto
-      ENDIF
+#else
+         trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorca (:,:,:) * zrfact2 * tmask(:,:,:)
+         trc3d(:,:,:,jp_pcs0_3d + 5)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:)
+         trc3d(:,:,:,jp_pcs0_3d + 6)  = zpronew (:,:,:) * zrfact2 * tmask(:,:,:)
+         trc3d(:,:,:,jp_pcs0_3d + 7)  = zpronewd(:,:,:) * zrfact2 * tmask(:,:,:)
+         trc3d(:,:,:,jp_pcs0_3d + 8)  = zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:)
+         trc3d(:,:,:,jp_pcs0_3d + 9)  = zprofed (:,:,:) * zrfact2 * tmask(:,:,:)
+#  if ! defined key_kriest
+         trc3d(:,:,:,jp_pcs0_3d + 10) = zprofen (:,:,:) * zrfact2 * tmask(:,:,:)
+#  endif
+#endif
 #endif
 
@@ -367,6 +391,8 @@
 
       IF(  wrk_not_released(2, 1,2,3)                          .OR.  &
-           wrk_not_released(3, 2,3,4,5,6,7,8,9,10,11,12,13,14)   )   &
+           wrk_not_released(3, 2,3,4,5,6,7,8,9,10,11,12,13,14,15)   )   &
            CALL ctl_stop('p4z_prod: failed to release workspace arrays')
+      !
+      DEALLOCATE( zysopt )
       !
    END SUBROUTINE p4z_prod
@@ -384,10 +410,11 @@
       !! ** input   :   Namelist nampisprod
       !!----------------------------------------------------------------------
-      NAMELIST/nampisprod/ pislope, pislope2, excret, excret2, chlcnm, chlcdm,   &
-         &              fecnm, fecdm, grosip
+      !
+      NAMELIST/nampisprod/ pislope, pislope2, bresp, excret, excret2,  &
+         &                 chlcnm, chlcdm, chlcmin, fecnm, fecdm, grosip
       !!----------------------------------------------------------------------
 
-      REWIND( numnat )                     ! read numnat
-      READ  ( numnat, nampisprod )
+      REWIND( numnatp )                     ! read numnat
+      READ  ( numnatp, nampisprod )
 
       IF(lwp) THEN                         ! control print
@@ -399,15 +426,18 @@
          WRITE(numout,*) '    excretion ratio of nanophytoplankton      excret    =', excret
          WRITE(numout,*) '    excretion ratio of diatoms                excret2   =', excret2
+         WRITE(numout,*) '    basal respiration in phytoplankton        bresp     =', bresp
          WRITE(numout,*) '    P-I slope  for diatoms                    pislope2  =', pislope2
          WRITE(numout,*) '    Minimum Chl/C in nanophytoplankton        chlcnm    =', chlcnm
          WRITE(numout,*) '    Minimum Chl/C in diatoms                  chlcdm    =', chlcdm
+         WRITE(numout,*) '    Maximum Chl/C in phytoplankton            chlcmin   =', chlcmin
          WRITE(numout,*) '    Maximum Fe/C in nanophytoplankton         fecnm     =', fecnm
          WRITE(numout,*) '    Minimum Fe/C in diatoms                   fecdm     =', fecdm
       ENDIF
       !
-      rday1     = 0.6 / rday 
-      texcret   = 1.0 - excret
-      texcret2  = 1.0 - excret2
-      tpp       = 0.
+      r1_rday   = 1._wp / rday 
+      r1_bresp  = bresp * r1_rday
+      texcret   = 1._wp - excret
+      texcret2  = 1._wp - excret2
+      tpp       = 0._wp
       !
    END SUBROUTINE p4z_prod_init
@@ -418,5 +448,5 @@
       !!                     ***  ROUTINE p4z_prod_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( prmax(jpi,jpj,jpk), STAT=p4z_prod_alloc )
+      ALLOCATE( prmax(jpi,jpj,jpk), quotan(jpi,jpj,jpk), quotad(jpi,jpj,jpk), STAT = p4z_prod_alloc )
       !
       IF( p4z_prod_alloc /= 0 ) CALL ctl_warn('p4z_prod_alloc : failed to allocate arrays.')
