Index: NEMO/trunk/src/TOP/PISCES/P4Z/p4zpoc.F90
===================================================================
--- NEMO/trunk/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 15090)
+++ NEMO/trunk/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 15459)
@@ -3,4 +3,5 @@
    !!                         ***  MODULE p4zpoc  ***
    !! TOP :   PISCES Compute remineralization of organic particles
+   !!         Same module for both PISCES and PISCES-QUOTA
    !!=========================================================================
    !! History :   1.0  !  2004     (O. Aumont) Original code
@@ -8,7 +9,9 @@
    !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Quota model for iron
    !!             3.6  !  2016-03  (O. Aumont) Quota model and diverse
+   !!             4.0  !  2018     (O. Aumont) Variable lability parameterization
    !!----------------------------------------------------------------------
    !!   p4z_poc       :  Compute remineralization/dissolution of organic compounds
    !!   p4z_poc_init  :  Initialisation of parameters for remineralisation
+   !!   alngam and gamain : computation of the incomplete gamma function
    !!----------------------------------------------------------------------
    USE oce_trc         !  shared variables between ocean and passive tracers
@@ -22,6 +25,6 @@
 
    PUBLIC   p4z_poc         ! called in p4zbio.F90
-   PUBLIC   p4z_poc_init    ! called in trcsms_pisces.F90
-   PUBLIC   alngam          !
+   PUBLIC   p4z_poc_init    ! called in trcini_pisces.F90
+   PUBLIC   alngam          ! 
    PUBLIC   gamain          !
 
@@ -33,6 +36,6 @@
    REAL(wp), PUBLIC ::   rshape     !: shape factor of the gamma distribution
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   alphan, reminp   !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   alphap           !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   alphan, reminp   !: variable lability of POC and initial distribution
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   alphap           !: lability distribution of small particles
 
 
@@ -52,6 +55,13 @@
       !!
       !! ** Purpose :   Compute remineralization of organic particles
+      !!                A reactivity-continuum parameterization is chosen
+      !!                to describe the lability of the organic particles
+      !!                As a consequence, the remineralisation rates of the 
+      !!                the different pools change with time as a function of 
+      !!                the lability distribution
       !!
-      !! ** Method  : - ???
+      !! ** Method  : - Computation of the remineralisation rates is performed
+      !!                according to reactivity continuum formalism described
+      !!                in Aumont et al. (2017). 
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, knt         ! ocean time step and ???
@@ -83,5 +93,5 @@
       solgoc = 0.04/ 2.56 * 1./ ( 1.-50**(-0.04) )
 
-      ! Initialisation of temprary arrys
+      ! Initialisation of temporary arrays
       IF( ln_p4z ) THEN
          zremipoc(:,:,:) = xremip
@@ -96,4 +106,6 @@
       zfolimi (:,:,:) = 0.
 
+      ! Initialisation of the lability distributions that are set to 
+      ! the distribution of newly produced organic particles
       DO jn = 1, jcpoc
         alphag(:,:,:,jn) = alphan(jn)
@@ -101,8 +113,8 @@
       END DO
 
-     ! -----------------------------------------------------------------------
      ! 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.
+     ! This lability parameterization is always active. However, if only one
+     ! lability class is specified in the namelist, this is equivalent to 
+     ! a standard parameterisation with a constant lability
      ! -----------------------------------------------------------------------
      ztremint(:,:,:) = zremigoc(:,:,:)
@@ -192,7 +204,10 @@
 
       IF( ln_p4z ) THEN
+         ! The standard PISCES part
          DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-            ! POC disaggregation by turbulence and bacterial activity. 
-            ! --------------------------------------------------------
+            ! POC degradation by bacterial activity. It is a function 
+            ! of the mean lability and of temperature. This also includes
+            ! shrinking of particles due to the bacterial activity
+            ! -----------------------------------------------------------
             zremig = zremigoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
             zorem2  = zremig * tr(ji,jj,jk,jpgoc,Kbb)
@@ -202,5 +217,5 @@
             zofer3 = zremig * solgoc * tr(ji,jj,jk,jpbfe,Kbb)
 
-            ! -------------------------------------
+            ! update of the TRA arrays
             tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zorem3(ji,jj,jk)
             tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) - zorem2 - zorem3(ji,jj,jk)
@@ -213,5 +228,7 @@
       ELSE
          DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-             ! POC disaggregation by turbulence and bacterial activity. 
+            ! POC degradation by bacterial activity. It is a function 
+            ! of the mean lability and of temperature. This also includes
+            ! shrinking of particles due to the bacterial activity
             ! --------------------------------------------------------
             zremig = zremigoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
@@ -223,5 +240,5 @@
             zofer2 = xremipn / xremipc * zremig * tr(ji,jj,jk,jpbfe,Kbb)
 
-            ! -------------------------------------
+            ! update of the TRA arrays
             tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zorem3(ji,jj,jk)
             tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) + solgoc * zopon2 
@@ -246,18 +263,16 @@
      ENDIF
 
-     ! ------------------------------------------------------------------
      ! 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.
-     ! -------------------------------------------------------------------
-     !
+     ! organic matter) in the MLD. 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_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
         zdep = hmld(ji,jj)
@@ -272,5 +287,6 @@
 
      ! Computation of the lability spectrum in the mixed layer. In the mixed 
-     ! layer, this spectrum is supposed to be uniform.
+     ! layer, this spectrum is supposed to be uniform as a result of intense
+     ! mixing.
      ! ---------------------------------------------------------------------
      ztremint(:,:,:) = zremipoc(:,:,:)
@@ -302,5 +318,4 @@
      ENDIF
 
-     ! -----------------------------------------------------------------------
      ! The lability parameterization is used here. The code is here 
      ! almost identical to what is done for big particles. The only difference
@@ -309,5 +324,4 @@
      ! should be determined before.
      ! -----------------------------------------------------------------------
-     !
      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1)
         IF (tmask(ji,jj,jk) == 1.) THEN
@@ -386,10 +400,12 @@
          DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
             IF (tmask(ji,jj,jk) == 1.) THEN
-              ! POC disaggregation by turbulence and bacterial activity. 
+              ! POC disaggregation by turbulence and bacterial activity.It is a function
+              ! of the mean lability and of temperature  
               ! --------------------------------------------------------
               zremip          = zremipoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
               zorem           = zremip * tr(ji,jj,jk,jppoc,Kbb)
               zofer           = zremip * tr(ji,jj,jk,jpsfe,Kbb)
-
+              
+              ! Update of the TRA arrays
               tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zorem
               orem(ji,jj,jk)      = orem(ji,jj,jk) + zorem
@@ -402,6 +418,7 @@
      ELSE
        DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-          ! POC disaggregation by turbulence and bacterial activity. 
-          ! --------------------------------------------------------
+          ! POC disaggregation by turbulence and bacterial activity.It is a function
+          ! of the mean lability and of temperature  
+          !--------------------------------------------------------
           zremip = zremipoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
           zopoc  = zremip * tr(ji,jj,jk,jppoc,Kbb)
@@ -410,5 +427,6 @@
           zopop  = xremipp / xremipc * zremip * tr(ji,jj,jk,jppop,Kbb)
           zofer  = xremipn / xremipc * zremip * tr(ji,jj,jk,jpsfe,Kbb)
-
+              
+          ! Update of the TRA arrays
           tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) - zopoc
           tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) - zopon
@@ -426,7 +444,7 @@
         IF( knt == nrdttrc ) THEN
           zrfact2 = 1.e3 * rfact2r
-          CALL iom_put( "REMINP" , zremipoc(:,:,:) * tmask(:,:,:) )  ! Remineralisation rate
-          CALL iom_put( "REMING" , zremigoc(:,:,:) * tmask(:,:,:) )  ! Remineralisation rate
-          CALL iom_put( "REMINF" , zfolimi(:,:,:)  * tmask(:,:,:)  * 1.e+9 * zrfact2 )  ! Remineralisation rate
+          CALL iom_put( "REMINP" , zremipoc(:,:,:) * tmask(:,:,:) )  ! Remineralisation rate of small particles
+          CALL iom_put( "REMING" , zremigoc(:,:,:) * tmask(:,:,:) )  ! Remineralisation rate of large particles
+          CALL iom_put( "REMINF" , zfolimi(:,:,:)  * tmask(:,:,:)  * 1.e+9 * zrfact2 )  ! Remineralisation of biogenic particulate iron
         ENDIF
      ENDIF
@@ -493,5 +511,5 @@
       ALLOCATE( alphan(jcpoc) , reminp(jcpoc) , alphap(jpi,jpj,jpk,jcpoc) )
       !
-      IF (jcpoc > 1) THEN
+      IF (jcpoc > 1) THEN  ! Case when more than one lability class is used
          !
          remindelta = LOG(4. * 1000. ) / REAL(jcpoc-1, wp)
@@ -517,5 +535,5 @@
          reminp(jcpoc) = reminp(jcpoc) * xremip / alphan(jcpoc)
 
-      ELSE
+      ELSE  ! Only one lability class is used
          alphan(jcpoc) = 1.
          reminp(jcpoc) = xremip
