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 6453)
+++ 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 6453)
+++ 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 6453)
+++ 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 6453)
+++ 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 6453)
+++ 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 6453)
+++ 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 6453)
+++ 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
       !!
