Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90	(revision 6841)
@@ -58,10 +58,4 @@
       IF( nn_timing == 1 )  CALL timing_start('p4z_agg')
 
-      ! Initialization of some global variables
-      ! ---------------------------------------
-      prodpoc(:,:,:) = 0.
-      conspoc(:,:,:) = 0.
-      prodgoc(:,:,:) = 0.
-      consgoc(:,:,:) = 0.
       !
       !  Exchange between organic matter compartments due to coagulation/disaggregation
@@ -100,5 +94,5 @@
                zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
                ! tranfer of DOC to POC due to brownian motion
-               zaggdoc3 =  ( 5095. * 0. * trb(ji,jj,jk,jppoc) + 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc)
+               zaggdoc3 =  ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc)
 
                !  Update the trends
@@ -109,6 +103,5 @@
                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3
                !
-               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zaggdoc + zaggdoc3
-               conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zagg
+               conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zagg + zaggdoc + zaggdoc3
                prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zagg + zaggdoc2
                !
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 6841)
@@ -150,5 +150,5 @@
       !!---------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   ztkel, ztkel1, zt , zt2   , zsal  , zsal2 , zbuf1 , zbuf2
+      REAL(wp) ::   ztkel, ztkel1, zt , zsal  , zsal2 , zbuf1 , zbuf2
       REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5
       REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2
@@ -452,5 +452,5 @@
       DO jk = 1, jpk
         DO jj = 1, jpj
-          DO ji = 1, jpj
+          DO ji = 1, jpi
             p_alkcb  = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn)
             p_dictot = trb(ji,jj,jk,jpdic) * 1000. / (rhop(ji,jj,jk) + rtrn)
@@ -533,6 +533,4 @@
    REAL(wp)  ::  znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4
    REAL(wp)  ::  znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil
-   REAL(wp)  ::  znumer_nh4, zdnumer_nh4, zdenom_nh4, zalk_nh4, zdalk_nh4
-   REAL(wp)  ::  znumer_h2s, zdnumer_h2s, zdenom_h2s, zalk_h2s, zdalk_h2s
    REAL(wp)  ::  znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4
    REAL(wp)  ::  znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu
@@ -556,5 +554,5 @@
    DO jk = 1, jpk
       DO jj = 1, jpj
-         DO ji = 1, jpj
+         DO ji = 1, jpi
             IF (rmask(ji,jj,jk) == 1.) THEN
                p_alktot = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn)
@@ -589,5 +587,5 @@
    DO jk = 1, jpk
       DO jj = 1, jpj
-         DO ji = 1, jpj
+         DO ji = 1, jpi
             IF (rmask(ji,jj,jk) == 1.) THEN
                zfact = rhop(ji,jj,jk) / 1000. + rtrn
@@ -789,9 +787,5 @@
       !!                     ***  ROUTINE p4z_che_alloc  ***
       !!----------------------------------------------------------------------
-#if defined key_ligand
       INTEGER ::   ierr(3)        ! Local variables
-#else
-      INTEGER ::   ierr(2)        ! Local variables
-#endif
       !!----------------------------------------------------------------------
 
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 6841)
@@ -73,5 +73,5 @@
       REAL(wp) ::   ztrc, zdust
 #if ! defined key_kriest
-      REAL(wp) ::   zdenom, zdenom2
+      REAL(wp) ::   zdenom2
 #endif
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig, precip
@@ -79,5 +79,5 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zstrn, zstrn2
       REAL(wp) ::   zzFeL1, zzFeL2, zzFe2, zzFeP, zzFe3, zzstrn2
-      REAL(wp) ::   zrum, zcodel, zargu, zval, zlight
+      REAL(wp) ::   zrum, zcodel, zargu, zlight
       REAL(wp) :: zkox, zkph1, zkph2, zph, zionic, ztligand
       REAL(wp) :: za, zb, zc, zkappa1, zkappa2, za0, za1, za2
@@ -172,5 +172,5 @@
                   zkox   = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn )
                   ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006)
-                  zkph2 = MAX( 0., 15. * zlight / ( zlight + 2. ) )
+                  zkph2 = MAX( 0., 15. * zlight / ( zlight + 2. ) ) * (1. - fr_i(ji,jj))
                   zkph1 = zkph2 / 5.
                   ! pass the dfe concentration from PISCES
@@ -325,5 +325,5 @@
                !  ----------------------------------------------------------------
                zlam1a  = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
-                   &   + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) + 5.09E3 * trb(ji,jj,jk,jppoc) )
+                   &   + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) + 0. * trb(ji,jj,jk,jppoc) )
                zaggdfea = zlam1a * zstep * zfecoll
 #if defined key_ligand
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90	(revision 6841)
@@ -89,5 +89,5 @@
                ! photochem loss of weak ligand
                zlablgw = MAX( 0., trn(ji,jj,jk, jpfer) * plig(ji,jj,jk) )
-               zlgwpr = prlgw * zstep * etot(ji,jj,jk) * trn(ji,jj,jk,jplgw)
+               zlgwpr = prlgw * zstep * etot(ji,jj,jk) * trn(ji,jj,jk,jplgw) * (1. - fr_i(ji,jj))
                tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zlgwp - zlgwr - zlgwpr
             END DO
@@ -104,5 +104,6 @@
                ! dissolution rate is maximal in the presence of light and 
                ! lower in the aphotici zone
-               zrfepa = rfep * (1- EXP(-1*etot(ji,jj,jk) / 25. ) ) ! 25 Wm-2 constant
+               ! ! 25 Wm-2 constant
+               zrfepa = rfep * (1- EXP(-1*etot(ji,jj,jk) / 25. ) ) * (1.- fr_i(ji,jj))
                zrfepa = MAX( (zrfepa / 10.0), zrfepa ) ! min of 10 days lifetime
                zfepr = rfep * zstep * trn(ji,jj,jk,jpfep)
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 6841)
@@ -187,5 +187,5 @@
                   &   / ( concnno3 * concnnh4 + concnnh4 * trb(ji,jj,jk,jpno3) + concnno3 * trb(ji,jj,jk,jpnh4) ) 
                zlim2  = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concnnh4 )
-               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  5.E-11   )
+               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  1.E-10   )
                ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) )
                ztem2  = tsn(ji,jj,jk,jp_tem) - 10.
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 6841)
@@ -65,5 +65,4 @@
       REAL(wp) ::   zdispot, zfact, zcalcon
       REAL(wp) ::   zomegaca, zexcess, zexcess0
-      REAL(wp) ::   zrfact2
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss, zhinit, zhi
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 6841)
@@ -244,5 +244,5 @@
                  &                + zgraztotf * unass2 - zfracfe
               zfracal = trb(ji,jj,jk,jpcal) / (trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + rtrn )
-              zgrazcal = zgrazffeg * (1. - part2) * zfracal
+              zgrazcal = (zgrazffeg + zgrazpoc) * (1. - part2) * zfracal
 #endif
 
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90	(revision 6841)
@@ -85,5 +85,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e-8 ), 0.e0 )
+               zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e-9 ), 0.e0 )
                zstep    = xstep
 # if defined key_degrad
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 6841)
@@ -78,6 +78,6 @@
       INTEGER  ::   irgb
       REAL(wp) ::   zchl
-      REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100
+      REAL(wp) ::   z1_dep
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100, zqsr_corr
 #if defined key_pisces_quota
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zetmp5
@@ -89,5 +89,5 @@
       !
       ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_alloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr_corr )
       CALL wrk_alloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 )
 #if defined key_pisces_quota
@@ -128,5 +128,6 @@
          zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr
          !
-         CALL p4z_opt_par( kt, qsr_mean, ze1, ze2, ze3 ) 
+         zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )
          !
          DO jk = 1, nksrp      
@@ -139,5 +140,6 @@
          END DO
          !
-         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 ) 
+         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )
          !
          DO jk = 1, nksrp      
@@ -149,5 +151,6 @@
          zqsr100(:,:) = 0.01 * qsr(:,:)
          !
-         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 ) 
+         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn )
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )
          !
          DO jk = 1, nksrp      
@@ -174,6 +177,7 @@
       ENDIF
       !                                        !* Euphotic depth and level
-      neln(:,:) = 1                            !  ------------------------
-      heup(:,:) = 300.
+      neln(:,:)    = 1                            !  ------------------------
+      heup(:,:)    = fsdepw(:,:,2)
+      heup_01(:,:) = fsdepw(:,:,2)
 
       DO jk = 2, nksrp
@@ -185,9 +189,14 @@
                  heup(ji,jj) = fsdepw(ji,jj,jk+1)      ! Euphotic layer depth
               ENDIF
+              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 )  THEN
+                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
+                 heup_01(ji,jj) = fsdepw(ji,jj,jk+1)      ! Euphotic layer depth
+              ENDIF
            END DO
         END DO
       END DO
       !
-      heup(:,:) = MIN( 300., heup(:,:) )
+      heup(:,:)    = MIN( 300., heup (:,:)   )
+      heup_01(:,:) = MIN( 300., heup_01(:,:) )
       !                                        !* mean light over the mixed layer
       zdepmoy(:,:)   = 0.e0                    !  -------------------------------
@@ -254,5 +263,5 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_dealloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr_corr )
       CALL wrk_dealloc( jpi, jpj, jpk, zpar,  ze0, ze1, ze2, ze3 )
 #if defined key_pisces_quota
@@ -349,5 +358,4 @@
       !! * local declarations
       INTEGER  :: ji,jj
-      REAL(wp) :: zcoef
       !!---------------------------------------------------------------------
       !
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90	(revision 6841)
@@ -30,9 +30,11 @@
    PUBLIC   p4z_poc         ! called in p4zbio.F90
    PUBLIC   p4z_poc_init    ! called in trcsms_pisces.F90
+   PUBLIC   alngam
+   PUBLIC   gamain
 
    !! * Shared module variables
    REAL(wp), PUBLIC ::  xremip     !: remineralisation rate of POC 
    INTEGER , PUBLIC ::  jcpoc      !: number of lability classes
-
+   REAL(wp), PUBLIC ::  rshape     !: shape factor of the gamma distribution
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::   alphan, reminp
@@ -62,6 +64,6 @@
       INTEGER  ::   ji, jj, jk, jn
       REAL(wp) ::   zremip, zremig, zdep, zorem2, zofer
-      REAL(wp) ::   ztem, zsizek, zsizek1, alphat, remint
-      REAL(wp) ::   solgoc, zpoc, zpoctot, zremif
+      REAL(wp) ::   zsizek, zsizek1, alphat, remint
+      REAL(wp) ::   solgoc, zpoc
 #if ! defined key_kriest
       REAL(wp) ::   zofer2, zofer3
@@ -100,4 +102,5 @@
       DO jn = 1, jcpoc
         alphag(:,:,:,jn) = alphan(jn)
+        alphap(:,:,:,jn) = alphan(jn)
       END DO
 
@@ -118,16 +121,23 @@
                 ! ------------------------------------------------------------
                 !
-                IF( fsdept(ji,jj,jk) >= zdep ) THEN
+                IF( fsdept(ji,jj,jk) > zdep ) THEN
                   alphat = 0.
                   remint = 0.
-                  ! 
-                  ! The first level just below the mixed layer needs a 
-                  ! specific treatment because lability is supposed constant
-                  ! everywhere within the mixed layer. This means that 
-                  ! change in lability in the bottom part of the previous cell
-                  ! should not be computed
-                  ! ----------------------------------------------------------
                   !
-                  IF ( fsdept(ji,jj,jk-1) < zdep ) THEN
+                  IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN
+                    ! 
+                    ! The first level just below the mixed layer needs a 
+                    ! specific treatment because lability is supposed constant
+                    ! everywhere within the mixed layer. This means that 
+                    ! change in lability in the bottom part of the previous cell
+                    ! should not be computed
+                    ! ----------------------------------------------------------
+                    !
+                    ! POC concentration is computed using the lagrangian 
+                    ! framework. It is only used for the lability param
+                    zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk) * rday / rfact2               &
+                    &   * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
+                    zpoc = max(0., zpoc)
+                    !
                     DO jn = 1, jcpoc
                        !
@@ -139,9 +149,4 @@
                        zsizek  = zdep / (wsbio2 + rtrn) * tgfunc(ji,jj,jk-1)
                        zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       ! POC concentration is computed using the lagrangian 
-                       ! framework. It is only used for the lability param
-                       zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk) * rday / rfact2               &
-                       &   * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
-                       zpoc = max(0., zpoc)
                        ! the concentration of each lability class is calculated
                        ! as the sum of the different sources and sinks
@@ -161,11 +166,12 @@
                     ! ---------------------------------------------------
                     !
+                    zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk-1) * rday / rfact2               &
+                    &   * fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk)   &
+                    &   * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
+                    zpoc = max(0., zpoc)
+                    !
                     DO jn = 1, jcpoc
                        zsizek  = fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
                        zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk-1) * rday / rfact2               &
-                       &   * fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk)   &
-                       &   * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
-                       zpoc = max(0., zpoc)
                        alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn)                &
                        &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1)  &
@@ -273,11 +279,11 @@
                       &                     * totthick(ji,jj) + totcons(ji,jj) + wsbio + rtrn )
                       alphat = alphat + alphap(ji,jj,jk,jn)
-                      remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
                    END DO
                    DO jn = 1, jcpoc
                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
+                      remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
                    END DO
                    ! Mean remineralization rate in the mixed layer
-                   zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint / ( alphat + rtrn) ))
+                   zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint ))
                 ENDIF
               ENDIF
@@ -308,12 +314,14 @@
                   !
                   IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN
+                    !
+                    ! Computation of the POC concentration using the 
+                    ! lagrangian algorithm
+                    zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk) * rday / rfact2               &
+                    &   * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
+                    zpoc = max(0., zpoc)
+                    ! 
                     DO jn = 1, jcpoc
                        ! the scale factor is corrected with temperature
                        zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       ! Computation of the POC concentration using the 
-                       ! lagrangian algorithm
-                       zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk) * rday / rfact2               &
-                       &   * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
-                       zpoc = max(0., zpoc)
                        ! computation of the lability spectrum applying the 
                        ! different sources and sinks
@@ -321,6 +329,6 @@
                        &   * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday / rfact2                   &
                        &   * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk)
+                       alphap(ji,jj,jk,jn) = MAX( 0., alphap(ji,jj,jk,jn) )
                        alphat = alphat + alphap(ji,jj,jk,jn)
-                       remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
                     END DO
                   ELSE
@@ -331,11 +339,12 @@
                     ! --------------------------------------------------------
                     !
+                    zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk-1) * rday / rfact2               &
+                    &   * fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk)   &
+                    &   * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
+                    zpoc = max(0., zpoc)
+                    !
                     DO jn = 1, jcpoc
                        zsizek  = fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
                        zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk-1) * rday / rfact2               &
-                       &   * fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk)   &
-                       &   * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
-                       zpoc = max(0., zpoc)
                        alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn)             &
                        &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1)  &
@@ -349,6 +358,6 @@
                        &  + zorem3(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday         &
                        &  / rfact2 * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk)
+                       alphap(ji,jj,jk,jn) = max(0., alphap(ji,jj,jk,jn) )
                        alphat = alphat + alphap(ji,jj,jk,jn)
-                       remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
                     END DO
                   ENDIF
@@ -357,7 +366,8 @@
                   DO jn = 1, jcpoc
                      alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
+                     remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
                   END DO
                    ! Mean remineralization rate in the water column
-                   zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint / ( alphat + rtrn) ))
+                  zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint ))
                 ENDIF
               ENDIF
@@ -436,5 +446,7 @@
       INTEGER :: jn
       REAL(wp) :: remindelta, reminup, remindown
-      NAMELIST/nampispoc/ xremip, jcpoc
+      INTEGER  :: ifault
+
+      NAMELIST/nampispoc/ xremip, jcpoc, rshape
       INTEGER :: ios                 ! Local integer output status for namelist read
 
@@ -454,4 +466,5 @@
          WRITE(numout,*) '    remineralisation rate of POC              xremip    =', xremip
          WRITE(numout,*) '    Number of lability classes for POC        jcpoc     =', jcpoc
+         WRITE(numout,*) '    Shape factor of the gamma distribution    rshape    =', rshape
       ENDIF
       !
@@ -463,18 +476,27 @@
       !
       IF (jcpoc > 1) THEN
+         !
          remindelta = log(4. * 1000. ) / float(jcpoc-1)
-         reminup = xremip/400. * exp(remindelta)
-         alphan(1) = 1. - exp(-reminup/xremip)
-         reminp(1) = xremip - reminup * exp(-reminup/xremip) / alphan(1)
+         reminup = 1./ 400. * exp(remindelta)
+         !
+         ! Discretization based on incomplete gamma functions
+         ! As incomplete gamma functions are not available in standard 
+         ! fortran 95, they have been coded as functions in this module (gamain)
+         ! ---------------------------------------------------------------------
+         !
+         alphan(1) = gamain(reminup, rshape, ifault)
+         reminp(1) = gamain(reminup, rshape+1.0, ifault) * xremip / alphan(1)
          DO jn = 2, jcpoc-1
-            reminup = xremip/400. * exp(float(jn) * remindelta)
-            remindown = xremip / 400. * exp(float(jn-1) * remindelta)
-            alphan(jn) = exp(-remindown /xremip) - exp(-reminup/xremip)
-            reminp(jn) = xremip + (remindown * exp(-remindown /xremip)    &
-            &            - reminup * exp(-reminup/xremip) ) / alphan(jn)
+            reminup = 1./ 400. * exp(float(jn) * remindelta)
+            remindown = 1. / 400. * exp(float(jn-1) * remindelta)
+            alphan(jn) = gamain(reminup, rshape, ifault) - gamain(remindown, rshape, ifault)
+            reminp(jn) = gamain(reminup, rshape+1.0, ifault) - gamain(remindown, rshape+1.0, ifault)
+            reminp(jn) = reminp(jn) * xremip / alphan(jn)
          END DO
-         remindown = xremip/400. * exp(float(jcpoc-1) * remindelta)
-         alphan(jcpoc) = exp(-remindown /xremip)
-         reminp(jcpoc) = xremip + remindown
+         remindown = 1. / 400. * exp(float(jcpoc-1) * remindelta)
+         alphan(jcpoc) = 1.0 - gamain(remindown, rshape, ifault)
+         reminp(jcpoc) = 1.0 - gamain(remindown, rshape+1.0, ifault)
+         reminp(jcpoc) = reminp(jcpoc) * xremip / alphan(jcpoc)
+
       ELSE
          alphan(jcpoc) = 1.
@@ -487,4 +509,387 @@
 
    END SUBROUTINE p4z_poc_init
+
+   REAL FUNCTION alngam( xvalue, ifault )
+
+!*****************************************************************************80
+!
+!! ALNGAM computes the logarithm of the gamma function.
+!
+!  Modified:
+!
+!    13 January 2008
+!
+!  Author:
+!
+!    Allan Macleod
+!    FORTRAN90 version by John Burkardt
+!
+!  Reference:
+!
+!    Allan Macleod,
+!    Algorithm AS 245,
+!    A Robust and Reliable Algorithm for the Logarithm of the Gamma Function,
+!    Applied Statistics,
+!    Volume 38, Number 2, 1989, pages 397-402.
+!
+!  Parameters:
+!
+!    Input, real ( kind = 8 ) XVALUE, the argument of the Gamma function.
+!
+!    Output, integer ( kind = 4 ) IFAULT, error flag.
+!    0, no error occurred.
+!    1, XVALUE is less than or equal to 0.
+!    2, XVALUE is too big.
+!
+!    Output, real ( kind = 8 ) ALNGAM, the logarithm of the gamma function of X.
+!
+  implicit none
+
+  real(wp), parameter :: alr2pi = 0.918938533204673E+00
+  integer:: ifault
+  real(wp), dimension ( 9 ) :: r1 = (/ &
+    -2.66685511495E+00, &
+    -24.4387534237E+00, &
+    -21.9698958928E+00, &
+     11.1667541262E+00, &
+     3.13060547623E+00, &
+     0.607771387771E+00, &
+     11.9400905721E+00, &
+     31.4690115749E+00, &
+     15.2346874070E+00 /)
+  real(wp), dimension ( 9 ) :: r2 = (/ &
+    -78.3359299449E+00, &
+    -142.046296688E+00, &
+     137.519416416E+00, &
+     78.6994924154E+00, &
+     4.16438922228E+00, &
+     47.0668766060E+00, &
+     313.399215894E+00, &
+     263.505074721E+00, &
+     43.3400022514E+00 /)
+  real(wp), dimension ( 9 ) :: r3 = (/ &
+    -2.12159572323E+05, &
+     2.30661510616E+05, &
+     2.74647644705E+04, &
+    -4.02621119975E+04, &
+    -2.29660729780E+03, &
+    -1.16328495004E+05, &
+    -1.46025937511E+05, &
+    -2.42357409629E+04, &
+    -5.70691009324E+02 /)
+  real(wp), dimension ( 5 ) :: r4 = (/ &
+     0.279195317918525E+00, &
+     0.4917317610505968E+00, &
+     0.0692910599291889E+00, &
+     3.350343815022304E+00, &
+     6.012459259764103E+00 /)
+  real (wp) :: x
+  real (wp) :: x1
+  real (wp) :: x2
+  real (wp), parameter :: xlge = 5.10E+05
+  real (wp), parameter :: xlgst = 1.0E+30
+  real (wp) :: xvalue
+  real (wp) :: y
+
+  x = xvalue
+  alngam = 0.0E+00
+!
+!  Check the input.
+!
+  if ( xlgst <= x ) then
+    ifault = 2
+    return
+  end if
+  if ( x <= 0.0E+00 ) then
+    ifault = 1
+    return
+  end if
+
+  ifault = 0
+!
+!  Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined.
+!
+  if ( x < 1.5E+00 ) then
+
+    if ( x < 0.5E+00 ) then
+      alngam = - log ( x )
+      y = x + 1.0E+00
+!
+!  Test whether X < machine epsilon.
+!
+      if ( y == 1.0E+00 ) then
+        return
+      end if
+
+    else
+
+      alngam = 0.0E+00
+      y = x
+      x = ( x - 0.5E+00 ) - 0.5E+00
+
+    end if
+
+    alngam = alngam + x * (((( &
+        r1(5)   * y &
+      + r1(4) ) * y &
+      + r1(3) ) * y &
+      + r1(2) ) * y &
+      + r1(1) ) / (((( &
+                  y &
+      + r1(9) ) * y &
+      + r1(8) ) * y &
+      + r1(7) ) * y &
+      + r1(6) )
+
+    return
+
+  end if
+!
+!  Calculation for 1.5 <= X < 4.0.
+!
+  if ( x < 4.0E+00 ) then
+
+    y = ( x - 1.0E+00 ) - 1.0E+00
+
+    alngam = y * (((( &
+        r2(5)   * x &
+      + r2(4) ) * x &
+      + r2(3) ) * x &
+      + r2(2) ) * x &
+      + r2(1) ) / (((( &
+                  x &
+      + r2(9) ) * x &
+      + r2(8) ) * x &
+      + r2(7) ) * x &
+      + r2(6) )
+!
+!  Calculation for 4.0 <= X < 12.0.
+!
+  else if ( x < 12.0E+00 ) then
+
+    alngam = (((( &
+        r3(5)   * x &
+      + r3(4) ) * x &
+      + r3(3) ) * x &
+      + r3(2) ) * x &
+      + r3(1) ) / (((( &
+                  x &
+      + r3(9) ) * x &
+      + r3(8) ) * x &
+      + r3(7) ) * x &
+      + r3(6) )
+!
+!  Calculation for 12.0 <= X.
+!
+  else
+
+    y = log ( x )
+    alngam = x * ( y - 1.0E+00 ) - 0.5E+00 * y + alr2pi
+
+    if ( x <= xlge ) then
+
+      x1 = 1.0E+00 / x
+      x2 = x1 * x1
+
+      alngam = alngam + x1 * ( ( &
+             r4(3)   * &
+        x2 + r4(2) ) * &
+        x2 + r4(1) ) / ( ( &
+        x2 + r4(5) ) * &
+        x2 + r4(4) )
+
+    end if
+
+   end if
+
+   END FUNCTION alngam
+
+   REAL FUNCTION gamain( x, p, ifault )
+
+!*****************************************************************************80
+!
+!! GAMAIN computes the incomplete gamma ratio.
+!
+!  Discussion:
+!
+!    A series expansion is used if P > X or X <= 1.  Otherwise, a
+!    continued fraction approximation is used.
+!
+!  Modified:
+!
+!    17 January 2008
+!
+!  Author:
+!
+!    G Bhattacharjee
+!    FORTRAN90 version by John Burkardt
+!
+!  Reference:
+!
+!    G Bhattacharjee,
+!    Algorithm AS 32:
+!    The Incomplete Gamma Integral,
+!    Applied Statistics,
+!    Volume 19, Number 3, 1970, pages 285-287.
+!
+!  Parameters:
+!
+!    Input, real ( kind = 8 ) X, P, the parameters of the incomplete 
+!    gamma ratio.  0 <= X, and 0 < P.
+!
+!    Output, integer ( kind = 4 ) IFAULT, error flag.
+!    0, no errors.
+!    1, P <= 0.
+!    2, X < 0.
+!    3, underflow.
+!    4, error return from the Log Gamma routine.
+!
+!    Output, real ( kind = 8 ) GAMAIN, the value of the incomplete
+!    gamma ratio.
+!
+  implicit none
+
+  real (wp) a
+  real (wp), parameter :: acu = 1.0E-08
+  real (wp) an
+  real (wp) arg
+  real (wp) b
+  real (wp) dif
+  real (wp) factor
+  real (wp) g
+  real (wp) gin
+  integer i
+  integer ifault
+  real (wp), parameter :: oflo = 1.0E+37
+  real (wp) p
+  real (wp) pn(6)
+  real (wp) rn
+  real (wp) term
+  real (wp), parameter :: uflo = 1.0E-37
+  real (wp) x
+!
+!  Check the input.
+!
+  if ( p <= 0.0E+00 ) then
+    ifault = 1
+    gamain = 0.0E+00
+    return
+  end if
+
+  if ( x < 0.0E+00 ) then
+    ifault = 2
+    gamain = 0.0E+00
+    return
+  end if
+
+  if ( x == 0.0E+00 ) then
+    ifault = 0
+    gamain = 0.0E+00
+    return
+  end if
+
+  g = alngam ( p, ifault )
+
+  if ( ifault /= 0 ) then
+    ifault = 4
+    gamain = 0.0E+00
+    return
+  end if
+
+  arg = p * log ( x ) - x - g
+
+  if ( arg < log ( uflo ) ) then
+    ifault = 3
+    gamain = 0.0E+00
+    return
+  end if
+
+  ifault = 0
+  factor = exp ( arg )
+!
+!  Calculation by series expansion.
+!
+  if ( x <= 1.0E+00 .or. x < p ) then
+
+    gin = 1.0E+00
+    term = 1.0E+00
+    rn = p
+
+    do
+
+      rn = rn + 1.0E+00
+      term = term * x / rn
+      gin = gin + term
+
+      if ( term <= acu ) then
+        exit
+      end if
+
+    end do
+
+    gamain = gin * factor / p
+    return
+
+  end if
+!
+!  Calculation by continued fraction.
+!
+  a = 1.0E+00 - p
+  b = a + x + 1.0E+00
+  term = 0.0E+00
+
+  pn(1) = 1.0E+00
+  pn(2) = x
+  pn(3) = x + 1.0E+00
+  pn(4) = x * b
+
+  gin = pn(3) / pn(4)
+
+  do
+
+    a = a + 1.0E+00
+    b = b + 2.0E+00
+    term = term + 1.0E+00
+    an = a * term
+    do i = 1, 2
+      pn(i+4) = b * pn(i+2) - an * pn(i)
+    end do
+
+    if ( pn(6) /= 0.0E+00 ) then
+
+      rn = pn(5) / pn(6)
+      dif = abs ( gin - rn )
+!
+!  Absolute error tolerance satisfied?
+!
+      if ( dif <= acu ) then
+!
+!  Relative error tolerance satisfied?
+!
+        if ( dif <= acu * rn ) then
+          gamain = 1.0E+00 - factor * gin
+          exit
+        end if
+
+      end if
+
+      gin = rn
+
+    end if
+
+    do i = 1, 4
+      pn(i) = pn(i+2)
+    end do
+    if ( oflo <= abs ( pn(5) ) ) then
+
+      do i = 1, 4
+        pn(i) = pn(i) / oflo
+      end do
+
+    end if
+
+  end do
+
+  END FUNCTION gamain
 
 #else
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 6841)
@@ -85,7 +85,8 @@
       REAL(wp) ::   zfact
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zmixnano, zmixdiat, zstrn, zw2d
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zstrn, zw2d
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt, zw3d   
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmxl_fac, zmxl_chl
       !!---------------------------------------------------------------------
       !
@@ -93,6 +94,7 @@
       !
       !  Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj,      zmixnano, zmixdiat, zstrn                                                  )
-      CALL wrk_alloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt            ) 
+      CALL wrk_alloc( jpi, jpj,      zstrn )
+      CALL wrk_alloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt ) 
+      CALL wrk_alloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
       CALL wrk_alloc( jpi, jpj, jpk, zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd )
       !
@@ -129,5 +131,5 @@
       END DO
 
-      ! Impact of the day duration on phytoplankton growth
+      ! Impact of the day duration and light intermittency on phytoplankton growth
       DO jk = 1, jpkm1
          DO jj = 1 ,jpj
@@ -135,7 +137,19 @@
                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                   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)
+!                  zval = 1.5 * zval / ( 12. + zval )
+!                  zmxl_fac(ji,jj,jk) = zval * ( 1. - fr_i(ji,jj))
+                  zmxl_fac(ji,jj,jk) = zval
+                  zmxl_chl(ji,jj,jk) = zval / 24. * ( 1. - fr_i(ji,jj))
+                  IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
+                     zval = MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
+                     zmxl_fac(ji,jj,jk) = zmxl_fac(ji,jj,jk) * zval
+                     zmxl_chl(ji,jj,jk) = zmxl_chl(ji,jj,jk) * zval
+                  ENDIF
+                  IF (mig(ji) == 60 .and. mjg(jj) == 25 .and. jk == 1) THEn
+                    write(0,*) 'plante ',zmxl_fac(ji,jj,jk),MAX( 1., zstrn(ji,jj) ),   &
+                    & -0.2 * zmxl_fac(ji,jj,jk)
+                  ENDIF 
+                  zmxl_fac(ji,jj,jk) = ( 1. - exp( -0.2 * zmxl_fac(ji,jj,jk) ) ) * ( 1. - fr_i(ji,jj) )
+                  zmxl_chl(ji,jj,jk) = zmxl_chl(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
                ENDIF
             END DO
@@ -143,81 +157,46 @@
       END DO
 
-      ! Maximum light intensity
-      WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
-      zstrn(:,:) = 24. / zstrn(:,:)
+      zprbio(:,:,:) = prmax(:,:,:) * zmxl_fac(:,:,:)
+      zprdia(:,:,:) = zprbio(:,:,:)
+
+      ! Computation of the P-I slope for nanos and diatoms
+      DO jk = 1, jpkm1
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+                  ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
+                  zadap       = xadap * ztn / ( 2.+ ztn )
+                  zconctemp   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
+                  zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp
+                  !
+                  zpislopead (ji,jj,jk) = pislope * ( 1.+ zadap  * EXP( -0.25 * enano(ji,jj,jk) ) )  &
+                  &                   * trb(ji,jj,jk,jpnch) /( trb(ji,jj,jk,jpphy) * 12. + rtrn)
+                  !
+                  zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp) / ( trb(ji,jj,jk,jpdia) + rtrn )   &
+                  &                   * trb(ji,jj,jk,jpdch) /( trb(ji,jj,jk,jpdia) * 12. + rtrn)
+               ENDIF
+            END DO
+         END DO
+      END DO
 
       IF( ln_newprod ) THEN
-!CDIR NOVERRCHK
          DO jk = 1, jpkm1
-!CDIR NOVERRCHK
             DO jj = 1, jpj
-!CDIR NOVERRCHK
                DO ji = 1, jpi
-                  ! Computation of the P-I slope for nanos and diatoms
                   IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                      ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
-                      zadap       = xadap * ztn / ( 2.+ ztn )
-                      zconctemp   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
-                      zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp
-                      znanotot    = enano(ji,jj,jk) * zstrn(ji,jj)
-                      zdiattot    = ediat(ji,jj,jk) * zstrn(ji,jj)
-                      !
-                      zpislopead (ji,jj,jk) = pislope * ( 1.+ zadap  * EXP( -znanotot ) )  &
-                         &                   * trb(ji,jj,jk,jpnch) /( trb(ji,jj,jk,jpphy) * 12. + rtrn)
-                      !
-                      zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp) / ( trb(ji,jj,jk,jpdia) + rtrn )   &
-                         &                   * trb(ji,jj,jk,jpdch) /( trb(ji,jj,jk,jpdia) * 12. + rtrn)
-
                       ! Computation of production function for Carbon
                       !  ---------------------------------------------
-                      zpislopen  = zpislopead (ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) * rday + rtrn)
-                      zpislope2n = zpislopead2(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) * 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 )  )
-
+                      zpislopen  = zpislopead (ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
+                      &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
+                      zpislope2n = zpislopead2(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
+                      &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
+                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) )  )
+                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) )  )
                       !  Computation of production function for Chlorophyll
                       !--------------------------------------------------
-                      zmaxday  = 1._wp / ( prmax(ji,jj,jk) * rday + rtrn )
-                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopead (ji,jj,jk) * zmaxday * znanotot ) )
-                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopead2(ji,jj,jk) * zmaxday * zdiattot ) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ELSE
-!CDIR NOVERRCHK
-         DO jk = 1, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
-               DO ji = 1, jpi
-
-                  ! Computation of the P-I slope for nanos and diatoms
-                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                      ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
-                      zadap       = ztn / ( 2.+ ztn )
-                      zconctemp   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
-                      zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp
-                      znanotot    = enano(ji,jj,jk) * zstrn(ji,jj)
-                      zdiattot    = ediat(ji,jj,jk) * zstrn(ji,jj)
-                      !
-                      zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * EXP( -znanotot ) )
-                      zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)  / ( trb(ji,jj,jk,jpdia) + rtrn )
-
-                      zpislopen =  zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch)                &
-                        &          / ( trb(ji,jj,jk,jpphy) * 12.                  + rtrn )   &
-                        &          / ( prmax(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
-
-                      zpislope2n = zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch)                &
-                        &          / ( trb(ji,jj,jk,jpdia) * 12.                  + rtrn )   &
-                        &          / ( prmax(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
-
-                      ! Computation of production function for Carbon
-                      !  ---------------------------------------------
-                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot ) )
-                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) )
-
-                      !  Computation of production function for Chlorophyll
-                      !--------------------------------------------------
+                      zpislopen  = zpislopead (ji,jj,jk) / ( prmax(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
+                      zpislope2n = zpislopead2(ji,jj,jk) / ( prmax(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
                       zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) ) )
                       zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) )
@@ -226,4 +205,25 @@
             END DO
          END DO
+      ELSE
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+                      ! Computation of production function for Carbon
+                      !  ---------------------------------------------
+                      zpislopen =  zpislopead(ji,jj,jk)  / ( zprbio(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
+                      zpislope2n = zpislopead2(ji,jj,jk) / ( zprdia(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
+                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) ) )
+                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) )
+                      !  Computation of production function for Chlorophyll
+                      !--------------------------------------------------
+                      zpislopen = zpislopen   * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
+                      zpislope2n = zpislope2n * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
+                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) ) )
+                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) )
+                  ENDIF
+               END DO
+            END DO
+         END DO
       ENDIF
 
@@ -231,9 +231,6 @@
       !  Computation of a proxy of the N/C ratio
       !  ---------------------------------------
-!CDIR NOVERRCHK
-      DO jk = 1, jpkm1
-!CDIR NOVERRCHK
+      DO jk = 1, jpkm1
          DO jj = 1, jpj
-!CDIR NOVERRCHK
             DO ji = 1, jpi
                 zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) )   &
@@ -269,26 +266,4 @@
                   zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2
               ENDIF
-            END DO
-         END DO
-      END DO
-
-      !  Computation of the limitation term due to a mixed layer deeper than the euphotic depth
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zmxltst = MAX( 0.e0, hmld(ji,jj) - heup(ji,jj) )
-            zmxlday = zmxltst * zmxltst * r1_rday
-            zmixnano(ji,jj) = 1. - zmxlday / ( 2. + zmxlday )
-            zmixdiat(ji,jj) = 1. - zmxlday / ( 4. + zmxlday )
-         END DO
-      END DO
- 
-      !  Mixed-layer effect on production                                                                               
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * zmixnano(ji,jj)
-                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj)
-               ENDIF
             END DO
          END DO
@@ -330,57 +305,24 @@
       END DO
 
-      IF( ln_newprod ) THEN
-!CDIR NOVERRCHK
-         DO jk = 1, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
-               DO ji = 1, jpi
-                  IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-                     zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj)
-                     zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj)
-                  ENDIF
-                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                     !  production terms for nanophyto. ( chlorophyll )
-                     znanotot = enano(ji,jj,jk) * zstrn(ji,jj)
-                     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-chlcmin) * 12. * zprod / &
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+                  !  production terms for nanophyto. ( chlorophyll )
+                  znanotot = enano(ji,jj,jk) / ( zmxl_chl(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-chlcmin) * 12. * zprod / &
                                         & (  zpislopead(ji,jj,jk) * znanotot +rtrn)
-                     !  production terms for diatomees ( chlorophyll )
-                     zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj)
-                     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-chlcmin) * 12. * zprod / &
-                                        & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ELSE
-!CDIR NOVERRCHK
-         DO jk = 1, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
-               DO ji = 1, jpi
-                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                     !  production terms for nanophyto. ( chlorophyll )
-                     znanotot = enano(ji,jj,jk)
-                     zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * trb(ji,jj,jk,jpphy) * xlimphy(ji,jj,jk)
-                     zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk)
-                     zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 144. * zprod            &
-                     &                    / ( zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch) * znanotot +rtrn )
-                     !  production terms for diatomees ( chlorophyll )
-                     zdiattot = ediat(ji,jj,jk)
-                     zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * trb(ji,jj,jk,jpdia) * xlimdia(ji,jj,jk)
-                     zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk)
-                     zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 144. * zprod             &
-                     &                    / ( zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch) * zdiattot +rtrn )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
+                  !  production terms for diatomees ( chlorophyll )
+                  zdiattot = ediat(ji,jj,jk) / ( zmxl_chl(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-chlcmin) * 12. * zprod / &
+                  & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn )
+               ENDIF
+            END DO
+         END DO
+      END DO
 
       !   Update the arrays TRA which contain the biological sources and sinks
@@ -550,6 +492,7 @@
      ENDIF
      !
-     CALL wrk_dealloc( jpi, jpj,      zmixnano, zmixdiat, zstrn                                                  )
-     CALL wrk_dealloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt            ) 
+     CALL wrk_dealloc( jpi, jpj,      zstrn )
+     CALL wrk_dealloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt ) 
+     CALL wrk_dealloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
      CALL wrk_dealloc( jpi, jpj, jpk, zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd )
      !
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 6841)
@@ -163,5 +163,6 @@
                ! below 2 umol/L. Inhibited at strong light 
                ! ----------------------------------------------------------
-               zonitr  =nitrif * zstep * trb(ji,jj,jk,jpnh4) / ( 1.+ emoy(ji,jj,jk) ) * ( 1.- nitrfac(ji,jj,jk) ) 
+               zonitr  = nitrif * zstep * trb(ji,jj,jk,jpnh4) * ( 1.- nitrfac(ji,jj,jk) )  &
+               &         / ( 1.+ emoy(ji,jj,jk) ) * ( 1. + fr_i(ji,jj) * emoy(ji,jj,jk) )
                denitnh4(ji,jj,jk) = nitrif * zstep * trb(ji,jj,jk,jpnh4) * nitrfac(ji,jj,jk) 
                ! Update of the tracers trends
@@ -229,5 +230,5 @@
                ! constant and specified in the namelist.
                ! ----------------------------------------------------------
-               zdep     = MAX( hmld(ji,jj), heup(ji,jj) ) 
+               zdep     = MAX( hmld(ji,jj), heup_01(ji,jj) ) 
                zdep     = MAX( 0., fsdept(ji,jj,jk) - zdep )
                ztem     = MAX( tsn(ji,jj,1,jp_tem), 0. )
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 6841)
@@ -393,5 +393,5 @@
                zlight =  ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) 
                nitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   &
-                 &         *  zfact * MIN( ztrfer, ztrpo4 ) * zlight
+                 &         *  zfact * MIN( ztrfer, ztrpo4 ) * zlight * (1. - fr_i(ji,jj))
                zsoufer(ji,jj,jk) = zlight * 2E-11 / (2E-11 + biron(ji,jj,jk))
             END DO
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 6841)
@@ -95,5 +95,5 @@
       INTEGER  ::   ji, jj, jk, jit
       INTEGER  ::   iiter1, iiter2
-      REAL(wp) ::   zfact, zwsmax, zmax, zstep
+      REAL(wp) ::   zfact, zwsmax, zmax
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
@@ -103,4 +103,10 @@
       IF( nn_timing == 1 )  CALL timing_start('p4z_sink')
 
+      ! Initialization of some global variables
+      ! ---------------------------------------
+      prodpoc(:,:,:) = 0.
+      conspoc(:,:,:) = 0.
+      prodgoc(:,:,:) = 0.
+      consgoc(:,:,:) = 0.
       !
       !    Sinking speeds of detritus is increased with depth as shown
@@ -110,5 +116,5 @@
          DO jj = 1, jpj
             DO ji = 1,jpi
-               zmax  = MAX( heup(ji,jj), hmld(ji,jj) )
+               zmax  = MAX( heup_01(ji,jj), hmld(ji,jj) )
                zfact = MAX( 0., fsdepw(ji,jj,jk+1) - zmax ) / wsbio2scale
                wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact
@@ -119,5 +125,5 @@
       ! limit the values of the sinking speeds to avoid numerical instabilities  
       wsbio3(:,:,:) = wsbio
-      wscal (:,:,:) = wsbio4(:,:,:)
+      wscal(:,:,:)  = wsbio4(:,:,:)
 #if defined key_ligand
       wsfep (:,:,:) = wfep
@@ -198,6 +204,6 @@
         CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 )
         CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 )
-        CALL p4z_sink2( wsbio4, sinksil , jpgsi, iiter2 )
-        CALL p4z_sink2( wscal , sinkcal , jpcal, iiter2 )
+        CALL p4z_sink2( wscal, sinksil , jpgsi, iiter2 )
+        CALL p4z_sink2( wscal, sinkcal , jpcal, iiter2 )
       END DO
 
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 6841)
@@ -455,5 +455,4 @@
       REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot
       CHARACTER(LEN=100)   ::   cltxt
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol
       INTEGER :: jk
       !!----------------------------------------------------------------------
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zagg.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zagg.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zagg.F90	(revision 6841)
@@ -60,12 +60,5 @@
       !
       IF( nn_timing == 1 )  CALL timing_start('p5z_agg')
-
-      ! Initialization of some global variables
-      ! ---------------------------------------
-      prodpoc(:,:,:) = 0.
-      conspoc(:,:,:) = 0.
-      prodgoc(:,:,:) = 0.
-      consgoc(:,:,:) = 0.
-
+      !
       !  Exchange between organic matter compartments due to coagulation/disaggregation
       !  ---------------------------------------------------
@@ -117,5 +110,5 @@
 
                ! tranfer of DOC to POC due to brownian motion
-               zaggtmp = ( 5095. * trb(ji,jj,jk,jppoc) + 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep
+               zaggtmp = ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep
                zaggdoc3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
                zaggdon3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
@@ -135,6 +128,5 @@
                tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) - zaggdop - zaggdop2 - zaggdop3
                !
-               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zaggdoc + zaggdoc3
-               conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zaggpoc
+               conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zagg + zaggdoc + zaggdoc3
                prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zaggpoc + zaggdoc2
                !
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zpoc.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zpoc.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zpoc.F90	(revision 6841)
@@ -32,4 +32,6 @@
    PUBLIC   p5z_poc         ! called in p4zbio.F90
    PUBLIC   p5z_poc_init    ! called in trcsms_pisces.F90
+   PUBLIC   alngam
+   PUBLIC   gamain
 
    !! * Shared module variables
@@ -38,4 +40,5 @@
    REAL(wp), PUBLIC ::  xremipp    !: remineralisation rate of DOP
    INTEGER , PUBLIC ::  jcpoc      !: number of lability classes
+   REAL(wp), PUBLIC ::  rshape     !: shape factor of the gamma distribution
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::   alphan, reminp
@@ -63,8 +66,8 @@
       !
       INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zremip, zremig, zdep, ztem, zstep
+      REAL(wp) ::   zremip, zremig, zdep, zstep
       REAL(wp) ::   zopon, zopop, zopoc2, zopon2, zopop2, zofer
       REAL(wp) ::   zsizek, zsizek1, alphat, remint
-      REAL(wp) ::   solgoc, zpoc, zpoctot, zremif
+      REAL(wp) ::   solgoc, zpoc
 #if ! defined key_kriest
       REAL(wp) ::   zofer2, zofer3
@@ -102,7 +105,13 @@
       DO jn = 1, jcpoc
         alphag(:,:,:,jn) = alphan(jn)
+        alphap(:,:,:,jn) = alphan(jn)
       END DO
 
 #if ! defined key_kriest
+! -----------------------------------------------------------------------
+! Lability parameterization. This is the big particles part (GOC)
+! This lability parameterization can be activated only with the standard
+! particle scheme. Does not work with Kriest parameterization.
+! -----------------------------------------------------------------------
      DO jk = 2, jpkm1
         DO jj = 1, jpj
@@ -110,14 +119,41 @@
               IF (tmask(ji,jj,jk) == 1.) THEN
                 zdep = hmld(ji,jj)
-                IF( fsdept(ji,jj,jk) >= zdep ) THEN
+                !
+                ! In the case of GOC, lability is constant in the mixed layer 
+                ! It is computed only below the mixed layer depth
+                ! ------------------------------------------------------------
+                !
+                IF( fsdept(ji,jj,jk) > zdep ) THEN
                   alphat = 0.
                   remint = 0.
-                  IF ( fsdept(ji,jj,jk-1) < zdep ) THEN
+                  !
+                  IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN
+                    ! 
+                    ! The first level just below the mixed layer needs a 
+                    ! specific treatment because lability is supposed constant
+                    ! everywhere within the mixed layer. This means that 
+                    ! change in lability in the bottom part of the previous cell
+                    ! should not be computed
+                    ! ----------------------------------------------------------
+                    !
+                    ! POC concentration is computed using the lagrangian 
+                    ! framework. It is only used for the lability param
+                    zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk) * rday / rfact2               &
+                    &   * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
+                    zpoc = max(0., zpoc)
+                    !
                     DO jn = 1, jcpoc
+                       !
+                       ! Lagrangian based algorithm. The fraction of each 
+                       ! lability class is computed starting from the previous
+                       ! level
+                       ! -----------------------------------------------------
+                       !
                        zsizek  = zdep / (wsbio2 + rtrn) * tgfunc(ji,jj,jk-1)
                        zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk) * rday / rfact2               &
-                       &   * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
-                       zpoc = max(0., zpoc)
+                       ! the concentration of each lability class is calculated
+                       ! as the sum of the different sources and sinks
+                       ! Please note that production of new GOC experiences
+                       ! degradation 
                        alphag(ji,jj,jk,jn) = alphan(jn) / (reminp(jn) * tgfunc(ji,jj,jk-1) )       &
                        &   * (1. - exp( -reminp(jn) * zsizek ) ) * exp( -reminp(jn) * zsizek1 )     &
@@ -128,11 +164,17 @@
                     END DO
                   ELSE
+                    !
+                    ! standard algorithm in the rest of the water column
+                    ! See the comments in the previous block.
+                    ! ---------------------------------------------------
+                    !
+                    zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk-1) * rday / rfact2               &
+                    &   * fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk)   &
+                    &   * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
+                    zpoc = max(0., zpoc)
+                    !
                     DO jn = 1, jcpoc
                        zsizek  = fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
                        zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk-1) * rday / rfact2               &
-                       &   * fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk)   &
-                       &   * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
-                       zpoc = max(0., zpoc)
                        alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn)                &
                        &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1)  &
@@ -146,6 +188,9 @@
                   ENDIF
                   DO jn = 1, jcpoc
+                     ! The contribution of each lability class at the current
+                     ! level is computed
                      alphag(ji,jj,jk,jn) = alphag(ji,jj,jk,jn) / ( alphat + rtrn)
                   END DO
+                  ! Computation of the mean remineralisation rate
                   zremigoc(ji,jj,jk) = MIN(xremipc, MAX(0., remint / ( alphat + rtrn) ))
                 ENDIF
@@ -195,23 +240,38 @@
       ENDIF
 
-      totprod(:,:) = 0.
-      totthick(:,:) = 0.
-      totcons(:,:) = 0.
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF (tmask(ji,jj,jk) == 1.) THEN
-                  zdep = hmld(ji,jj)
-                  IF( fsdept(ji,jj,jk) <= zdep ) THEN
-                     totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2
-                     totthick(ji,jj) = totthick(ji,jj) + fse3t(ji,jj,jk)* tgfunc(ji,jj,jk)
-                     totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2    &
+     ! ------------------------------------------------------------------
+     ! Lability parameterization for the small OM particles. This param 
+     ! is based on the same theoretical background as the big particles.
+     ! However, because of its low sinking speed, lability is not supposed
+     ! to be equal to its initial value (the value of the freshly produced
+     ! organic matter). It is however uniform in the mixed layer.
+     ! -------------------------------------------------------------------
+     !
+     totprod(:,:) = 0.
+     totthick(:,:) = 0.
+     totcons(:,:) = 0.
+     ! intregrated production and consumption of POC in the mixed layer
+     ! ----------------------------------------------------------------
+     ! 
+     DO jk = 1, jpkm1
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              IF (tmask(ji,jj,jk) == 1.) THEN
+                zdep = hmld(ji,jj)
+                IF( fsdept(ji,jj,jk) <= zdep ) THEN
+                  totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2
+                  ! The temperature effect is included here
+                  totthick(ji,jj) = totthick(ji,jj) + fse3t(ji,jj,jk)* tgfunc(ji,jj,jk)
+                  totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2    &
                      &                / ( trb(ji,jj,jk,jppoc) + rtrn )
-                  ENDIF
-               ENDIF
-            END DO
-         END DO
-      END DO
-
+                ENDIF
+              ENDIF
+           END DO
+        END DO
+     END DO
+
+     ! Computation of the lability spectrum in the mixed layer. In the mixed 
+     ! layer, this spectrum is supposed to be uniform.
+     ! ---------------------------------------------------------------------
      DO jk = 1, jpkm1
         DO jj = 1, jpj
@@ -223,13 +283,16 @@
                 IF( fsdept(ji,jj,jk) <= zdep ) THEN
                    DO jn = 1, jcpoc
+                      ! For each lability class, the system is supposed to be 
+                      ! at equilibrium: Prod - Sink - w alphap = 0.
                       alphap(ji,jj,jk,jn) = totprod(ji,jj) * alphan(jn) / ( reminp(jn)    &
                       &                     * totthick(ji,jj) + totcons(ji,jj) + wsbio + rtrn )
                       alphat = alphat + alphap(ji,jj,jk,jn)
-                      remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
                    END DO
                    DO jn = 1, jcpoc
                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
+                      remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
                    END DO
-                   zremipoc(ji,jj,jk) = MIN(xremipc, MAX(0., remint / ( alphat + rtrn) ))
+                   ! Mean remineralization rate in the mixed layer
+                   zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint ))
                 ENDIF
               ENDIF
@@ -237,5 +300,13 @@
         END DO
      END DO
-
+     !
+     ! -----------------------------------------------------------------------
+     ! The lability parameterization is used here. The code is here 
+     ! almost identical to what is done for big particles. The only difference
+     ! is that an additional source from GOC to POC is included. This means 
+     ! that since we need the lability spectrum of GOC, GOC spectrum 
+     ! should be determined before.
+     ! -----------------------------------------------------------------------
+     !
      DO jk = 2, jpkm1
         DO jj = 1, jpj
@@ -246,24 +317,43 @@
                   alphat = 0.
                   remint = 0.
+                  !
+                  ! Special treatment of the level just below the MXL
+                  ! See the comments in the GOC section
+                  ! ---------------------------------------------------
+                  !
                   IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN
+                    !
+                    ! Computation of the POC concentration using the 
+                    ! lagrangian algorithm
+                    zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk) * rday / rfact2               &
+                    &   * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
+                    zpoc = max(0., zpoc)
+                    !
                     DO jn = 1, jcpoc
+                       ! the scale factor is corrected with temperature
                        zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk) * rday / rfact2               &
-                       &   * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
-                       zpoc = max(0., zpoc)
+                       ! computation of the lability spectrum applying the 
+                       ! different sources and sinks
                        alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) + zorem3(ji,jj,jk)                 &
                        &   * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday / rfact2                   &
                        &   * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk)
+                       alphap(ji,jj,jk,jn) = MAX( 0., alphap(ji,jj,jk,jn) )
                        alphat = alphat + alphap(ji,jj,jk,jn)
-                       remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
                     END DO
                   ELSE
+                    !
+                    ! Lability parameterization for the interior of the ocean
+                    ! This is very similar to what is done in the previous 
+                    ! block
+                    ! --------------------------------------------------------
+                    !
+                    zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk-1) * rday / rfact2               &
+                    &   * fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk)   &
+                    &   * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
+                    zpoc = max(0., zpoc)
+                    !
                     DO jn = 1, jcpoc
                        zsizek  = fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
                        zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk-1) * rday / rfact2               &
-                       &   * fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk)   &
-                       &   * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
-                       zpoc = max(0., zpoc)
                        alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn)             &
                        &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1)  &
@@ -277,12 +367,16 @@
                        &  + zorem3(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday         &
                        &  / rfact2 * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk)
+                       alphap(ji,jj,jk,jn) = max(0., alphap(ji,jj,jk,jn) )
                        alphat = alphat + alphap(ji,jj,jk,jn)
-                       remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
                     END DO
                   ENDIF
+                  ! Normalization of the lability spectrum so that the 
+                  ! integral is equal to 1
                   DO jn = 1, jcpoc
                      alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
+                     remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
                   END DO
-                  zremipoc(ji,jj,jk) = MIN(xremipc, MAX(0., remint / ( alphat + rtrn) ))
+                  ! Mean remineralization rate in the water column
+                  zremipoc(ji,jj,jk) = MIN(xremip, MAX(0., remint ))
                 ENDIF
               ENDIF
@@ -364,5 +458,7 @@
       INTEGER :: jn
       REAL(wp) :: remindelta, reminup, remindown
-      NAMELIST/nampispoc/ xremipc, xremipn, xremipp, jcpoc
+      INTEGER  :: ifault
+
+      NAMELIST/nampispoc/ xremipc, xremipn, xremipp, jcpoc, rshape
       INTEGER :: ios                 ! Local integer output status for namelist read
 
@@ -384,4 +480,5 @@
          WRITE(numout,*) '    remineralisation rate of POP              xremipp   =', xremipp
          WRITE(numout,*) '    Number of lability classes for POC        jcpoc     =', jcpoc
+         WRITE(numout,*) '    Shape factor of the gamma distribution    rshape    =', rshape
       ENDIF
       !
@@ -394,21 +491,28 @@
       !
       IF (jcpoc > 1) THEN
+         !
          remindelta = log(4. * 1000. ) / float(jcpoc-1)
-         reminup = xremipc/400. * exp(remindelta)
-         alphan(1) = 1. - exp(-reminup/xremipc)
-         reminp(1) = xremipc - reminup * exp(-reminup/xremipc) / alphan(1)
+         reminup = 1./ 400. * exp(remindelta)
+         !
+         ! Discretization based on incomplete gamma functions
+         ! As incomplete gamma functions are not available in standard 
+         ! fortran 95, they have been coded as functions in this module (gamain)
+         ! ---------------------------------------------------------------------
+         !
+         alphan(1) = gamain(reminup, rshape, ifault)
+         reminp(1) = gamain(reminup, rshape+1.0, ifault) * xremip / alphan(1)
          DO jn = 2, jcpoc-1
-            reminup = xremipc/400. * exp(float(jn) * remindelta)
-            remindown = xremipc / 400. * exp(float(jn-1) * remindelta)
-            alphan(jn) = exp(-remindown /xremipc) - exp(-reminup/xremipc)
-            reminp(jn) = xremipc + (remindown * exp(-remindown /xremipc)    &
-            &            - reminup * exp(-reminup/xremipc) ) / alphan(jn)
+            reminup = 1./ 400. * exp(float(jn) * remindelta)
+            remindown = 1. / 400. * exp(float(jn-1) * remindelta)
+            alphan(jn) = gamain(reminup, rshape, ifault) - gamain(remindown, rshape, ifault)
+            reminp(jn) = gamain(reminup, rshape+1.0, ifault) - gamain(remindown, rshape+1.0, ifault)
+            reminp(jn) = reminp(jn) * xremip / alphan(jn)
          END DO
-         remindown = xremipc/400. * exp(float(jcpoc-1) * remindelta)
-         alphan(jcpoc) = exp(-remindown /xremipc)
-         reminp(jcpoc) = xremipc + remindown
+         remindown = 1. / 400. * exp(float(jcpoc-1) * remindelta)
+         alphan(jcpoc) = 1.0 - gamain(remindown, rshape, ifault)
+         reminp(jcpoc) = reminp(jcpoc) * xremip / alphan(jcpoc)
       ELSE
          alphan(jcpoc) = 1.
-         reminp(jcpoc) = xremipc
+         reminp(jcpoc) = xremip
       ENDIF
 
@@ -418,4 +522,382 @@
 
    END SUBROUTINE p5z_poc_init
+
+   REAL FUNCTION alngam( xvalue, ifault )
+
+!*****************************************************************************80
+!
+!! ALNGAM computes the logarithm of the gamma function.
+!
+!  Modified:
+!
+!    13 January 2008
+!
+!  Author:
+!
+!    Allan Macleod
+!    FORTRAN90 version by John Burkardt
+!
+!  Reference:
+!
+!    Allan Macleod,
+!    Algorithm AS 245,
+!    A Robust and Reliable Algorithm for the Logarithm of the Gamma Function,
+!    Applied Statistics,
+!    Volume 38, Number 2, 1989, pages 397-402.
+!
+!  Parameters:
+!
+!    Input, real ( kind = 8 ) XVALUE, the argument of the Gamma function.
+!
+!    Output, integer ( kind = 4 ) IFAULT, error flag.
+!    0, no error occurred.
+!    1, XVALUE is less than or equal to 0.
+!    2, XVALUE is too big.
+!
+!    Output, real ( kind = 8 ) ALNGAM, the logarithm of the gamma function of X.
+!
+  implicit none
+
+  real(wp), parameter :: alr2pi = 0.918938533204673E+00
+  integer:: ifault
+  real(wp), dimension ( 9 ) :: r1 = (/ &
+    -2.66685511495E+00, &
+    -24.4387534237E+00, &
+    -21.9698958928E+00, &
+     11.1667541262E+00, &
+     3.13060547623E+00, &
+     0.607771387771E+00, &
+     11.9400905721E+00, &
+     31.4690115749E+00, &
+     15.2346874070E+00 /)
+  real(wp), dimension ( 9 ) :: r2 = (/ &
+    -78.3359299449E+00, &
+    -142.046296688E+00, &
+     137.519416416E+00, &
+     78.6994924154E+00, &
+     4.16438922228E+00, &
+     47.0668766060E+00, &
+     313.399215894E+00, &
+     263.505074721E+00, &
+     43.3400022514E+00 /)
+  real(wp), dimension ( 9 ) :: r3 = (/ &
+    -2.12159572323E+05, &
+     2.30661510616E+05, &
+     2.74647644705E+04, &
+    -4.02621119975E+04, &
+    -2.29660729780E+03, &
+    -1.16328495004E+05, &
+    -1.46025937511E+05, &
+    -2.42357409629E+04, &
+    -5.70691009324E+02 /)
+  real(wp), dimension ( 5 ) :: r4 = (/ &
+     0.279195317918525E+00, &
+     0.4917317610505968E+00, &
+     0.0692910599291889E+00, &
+     3.350343815022304E+00, &
+     6.012459259764103E+00 /)
+  real (wp) :: x
+  real (wp) :: x1
+  real (wp) :: x2
+  real (wp), parameter :: xlge = 5.10E+05
+  real (wp), parameter :: xlgst = 1.0E+30
+  real (wp) :: xvalue
+  real (wp) :: y
+
+  x = xvalue
+  alngam = 0.0E+00
+!
+!  Check the input.
+!
+  if ( xlgst <= x ) then
+    ifault = 2
+    return
+  end if
+  if ( x <= 0.0E+00 ) then
+    ifault = 1
+    return
+  end if
+  ifault = 0
+!
+!  Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined.
+!
+  if ( x < 1.5E+00 ) then
+
+    if ( x < 0.5E+00 ) then
+      alngam = - log ( x )
+      y = x + 1.0E+00
+!
+!  Test whether X < machine epsilon.
+!
+      if ( y == 1.0E+00 ) then
+        return
+      end if
+
+    else
+
+      alngam = 0.0E+00
+      y = x
+      x = ( x - 0.5E+00 ) - 0.5E+00
+
+    end if
+
+    alngam = alngam + x * (((( &
+        r1(5)   * y &
+      + r1(4) ) * y &
+      + r1(3) ) * y &
+      + r1(2) ) * y &
+      + r1(1) ) / (((( &
+                  y &
+      + r1(9) ) * y &
+      + r1(8) ) * y &
+      + r1(7) ) * y &
+      + r1(6) )
+
+    return
+
+  end if
+!
+!  Calculation for 1.5 <= X < 4.0.
+!
+  if ( x < 4.0E+00 ) then
+
+    y = ( x - 1.0E+00 ) - 1.0E+00
+
+    alngam = y * (((( &
+        r2(5)   * x &
+      + r2(4) ) * x &
+      + r2(3) ) * x &
+      + r2(2) ) * x &
+      + r2(1) ) / (((( &
+                  x &
+      + r2(9) ) * x &
+      + r2(8) ) * x &
+      + r2(7) ) * x &
+      + r2(6) )
+!
+!  Calculation for 4.0 <= X < 12.0.
+!
+  else if ( x < 12.0E+00 ) then
+
+    alngam = (((( &
+        r3(5)   * x &
+      + r3(4) ) * x &
+      + r3(3) ) * x &
+      + r3(2) ) * x &
+      + r3(1) ) / (((( &
+                  x &
+      + r3(9) ) * x &
+      + r3(8) ) * x &
+      + r3(7) ) * x &
+      + r3(6) )
+!
+!  Calculation for 12.0 <= X.
+!
+  else
+
+    y = log ( x )
+    alngam = x * ( y - 1.0E+00 ) - 0.5E+00 * y + alr2pi
+
+    if ( x <= xlge ) then
+
+      x1 = 1.0E+00 / x
+      x2 = x1 * x1
+
+      alngam = alngam + x1 * ( ( &
+             r4(3)   * &
+        x2 + r4(2) ) * &
+        x2 + r4(1) ) / ( ( &
+        x2 + r4(5) ) * &
+        x2 + r4(4) )
+
+    end if
+  end if
+
+   END FUNCTION alngam
+
+   REAL FUNCTION gamain( x, p, ifault )
+
+!*****************************************************************************80
+!
+!! GAMAIN computes the incomplete gamma ratio.
+!
+!  Discussion:
+!
+!    A series expansion is used if P > X or X <= 1.  Otherwise, a
+!    continued fraction approximation is used.
+!
+!  Modified:
+!
+!    17 January 2008
+!
+!  Author:
+!
+!    G Bhattacharjee
+!    FORTRAN90 version by John Burkardt
+!
+!  Reference:
+!
+!    G Bhattacharjee,
+!    Algorithm AS 32:
+!    The Incomplete Gamma Integral,
+!    Applied Statistics,
+!    Volume 19, Number 3, 1970, pages 285-287.
+!
+!  Parameters:
+!
+!    Input, real ( kind = 8 ) X, P, the parameters of the incomplete 
+!    gamma ratio.  0 <= X, and 0 < P.
+!
+!    Output, integer ( kind = 4 ) IFAULT, error flag.
+!    0, no errors.
+!    1, P <= 0.
+!    2, X < 0.
+!    3, underflow.
+!    4, error return from the Log Gamma routine.
+!
+!    Output, real ( kind = 8 ) GAMAIN, the value of the incomplete
+!    gamma ratio.
+!
+  implicit none
+  real (wp) a
+  real (wp), parameter :: acu = 1.0E-08
+  real (wp) an
+  real (wp) arg
+  real (wp) b
+  real (wp) dif
+  real (wp) factor
+  real (wp) g
+  real (wp) gin
+  integer i
+  integer ifault
+  real (wp), parameter :: oflo = 1.0E+37
+  real (wp) p
+  real (wp) pn(6)
+  real (wp) rn
+  real (wp) term
+  real (wp), parameter :: uflo = 1.0E-37
+  real (wp) x
+!
+!  Check the input.
+!
+  if ( p <= 0.0E+00 ) then
+    ifault = 1
+    gamain = 0.0E+00
+    return
+  end if
+
+  if ( x < 0.0E+00 ) then
+    ifault = 2
+    gamain = 0.0E+00
+    return
+  end if
+
+  if ( x == 0.0E+00 ) then
+    ifault = 0
+    gamain = 0.0E+00
+    return
+  end if
+
+  g = alngam ( p, ifault )
+
+  if ( ifault /= 0 ) then
+    ifault = 4
+    gamain = 0.0E+00
+    return
+  end if
+
+  arg = p * log ( x ) - x - g
+  if ( arg < log ( uflo ) ) then
+    ifault = 3
+    gamain = 0.0E+00
+    return
+  end if
+
+  ifault = 0
+  factor = exp ( arg )
+!
+!  Calculation by series expansion.
+!
+  if ( x <= 1.0E+00 .or. x < p ) then
+
+    gin = 1.0E+00
+    term = 1.0E+00
+    rn = p
+
+    do
+
+      rn = rn + 1.0E+00
+      term = term * x / rn
+      gin = gin + term
+
+      if ( term <= acu ) then
+        exit
+      end if
+
+    end do
+
+    gamain = gin * factor / p
+    return
+
+  end if
+!
+!  Calculation by continued fraction.
+!
+  a = 1.0E+00 - p
+  b = a + x + 1.0E+00
+  term = 0.0E+00
+
+  pn(1) = 1.0E+00
+  pn(2) = x
+  pn(3) = x + 1.0E+00
+  pn(4) = x * b
+
+  gin = pn(3) / pn(4)
+  do
+
+    a = a + 1.0E+00
+    b = b + 2.0E+00
+    term = term + 1.0E+00
+    an = a * term
+    do i = 1, 2
+      pn(i+4) = b * pn(i+2) - an * pn(i)
+    end do
+
+    if ( pn(6) /= 0.0E+00 ) then
+
+      rn = pn(5) / pn(6)
+      dif = abs ( gin - rn )
+!
+!  Absolute error tolerance satisfied?
+!
+      if ( dif <= acu ) then
+!
+!  Relative error tolerance satisfied?
+!
+        if ( dif <= acu * rn ) then
+          gamain = 1.0E+00 - factor * gin
+          exit
+        end if
+
+      end if
+
+      gin = rn
+
+    end if
+
+    do i = 1, 4
+      pn(i) = pn(i+2)
+    end do
+    if ( oflo <= abs ( pn(5) ) ) then
+
+      do i = 1, 4
+        pn(i) = pn(i) / oflo
+      end do
+
+    end if
+
+  end do
+
+  END FUNCTION gamain
 
 #else
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zprod.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zprod.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zprod.F90	(revision 6841)
@@ -161,5 +161,5 @@
                IF( etot(ji,jj,jk) > 1.E-3 ) THEN
                   zval = MAX( 1., zstrn(ji,jj) )
-                  zval = 1.5 * zval / (12. + zval)
+                  zval = 1.5 * zval / (12. + zval) * (1. - fr_i(ji,jj))
                   zprbio(ji,jj,jk) = prmaxn(ji,jj,jk) * zval
                   zprpic(ji,jj,jk) = prmaxp(ji,jj,jk) * zval
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zrem.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zrem.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zrem.F90	(revision 6841)
@@ -187,8 +187,8 @@
                ! below 2 umol/L. Inhibited at strong light 
                ! ----------------------------------------------------------
-               zonitr   = nitrif * zstep * trb(ji,jj,jk,jpnh4) * (0.2 + 0.8 / ( 1.+ emoy(ji,jj,jk) ) ) &
-               &          * ( 1.- nitrfac(ji,jj,jk) ) 
+               zonitr  = nitrif * zstep * trb(ji,jj,jk,jpnh4) * ( 1.- nitrfac(ji,jj,jk) )  &
+               &         * ( 0.2 + 0.8 / ( 1.+ emoy(ji,jj,jk) ) * ( 1. + fr_i(ji,jj) * emoy(ji,jj,jk) ) )
                zdenitnh4 = nitrif * zstep * trb(ji,jj,jk,jpnh4) * nitrfac(ji,jj,jk) 
-
+               !
                ! Update of the tracers trends
                ! ----------------------------
@@ -255,5 +255,5 @@
                ! constant and specified in the namelist.
                ! ----------------------------------------------------------
-               zdep     = MAX( hmld(ji,jj), heup(ji,jj) ) 
+               zdep     = MAX( hmld(ji,jj), heup_01(ji,jj) ) 
                zdep     = MAX( 0., fsdept(ji,jj,jk) - zdep )
                ztem     = MAX( tsn(ji,jj,1,jp_tem), 0. )
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsed.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsed.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsed.F90	(revision 6841)
@@ -414,5 +414,5 @@
                ztrdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( 1E-6 + trb(ji,jj,jk,jpdop) ) * (1. - ztrpo4(ji,jj,jk))
                ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk)
-               zlight =  ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) ) 
+               zlight =  ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) ) * (1. - fr_i(ji,jj)) 
                nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight
                zsoufer(ji,jj,jk) = zlight * 2E-11 / (2E-11 + biron(ji,jj,jk))
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsink.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsink.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsink.F90	(revision 6841)
@@ -98,5 +98,5 @@
       INTEGER  ::   ji, jj, jk, jit
       INTEGER  ::   iiter1, iiter2
-      REAL(wp) ::   zfact, zwsmax, zmax, zstep
+      REAL(wp) ::   zfact, zwsmax, zmax
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
@@ -106,4 +106,10 @@
       IF( nn_timing == 1 )  CALL timing_start('p5z_sink')
 
+      ! Initialization of some global variables
+      ! ---------------------------------------
+      prodpoc(:,:,:) = 0.
+      conspoc(:,:,:) = 0.
+      prodgoc(:,:,:) = 0.
+      consgoc(:,:,:) = 0.
       !
       !    Sinking speeds of detritus is increased with depth as shown
@@ -113,5 +119,5 @@
          DO jj = 1, jpj
             DO ji = 1,jpi
-               zmax  = MAX( heup(ji,jj), hmld(ji,jj) )
+               zmax  = MAX( heup_01(ji,jj), hmld(ji,jj) )
                zfact = MAX( 0., fsdepw(ji,jj,jk+1) - zmax ) / wsbio2scale
                wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsms.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsms.F90	(revision 6840)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsms.F90	(revision 6841)
@@ -466,5 +466,4 @@
       REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot
       CHARACTER(LEN=100)   ::   cltxt
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol
       INTEGER :: jk
       !!
