Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zagg.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zagg.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zagg.F90	(revision 12537)
@@ -2,5 +2,6 @@
    !!======================================================================
    !!                         ***  MODULE p4zagg  ***
-   !! TOP :  PISCES  aggregation of particles
+   !! TOP :  PISCES  aggregation of particles (DOC, POC, GOC)
+   !!        This module is the same for both PISCES and PISCES-QUOTA
    !!======================================================================
    !! History :   1.0  !  2004     (O. Aumont) Original code
@@ -35,7 +36,14 @@
       !!                     ***  ROUTINE p4z_agg  ***
       !!
-      !! ** Purpose :   Compute aggregation of particles
+      !! ** Purpose :   Compute aggregation of particle. Aggregation by 
+      !!                brownian motion, differential settling and shear
+      !!                are considered.
       !!
-      !! ** Method  : - ???
+      !! ** Method  : - Aggregation rates are computed assuming a fixed and 
+      !!                constant size spectrum in the different particulate 
+      !!                pools. The coagulation rates have been computed 
+      !!                externally using dedicated programs (O. Aumont). They 
+      !!                are hard-coded because they can't be changed 
+      !!                independently of each other. 
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, knt   !
@@ -53,6 +61,9 @@
       IF( ln_timing )   CALL timing_start('p4z_agg')
       !
-      !  Exchange between organic matter compartments due to coagulation/disaggregation
+      !  Exchange between organic matter compartments due to 
+      !  coagulation/disaggregation
       !  ---------------------------------------------------
+
+      ! PISCES part
       IF( ln_p4z ) THEN
          !
@@ -62,11 +73,12 @@
                   !
                   zfact = xstep * xdiss(ji,jj,jk)
-                  !  Part I : Coagulation dependent on turbulence
+                  ! Part I : Coagulation dependent on turbulence
+                  ! The stickiness has been assumed to be 0.1
                   zagg1 = 25.9  * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
                   zagg2 = 4452. * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
 
                   ! Part II : Differential settling
-
-                  !  Aggregation of small into large particles
+                  ! Aggregation of small into large particles
+                  ! The stickiness has been assumed to be 0.1
                   zagg3 =  47.1 * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
                   zagg4 =  3.3  * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
@@ -79,4 +91,5 @@
                   ! 2nd term is shear aggregation of DOC-POC
                   ! 3rd term is differential settling of DOC-POC
+                  ! 1/3 of DOC is supposed to experience aggregation (HMW)
                   zaggdoc  = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
                   &            + 2.4 * xstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc)
@@ -84,4 +97,5 @@
                   ! 1st term is shear aggregation
                   ! 2nd term is differential settling 
+                  ! 1/3 of DOC is supposed to experience aggregation (HMW)
                   zaggdoc2 = ( 3.53E3 * zfact + 0.1 * xstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
                   ! tranfer of DOC to POC due to brownian motion
@@ -103,4 +117,5 @@
       ELSE    ! ln_p5z
         !
+        ! PISCES-QUOTA part
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -109,4 +124,5 @@
                   zfact = xstep * xdiss(ji,jj,jk)
                   !  Part I : Coagulation dependent on turbulence
+                  ! The stickiness has been assumed to be 0.1
                   zaggtmp = 25.9  * zfact * trb(ji,jj,jk,jppoc)
                   zaggpoc1 = zaggtmp * trb(ji,jj,jk,jppoc)
@@ -115,4 +131,5 @@
 
                   ! Part II : Differential settling
+                  ! The stickiness has been assumed to be 0.1
    
                   !  Aggregation of small into large particles
@@ -131,4 +148,5 @@
                   ! 2nd term is shear aggregation of DOC-POC
                   ! 3rd term is differential settling of DOC-POC
+                  ! 1/3 of DOC is supposed to experience aggregation (HMW)
                   zaggtmp = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
                   &            + 2.4 * xstep * trb(ji,jj,jk,jppoc) )
@@ -140,4 +158,5 @@
                   ! 1st term is shear aggregation
                   ! 2nd term is differential settling 
+                  ! 1/3 of DOC is supposed to experience aggregation (HMW)
                   zaggtmp = ( 3.53E3 * zfact + 0.1 * xstep ) * trb(ji,jj,jk,jpgoc)
                   zaggdoc2 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
@@ -146,4 +165,5 @@
 
                   ! tranfer of DOC to POC due to brownian motion
+                  ! 1/3 of DOC is supposed to experience aggregation (HMW)
                   zaggtmp = ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) * xstep
                   zaggdoc3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zbio.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 12537)
@@ -2,8 +2,10 @@
    !!======================================================================
    !!                         ***  MODULE p4zbio  ***
-   !! TOP :   PISCES bio-model
+   !! TOP :   PISCES biogeochemical model
+   !!         This module is for both PISCES and PISCES-QUOTA
    !!======================================================================
    !! History :   1.0  !  2004     (O. Aumont) Original code
    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!             3.6  ! 2015 (O. Aumont) PISCES-QUOTA
    !!----------------------------------------------------------------------
    !!   p4z_bio        :   computes the interactions between the different
@@ -50,6 +52,6 @@
       !!
       !! ** Purpose :   Ecosystem model in the whole ocean: computes the
-      !!              different interactions between the different compartments
-      !!              of PISCES
+      !!                different interactions between the different compartments
+      !!                of PISCES
       !!
       !! ** Method  : - ???
@@ -62,14 +64,11 @@
       !
       IF( ln_timing )   CALL timing_start('p4z_bio')
-      !
-      !     ASSIGN THE SHEAR RATE THAT IS USED FOR AGGREGATION
-      !     OF PHYTOPLANKTON AND DETRITUS
 
+      ! ASSIGN THE SHEAR RATE THAT IS USED FOR AGGREGATION
+      ! OF PHYTOPLANKTON AND DETRITUS
       xdiss(:,:,:) = 1.
-!!gm the use of nmld should be better here?
       DO jk = 2, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
-!!gm  :  use nmln  and test on jk ...  less memory acces
                IF( gdepw_n(ji,jj,jk+1) > hmld(ji,jj) )   xdiss(ji,jj,jk) = 0.01
             END DO 
@@ -77,5 +76,5 @@
       END DO
 
-      CALL p4z_opt     ( kt, knt )     ! Optic: PAR in the water column
+      CALL p4z_opt     ( kt, knt )     ! Optics: PAR in the water column
       CALL p4z_sink    ( kt, knt )     ! vertical flux of particulate organic matter
       CALL p4z_fechem  ( kt, knt )     ! Iron chemistry/scavenging
@@ -87,14 +86,14 @@
          CALL p4z_mort ( kt      )     ! phytoplankton mortality
          !                             ! zooplankton sources/sinks routines 
-         CALL p4z_micro( kt, knt )           ! microzooplankton
-         CALL p4z_meso ( kt, knt )           ! mesozooplankton
+         CALL p4z_micro( kt, knt )     ! microzooplankton
+         CALL p4z_meso ( kt, knt )     ! mesozooplankton
       ELSE
          CALL p5z_lim  ( kt, knt )     ! co-limitations by the various nutrients
          CALL p5z_prod ( kt, knt )     ! phytoplankton growth rate over the global ocean. 
-         !                             ! (for each element : C, Si, Fe, Chl )
+         !                             ! (for each element : C, N, P, Si, Fe, Chl )
          CALL p5z_mort ( kt      )     ! phytoplankton mortality
          !                             ! zooplankton sources/sinks routines 
-         CALL p5z_micro( kt, knt )           ! microzooplankton
-         CALL p5z_meso ( kt, knt )           ! mesozooplankton
+         CALL p5z_micro( kt, knt )     ! microzooplankton
+         CALL p5z_meso ( kt, knt )     ! mesozooplankton
       ENDIF
       !
@@ -103,7 +102,9 @@
       CALL p4z_poc     ( kt, knt )     ! Remineralization of organic particles
       !
+      ! Ligand production. ln_ligand should be set .true. to activate
       IF( ln_ligand )  &
       & CALL p4z_ligand( kt, knt )
-      ! Update of the size of the organisms
+
+      ! Update of the size of the different phytoplankton groups
       IF (ln_p5z) THEN
          sized(:,:,:) = sizeda(:,:,:)
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zint.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zint.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zint.F90	(revision 12537)
@@ -1,7 +1,7 @@
 MODULE p4zint
-   !!======================================================================
+   !!=========================================================================
    !!                         ***  MODULE p4zint  ***
    !! TOP :   PISCES interpolation and computation of various accessory fields
-   !!======================================================================
+   !!=========================================================================
    !! History :   1.0  !  2004-03 (O. Aumont) Original code
    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
@@ -43,10 +43,14 @@
       ! Computation of phyto and zoo metabolic rate
       ! -------------------------------------------
+      ! Generic temperature dependence (Eppley, 1972)
       tgfunc (:,:,:) = EXP( 0.063913 * tsn(:,:,:,jp_tem) )
+      ! Temperature dependence of mesozooplankton (Buitenhuis et al. (2005))
       tgfunc2(:,:,:) = EXP( 0.07608  * tsn(:,:,:,jp_tem) )
+      ! Temperature dependence of picophytoplankton (Stawiarsky et al., 2016)
       tgfunc3(:,:,:) = EXP( 0.0825   * tsn(:,:,:,jp_tem) )
 
       ! Computation of the silicon dependant half saturation  constant for silica uptake
-      ! ---------------------------------------------------
+      ! This is based on an old study by Pondaven et al. (1998)
+      ! --------------------------------------------------------------------------------
       DO ji = 1, jpi
          DO jj = 1, jpj
@@ -56,4 +60,8 @@
       END DO
       !
+      ! At the end of each year, the half saturation constant for silica is 
+      ! updated as this is based on the highest concentration reached over 
+      ! the year
+      ! -------------------------------------------------------------------
       IF( nday_year == nyear_len(1) ) THEN
          xksi   (:,:) = xksimax(:,:)
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlim.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlim.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlim.F90	(revision 12537)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                         ***  MODULE p4zlim  ***
-   !! TOP :   PISCES 
+   !! TOP :   Computes the nutrient limitation terms of phytoplankton
    !!======================================================================
    !! History :   1.0  !  2004     (O. Aumont) Original code
@@ -14,17 +14,17 @@
    USE trc             ! Tracers defined
    USE sms_pisces      ! PISCES variables
-   USE iom             !  I/O manager
+   USE iom             ! I/O manager
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC p4z_lim    
-   PUBLIC p4z_lim_init    
-   PUBLIC p4z_lim_alloc
+   PUBLIC p4z_lim           ! called in p4zbio.F90 
+   PUBLIC p4z_lim_init      ! called in trcsms_pisces.F90 
+   PUBLIC p4z_lim_alloc     ! called in trcini_pisces.F90
 
    !! * Shared module variables
    REAL(wp), PUBLIC ::  concnno3    !:  NO3, PO4 half saturation   
    REAL(wp), PUBLIC ::  concdno3    !:  Phosphate half saturation for diatoms  
-   REAL(wp), PUBLIC ::  concnnh4    !:  NH4 half saturation for phyto  
+   REAL(wp), PUBLIC ::  concnnh4    !:  NH4 half saturation for nanophyto  
    REAL(wp), PUBLIC ::  concdnh4    !:  NH4 half saturation for diatoms
    REAL(wp), PUBLIC ::  concnfer    !:  Iron half saturation for nanophyto 
@@ -46,21 +46,21 @@
 
    !!* Phytoplankton limitation terms
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanono3   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatno3   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanonh4   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatnh4   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanopo4   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatpo4   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimphy    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdia    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimnfe    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdfe    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimsi     !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimbac    !: ??
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimbacl   !: ??
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concdfe    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concnfe    !: ???
-
-   ! Coefficient for iron limitation
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanono3   !: Nanophyto limitation by NO3
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatno3   !: Diatoms limitation by NO3
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanonh4   !: Nanophyto limitation by NH4
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatnh4   !:  Diatoms limitation by NH4
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanopo4   !: Nanophyto limitation by PO4
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatpo4   !: Diatoms limitation by PO4
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimphy    !: Nutrient limitation term of nanophytoplankton
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdia    !: Nutrient limitation term of diatoms
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimnfe    !: Nanophyto limitation by Iron
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdfe    !: Diatoms limitation by iron
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimsi     !: Diatoms limitation by Si
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimbac    !: Bacterial limitation term
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimbacl   !: Bacterial limitation term
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concdfe    !: Limitation of diatoms uptake of Fe
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concnfe    !: Limitation of Nano uptake of Fe
+
+   ! Coefficient for iron limitation following Flynn and Hipkin (1999)
    REAL(wp) ::  xcoef1   = 0.0016  / 55.85  
    REAL(wp) ::  xcoef2   = 1.21E-5 * 14. / 55.85 / 7.625 * 0.5 * 1.5
@@ -79,7 +79,7 @@
       !!
       !! ** Purpose :   Compute the co-limitations by the various nutrients
-      !!              for the various phytoplankton species
-      !!
-      !! ** Method  : - ???
+      !!                for the various phytoplankton species
+      !!
+      !! ** Method  : - Limitation follows the Liebieg law of the minimum
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in)  :: kt, knt
@@ -99,6 +99,7 @@
             DO ji = 1, jpi
                
-               ! Tuning of the iron concentration to a minimum level that is set to the detection limit
-               !-------------------------------------
+               ! Tuning of the iron concentration to a minimum level that 
+               ! is set to the detection limit
+               ! --------------------------------------------------------
                zno3    = trb(ji,jj,jk,jpno3) / 40.e-6
                zferlim = MAX( 3e-11 * zno3 * zno3, 5e-12 )
@@ -106,7 +107,7 @@
                trb(ji,jj,jk,jpfer) = MAX( trb(ji,jj,jk,jpfer), zferlim )
 
-               ! Computation of a variable Ks for iron on diatoms taking into account
+               ! Computation of a variable Ks of diatoms taking into account
                ! that increasing biomass is made of generally bigger cells
-               !------------------------------------------------
+               !------------------------------------------------------------
                zconcd   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
                zconcd2  = trb(ji,jj,jk,jpdia) - zconcd
@@ -124,6 +125,7 @@
                zconc0nnh4        = MAX( concnnh4, ( zconcn2 * concnnh4 + concnnh4 * xsizern * zconcn ) * z1_trbphy )
 
-               ! Michaelis-Menten Limitation term for nutrients Small bacteria
-               ! -------------------------------------------------------------
+               ! Michaelis-Menten Limitation term by nutrients of
+               ! heterotrophic bacteria
+               ! -------------------------------------------------
                zdenom = 1. /  ( concbno3 * concbnh4 + concbnh4 * trb(ji,jj,jk,jpno3) + concbno3 * trb(ji,jj,jk,jpnh4) )
                xnanono3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * concbnh4 * zdenom
@@ -134,9 +136,12 @@
                zlim3    = trb(ji,jj,jk,jpfer) / ( concbfe + trb(ji,jj,jk,jpfer) )
                zlim4    = trb(ji,jj,jk,jpdoc) / ( xkdoc   + trb(ji,jj,jk,jpdoc) )
+               ! Xlimbac is used for DOC solubilization whereas xlimbacl
+               ! is used for all the other bacterial-dependent terms
+               ! -------------------------------------------------------
                xlimbacl(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 )
                xlimbac (ji,jj,jk) = MIN( zlim1, zlim2, zlim3 ) * zlim4
 
-               ! Michaelis-Menten Limitation term for nutrients Small flagellates
-               ! -----------------------------------------------
+               ! Michaelis-Menten Limitation term by nutrients: Nanophyto
+               ! --------------------------------------------------------
                zdenom = 1. /  ( zconc0n * zconc0nnh4 + zconc0nnh4 * trb(ji,jj,jk,jpno3) + zconc0n * trb(ji,jj,jk,jpnh4) )
                xnanono3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * zconc0nnh4 * zdenom
@@ -146,4 +151,8 @@
                zlim2    = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zconc0nnh4 )
                zratio   = trb(ji,jj,jk,jpnfe) * z1_trbphy 
+
+               ! The minimum iron quota depends on the size of PSU, respiration
+               ! and the reduction of nitrate following the parameterization 
+               ! proposed by Flynn and Hipkin (1999)
                zironmin = xcoef1 * trb(ji,jj,jk,jpnch) * z1_trbphy + xcoef2 * zlim1 + xcoef3 * xnanono3(ji,jj,jk)
                zlim3    = MAX( 0.,( zratio - zironmin ) / qnfelim )
@@ -151,7 +160,7 @@
                xlimnfe (ji,jj,jk) = MIN( 1., zlim3 )
                xlimphy (ji,jj,jk) = MIN( zlim1, zlim2, zlim3 )
-               !
-               !   Michaelis-Menten Limitation term for nutrients Diatoms
-               !   ----------------------------------------------
+               
+               !   Michaelis-Menten Limitation term by nutrients : Diatoms
+               !   -------------------------------------------------------
                zdenom   = 1. / ( zconc1d * zconc1dnh4 + zconc1dnh4 * trb(ji,jj,jk,jpno3) + zconc1d * trb(ji,jj,jk,jpnh4) )
                xdiatno3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * zconc1dnh4 * zdenom
@@ -162,4 +171,8 @@
                zlim3    = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi(ji,jj) )
                zratio   = trb(ji,jj,jk,jpdfe) * z1_trbdia
+
+               ! The minimum iron quota depends on the size of PSU, respiration
+               ! and the reduction of nitrate following the parameterization 
+               ! proposed by Flynn and Hipkin (1999)
                zironmin = xcoef1 * trb(ji,jj,jk,jpdch) * z1_trbdia + xcoef2 * zlim1 + xcoef3 * xdiatno3(ji,jj,jk)
                zlim4    = MAX( 0., ( zratio - zironmin ) / qdfelim )
@@ -173,4 +186,6 @@
 
       ! Compute the fraction of nanophytoplankton that is made of calcifiers
+      ! This is a purely adhoc formulation described in Aumont et al. (2015)
+      ! This fraction depends on nutrient limitation, light, temperature
       ! --------------------------------------------------------------------
       DO jk = 1, jpkm1
@@ -206,5 +221,5 @@
                nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
                !
-               ! denitrification factor computed from NO3 levels
+               ! redox factor computed from NO3 levels
                nitrfac2(ji,jj,jk) = MAX( 0.e0,       ( 1.E-6 - trb(ji,jj,jk,jpno3) )  &
                   &                                / ( 1.E-6 + trb(ji,jj,jk,jpno3) ) )
@@ -215,5 +230,5 @@
       !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics
-        IF( iom_use( "xfracal" ) )   CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! euphotic layer deptht
+        IF( iom_use( "xfracal" ) )   CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! Faction of nanophytoplankton that is calcifiers
         IF( iom_use( "LNnut"   ) )   CALL iom_put( "LNnut"  , xlimphy(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
         IF( iom_use( "LDnut"   ) )   CALL iom_put( "LDnut"  , xlimdia(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
@@ -231,10 +246,10 @@
       !!                  ***  ROUTINE p4z_lim_init  ***
       !!
-      !! ** Purpose :   Initialization of nutrient limitation parameters
-      !!
-      !! ** Method  :   Read the nampislim namelist and check the parameters
+      !! ** Purpose :   Initialization of the nutrient limitation parameters
+      !!
+      !! ** Method  :   Read the namp4zlim namelist and check the parameters
       !!      called at the first timestep (nittrc000)
       !!
-      !! ** input   :   Namelist nampislim
+      !! ** input   :   Namelist namp4zlim
       !!
       !!----------------------------------------------------------------------
@@ -252,8 +267,8 @@
       ENDIF
       !
-      REWIND( numnatp_ref )              ! Namelist nampislim in reference namelist : Pisces nutrient limitation parameters
+      REWIND( numnatp_ref )              ! Namelist namp4zlim in reference namelist : Pisces nutrient limitation parameters
       READ  ( numnatp_ref, namp4zlim, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zlim in reference namelist' )
-      REWIND( numnatp_cfg )              ! Namelist nampislim in configuration namelist : Pisces nutrient limitation parameters 
+      REWIND( numnatp_cfg )              ! Namelist namp4zlim in configuration namelist : Pisces nutrient limitation parameters 
       READ  ( numnatp_cfg, namp4zlim, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zlim in configuration namelist' )
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmeso.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 12537)
@@ -8,6 +8,7 @@
    !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Quota model for iron
    !!----------------------------------------------------------------------
-   !!   p4z_meso       : Compute the sources/sinks for mesozooplankton
-   !!   p4z_meso_init  : Initialization of the parameters for mesozooplankton
+   !!   p4z_meso        : Compute the sources/sinks for mesozooplankton
+   !!   p4z_meso_init   : Initialization of the parameters for mesozooplankton
+   !!   p4z_meso_alloc  : Allocate variables for mesozooplankton 
    !!----------------------------------------------------------------------
    USE oce_trc         ! shared variables between ocean and passive tracers
@@ -23,6 +24,7 @@
    PUBLIC   p4z_meso              ! called in p4zbio.F90
    PUBLIC   p4z_meso_init         ! called in trcsms_pisces.F90
-   PUBLIC   p4z_meso_alloc
-
+   PUBLIC   p4z_meso_alloc        ! called in trcini_pisces.F90
+
+   !! * Shared module variables
    REAL(wp), PUBLIC ::  part2        !: part of calcite not dissolved in mesozoo guts
    REAL(wp), PUBLIC ::  xpref2d      !: mesozoo preference for diatoms
@@ -46,6 +48,6 @@
    REAL(wp), PUBLIC ::  xfracmig     !: Fractional biomass of meso that performs DVM
    LOGICAL , PUBLIC ::  ln_dvm_meso  !: Boolean to activate DVM of mesozooplankton
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: depmig
-   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: kmig 
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: depmig  !: DVM of mesozooplankton : migration depth
+   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: kmig    !: Vertical indice of the the migration depth
 
    !!----------------------------------------------------------------------
@@ -61,4 +63,10 @@
       !!
       !! ** Purpose :   Compute the sources/sinks for mesozooplankton
+      !!                This includes ingestion and assimilation, flux feeding
+      !!                and mortality. We use a passive prey switching  
+      !!                parameterization.
+      !!                All living compartments smaller than mesozooplankton
+      !!                are potential preys of mesozooplankton as well as small
+      !!                sinking particles 
       !!
       !! ** Method  : - ???
@@ -68,13 +76,10 @@
       INTEGER  :: ji, jj, jk, jkt
       REAL(wp) :: zcompadi, zcompaph, zcompapoc, zcompaz, zcompam
-      REAL(wp) :: zgraze2 , zdenom, zdenom2
-      REAL(wp) :: zfact   , zfood, zfoodlim, zproport, zbeta
+      REAL(wp) :: zgraze2 , zdenom, zdenom2, zfact   , zfood, zfoodlim, zproport, zbeta
       REAL(wp) :: zmortzgoc, zfrac, zfracfe, zratio, zratio2, zfracal, zgrazcal
       REAL(wp) :: zepsherf, zepshert, zepsherq, zepsherv, zgrarsig, zgraztotc, zgraztotn, zgraztotf
       REAL(wp) :: zmigreltime, zprcaca, zmortz, zgrasrat, zgrasratn
-      REAL(wp) :: zrespz, ztortz, zgrazd, zgrazz, zgrazpof
-      REAL(wp) :: zgrazn, zgrazpoc, zgraznf, zgrazf
-      REAL(wp) :: zgrazfffp, zgrazfffg, zgrazffep, zgrazffeg
-      REAL(wp) :: zrum, zcodel, zargu, zval
+      REAL(wp) :: zrespz, ztortz, zgrazd, zgrazz, zgrazpof, zgrazn, zgrazpoc, zgraznf, zgrazf
+      REAL(wp) :: zgrazfffp, zgrazfffg, zgrazffep, zgrazffeg, zrum, zcodel, zargu, zval
       CHARACTER (len=25) :: charout
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing, zfezoo2
@@ -95,4 +100,5 @@
       !
       ! Diurnal vertical migration of mesozooplankton
+      ! Computation of the migration depth
       ! ---------------------------------------------
       IF (ln_dvm_meso) CALL p4z_meso_depmig
@@ -104,19 +110,29 @@
                zfact     = xstep * tgfunc2(ji,jj,jk) * zcompam
 
-               !  Respiration rates of both zooplankton
-               !  -------------------------------------
+               !  linear mortality of mesozooplankton
+               !  A michaelis menten modulation term is used to avoid extinction of 
+               !  mesozooplankton at very low food concentration. Mortality is
+ 
+               !  enhanced in low O2 waters
+               !  -----------------------------------------------------------------
                zrespz    = resrat2 * zfact * ( trb(ji,jj,jk,jpmes) / ( xkmort + trb(ji,jj,jk,jpmes) )  &
                &           + 3. * nitrfac(ji,jj,jk) )
 
-               !  Zooplankton mortality. A square function has been selected with
-               !  no real reason except that it seems to be more stable and may mimic predation
-               !  ---------------------------------------------------------------
+               !  Zooplankton quadratic mortality. A square function has been selected with
+               !  to mimic predation and disease (density dependent mortality). It also tends
+               !  to stabilise the model
+               !  -------------------------------------------------------------------------
                ztortz    = mzrat2 * 1.e6 * zfact * trb(ji,jj,jk,jpmes)  * (1. - nitrfac(ji,jj,jk) )
-               !
+
+               !   Computation of the abundance of the preys
+               !   A threshold can be specified in the namelist
+               !   --------------------------------------------
                zcompadi  = MAX( ( trb(ji,jj,jk,jpdia) - xthresh2dia ), 0.e0 )
                zcompaz   = MAX( ( trb(ji,jj,jk,jpzoo) - xthresh2zoo ), 0.e0 )
                zcompapoc = MAX( ( trb(ji,jj,jk,jppoc) - xthresh2poc ), 0.e0 )
                ! Size effect of nanophytoplankton on grazing : the smaller it is, the less prone
-               ! it is to predation by mesozooplankton
+               ! it is to predation by mesozooplankton. We use a quota dependant parameterization
+               ! as a low quota indicates oligotrophic conditions which are charatcerized by
+               ! small cells
                ! -------------------------------------------------------------------------------
                zcompaph  = MAX( ( trb(ji,jj,jk,jpphy) - xthresh2phy ), 0.e0 ) &
@@ -124,5 +140,10 @@
 
                ! Mesozooplankton grazing
-               ! ------------------------
+               ! The total amount of food is the sum of all preys accessible to mesozooplankton 
+               ! multiplied by their food preference
+               ! A threshold can be specified in the namelist (xthresh2). However, when food 
+               ! concentration is close to this threshold, it is decreased to avoid the 
+               ! accumulation of food in the mesozoopelagic domain
+               ! -------------------------------------------------------------------------------
                zfood     = xpref2d * zcompadi + xpref2z * zcompaz + xpref2n * zcompaph + xpref2c * zcompapoc 
                zfoodlim  = MAX( 0., zfood - MIN( 0.5 * zfood, xthresh2 ) )
@@ -131,4 +152,8 @@
                zgraze2   = grazrat2 * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpmes) * (1. - nitrfac(ji,jj,jk)) 
 
+               ! The grazing pressure on each prey is computed assuming passive switching. This 
+               ! is equivalent to assuming that mesozooplankton have an opportunistic feeding
+               ! behaviour.
+               ! -----------------------------------------------------------------------------
                zgrazd    = zgraze2  * xpref2d  * zcompadi  * zdenom2 
                zgrazz    = zgraze2  * xpref2z  * zcompaz   * zdenom2 
@@ -140,6 +165,7 @@
                zgrazpof  = zgrazpoc * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn)
 
-               !  Mesozooplankton flux feeding on GOC
-               !  ----------------------------------
+               !  Mesozooplankton flux feeding on GOC and POC. The feeding pressure
+               ! is proportional to the flux
+               !  ------------------------------------------------------------------
                zgrazffeg = grazflux  * xstep * wsbio4(ji,jj,jk)      &
                &           * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpgoc) * trb(ji,jj,jk,jpmes) &
@@ -150,11 +176,14 @@
                &           * (1. - nitrfac(ji,jj,jk))
                zgrazfffp = zgrazffep * trb(ji,jj,jk,jpsfe) / (trb(ji,jj,jk,jppoc) + rtrn)
-               !
+               
                zgraztotc = zgrazd + zgrazz + zgrazn + zgrazpoc + zgrazffep + zgrazffeg
-               ! Compute the proportion of filter feeders
+               ! Compute the proportion of filter feeders. It is assumed steady state.
+               ! ---------------------------------------------------------------------  
                zproport  = (zgrazffep + zgrazffeg)/(rtrn + zgraztotc)
+
                ! Compute fractionation of aggregates. It is assumed that 
                ! diatoms based aggregates are more prone to fractionation
                ! since they are more porous (marine snow instead of fecal pellets)
+               ! -----------------------------------------------------------------
                zratio    = trb(ji,jj,jk,jpgsi) / ( trb(ji,jj,jk,jpgoc) + rtrn )
                zratio2   = zratio * zratio
@@ -164,4 +193,5 @@
                zfracfe   = zfrac * trb(ji,jj,jk,jpbfe) / (trb(ji,jj,jk,jpgoc) + rtrn)
 
+               ! Flux feeding is multiplied by the fractional biomass of flux feeders
                zgrazffep = zproport * zgrazffep
                zgrazffeg = zproport * zgrazffeg
@@ -176,6 +206,12 @@
                zgrazing(ji,jj,jk) = zgraztotc
 
-               !    Mesozooplankton efficiency
-               !    --------------------------
+               ! Mesozooplankton efficiency. 
+               ! We adopt a formulation proposed by Mitra et al. (2007)
+               ! The gross growth efficiency is controled by the most limiting nutrient.
+               ! Growth is also further decreased when the food quality is poor. This is currently
+               ! hard coded : it can be decreased by up to 50% (zepsherq)
+               ! GGE can also be decreased when food quantity is high, zepsherf (Montagnes and 
+               ! Fulton, 2012)
+               ! -----------------------------------------------------------------------------------
                zgrasrat  =  ( zgraztotf + rtrn )/ ( zgraztotc + rtrn )
                zgrasratn =  ( zgraztotn + rtrn )/ ( zgraztotc + rtrn )
@@ -189,4 +225,6 @@
                ! ---------------------------------------------
                zmortz = ztortz + zrespz
+               ! Mortality induced by the upper trophic levels, ztortz, is allocated 
+               ! according to a infinite chain of predators (ANderson et al., 2013)
                zmortzgoc = unass2 / ( 1. - epsher2 ) * ztortz + zrespz
                tra(ji,jj,jk,jpmes) = tra(ji,jj,jk,jpmes) - zmortz + zepsherv * zgraztotc
@@ -208,4 +246,5 @@
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) - zgrazfffg - zfracfe
                ! Calcite remineralization due to zooplankton activity
+               ! part2 of the ingested calcite is dissolving in the acidic gut
                zfracal = trb(ji,jj,jk,jpcal) / (trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + rtrn )
                zgrazcal = (zgrazffeg + zgrazpoc) * (1. - part2) * zfracal
@@ -218,7 +257,6 @@
                tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * ( zgrazcal + zprcaca )
                tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) - zgrazcal + zprcaca
-                
-               ! Correct the fluxes for the effect of DVM
-               ! A fixed fraction of mesozooplankton is assumed to migrate
+ 
+               ! Computation of total excretion and egestion by mesozoo. 
                ! ---------------------------------------------------------
                zgrarem(ji,jj,jk) = zgraztotc * ( 1. - zepsherv - unass2 ) &
@@ -232,4 +270,8 @@
       END DO
 
+      ! Computation of the effect of DVM by mesozooplankton
+      ! This part is only activated if ln_dvm_meso is set to true
+      ! The parameterization has been published in Gorgues et al. (2019).
+      ! -----------------------------------------------------------------
       IF (ln_dvm_meso) THEN
          ALLOCATE( zgramigrem(jpi,jpj), zgramigref(jpi,jpj), zgramigpoc(jpi,jpj), zgramigpof(jpi,jpj) )
@@ -253,10 +295,11 @@
          END DO
 
-
+        ! Compute the amount of materials that will go into vertical migration
+        ! This fraction is sumed over the euphotic zone and is removed from 
+        ! the fluxes driven by mesozooplankton in the euphotic zone.
+        ! --------------------------------------------------------------------
          DO jk = 1, jpk
             DO jj = 1, jpj
                DO ji = 1, jpi
-
-                  !   Compute the amount of materials that will go into vertical migration
                   zmigreltime = (1. - zstrn(ji,jj))
                   IF ( gdept_n(ji,jj,jk) <= heup(ji,jj) ) THEN
@@ -279,4 +322,7 @@
          END DO
       
+         ! The inorganic and organic fluxes induced by migrating organisms are added at the 
+         ! the migration depth (corresponding indice is set by kmig)
+         ! --------------------------------------------------------------------------------
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -296,4 +342,5 @@
          DEALLOCATE( zstrn )
 
+      ! End of the ln_dvm_meso part
       ENDIF
 
@@ -302,4 +349,6 @@
             DO ji = 1, jpi
                !   Update the arrays TRA which contain the biological sources and sinks
+               !   This only concerns the variables which are affected by DVM (inorganic 
+               !   nutrients, DOC agands, and particulate organic carbon). 
                zgrarsig  = zgrarem(ji,jj,jk) * sigma2
                tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zgrarsig
@@ -324,4 +373,5 @@
       END DO
       !
+      ! Write the output
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
          ALLOCATE( zw3d(jpi,jpj,jpk) )
@@ -364,5 +414,5 @@
       !! ** Purpose :   Initialization of mesozooplankton parameters
       !!
-      !! ** Method  :   Read the nampismes namelist and check the parameters
+      !! ** Method  :   Read the namp4zmes namelist and check the parameters
       !!      called at the first timestep (nittrc000)
       !!
@@ -383,8 +433,8 @@
       ENDIF
       !
-      REWIND( numnatp_ref )              ! Namelist nampismes in reference namelist : Pisces mesozooplankton
+      REWIND( numnatp_ref )              ! Namelist namp4zmes in reference namelist : Pisces mesozooplankton
       READ  ( numnatp_ref, namp4zmes, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zmes in reference namelist' )
-      REWIND( numnatp_cfg )              ! Namelist nampismes in configuration namelist : Pisces mesozooplankton
+      REWIND( numnatp_cfg )              ! Namelist namp4zmes in configuration namelist : Pisces mesozooplankton
       READ  ( numnatp_cfg, namp4zmes, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zmes in configuration namelist' )
@@ -427,6 +477,4 @@
       !!      temperature and chlorophylle following the parameterization 
       !!      proposed by Bianchi et al. (2013)
-      !!
-      !! ** input   :   
       !!----------------------------------------------------------------------
       INTEGER  :: ji, jj, jk
@@ -447,5 +495,5 @@
       ! Compute the averaged values of oxygen, temperature over the domain 
       ! 150m to 500 m depth.
-      !
+      ! ------------------------------------------------------------------
       DO jk =1, jpk
          DO jj = 1, jpj
@@ -462,4 +510,7 @@
       END DO
 
+      ! Compute the difference between surface values and the mean values in the mesopelagic
+      ! domain
+      ! ------------------------------------------------------------------------------------
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -472,5 +523,4 @@
       ! Bianchi et al. (2013)
       ! -------------------------------------------------------------------
-      !
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -484,5 +534,4 @@
       ! Computation of the corresponding jk indice 
       ! ------------------------------------------
-      ! 
       DO jk = 1, jpk-1
          DO jj = 1, jpj
@@ -501,5 +550,4 @@
       ! that it falls above the OMZ
       ! -----------------------------------------------------------------------
-      !
       DO ji =1, jpi
          DO jj = 1, jpj
@@ -530,5 +578,4 @@
    END FUNCTION p4z_meso_alloc
 
-
    !!======================================================================
 END MODULE p4zmeso
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmicro.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 12537)
@@ -22,4 +22,5 @@
    PRIVATE
 
+   !! * Shared module variables
    PUBLIC   p4z_micro         ! called in p4zbio.F90
    PUBLIC   p4z_micro_init    ! called in trcsms_pisces.F90
@@ -54,4 +55,9 @@
       !!
       !! ** Purpose :   Compute the sources/sinks for microzooplankton
+      !!                This includes ingestion and assimilation, flux feeding
+      !!                and mortality. We use a passive prey switching  
+      !!                parameterization.
+      !!                All living compartments smaller than microzooplankton
+      !!                are potential preys of mesozooplankton
       !!
       !! ** Method  : - ???
@@ -87,20 +93,35 @@
                zfact   = xstep * tgfunc2(ji,jj,jk) * zcompaz
 
-               !  Respiration rates of both zooplankton
-               !  -------------------------------------
+               !  linear mortality of mesozooplankton
+               !  A michaelis menten modulation term is used to avoid extinction of 
+               !  microzooplankton at very low food concentrations. Mortality is 
+               !  enhanced in low O2 waters
+               !  -----------------------------------------------------------------
                zrespz = resrat * zfact * trb(ji,jj,jk,jpzoo) / ( xkmort + trb(ji,jj,jk,jpzoo) )  &
                   &   + resrat * zfact * 3. * nitrfac(ji,jj,jk)
 
-               !  Zooplankton mortality. A square function has been selected with
-               !  no real reason except that it seems to be more stable and may mimic predation.
-               !  ---------------------------------------------------------------
+               !  Zooplankton quadratic mortality. A square function has been selected with
+               !  to mimic predation and disease (density dependent mortality). It also tends
+               !  to stabilise the model
+               !  -------------------------------------------------------------------------
                ztortz = mzrat * 1.e6 * zfact * trb(ji,jj,jk,jpzoo) * (1. - nitrfac(ji,jj,jk))
 
+               !   Computation of the abundance of the preys
+               !   A threshold can be specified in the namelist
+               !   Diatoms have a specific treatment. WHen concentrations 
+               !   exceed a certain value, diatoms are suppposed to be too 
+               !   big for microzooplankton.
+               !   --------------------------------------------------------
                zcompadi  = MIN( MAX( ( trb(ji,jj,jk,jpdia) - xthreshdia ), 0.e0 ), xsizedia )
                zcompaph  = MAX( ( trb(ji,jj,jk,jpphy) - xthreshphy ), 0.e0 )
                zcompapoc = MAX( ( trb(ji,jj,jk,jppoc) - xthreshpoc ), 0.e0 )
                
-               !     Microzooplankton grazing
-               !     ------------------------
+               ! Microzooplankton grazing
+               ! The total amount of food is the sum of all preys accessible to mesozooplankton 
+               ! multiplied by their food preference
+               ! A threshold can be specified in the namelist (xthresh). However, when food 
+               ! concentration is close to this threshold, it is decreased to avoid the 
+               ! accumulation of food in the mesozoopelagic domain
+               ! -------------------------------------------------------------------------------
                zfood     = xprefn * zcompaph + xprefc * zcompapoc + xprefd * zcompadi
                zfoodlim  = MAX( 0. , zfood - min(xthresh,0.5*zfood) )
@@ -109,4 +130,8 @@
                zgraze    = grazrat * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpzoo) * (1. - nitrfac(ji,jj,jk))
 
+               ! The grazing pressure on each prey is computed assuming passive switching. This 
+               ! is equivalent to assuming that microzooplankton have an opportunistic feeding
+               ! behaviour.
+               ! -----------------------------------------------------------------------------
                zgrazp    = zgraze  * xprefn * zcompaph  * zdenom2 
                zgrazm    = zgraze  * xprefc * zcompapoc * zdenom2 
@@ -124,6 +149,12 @@
                zgrazing(ji,jj,jk) = zgraztotc
 
-               !    Various remineralization and excretion terms
-               !    --------------------------------------------
+               ! Microzooplankton efficiency. 
+               ! We adopt a formulation proposed by Mitra et al. (2007)
+               ! The gross growth efficiency is controled by the most limiting nutrient.
+               ! Growth is also further decreased when the food quality is poor. This is currently
+               ! hard coded : it can be decreased by up to 50% (zepsherq)
+               ! GGE can also be decreased when food quantity is high, zepsherf (Montagnes and 
+               ! Fulton, 2012)
+               ! -----------------------------------------------------------------------------
                zgrasrat  = ( zgraztotf + rtrn ) / ( zgraztotc + rtrn )
                zgrasratn = ( zgraztotn + rtrn ) / ( zgraztotc + rtrn )
@@ -157,6 +188,4 @@
                tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarsig
                tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zgrarsig
-               !   Update the arrays TRA which contain the biological sources and sinks
-               !   --------------------------------------------------------------------
                zmortz = ztortz + zrespz
                tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) - zmortz + zepsherv * zgraztotc 
@@ -174,5 +203,6 @@
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz - zgrazmf
                !
-               ! calcite production
+               ! Calcite remineralization due to zooplankton activity
+               ! part of the ingested calcite is dissolving in the acidic gut
                zprcaca = xfracal(ji,jj,jk) * zgrazp
                prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
@@ -224,8 +254,8 @@
       !! ** Purpose :   Initialization of microzooplankton parameters
       !!
-      !! ** Method  :   Read the nampiszoo namelist and check the parameters
+      !! ** Method  :   Read the namp4zzoo namelist and check the parameters
       !!                called at the first timestep (nittrc000)
       !!
-      !! ** input   :   Namelist nampiszoo
+      !! ** input   :   Namelist namp4zzoo
       !!
       !!----------------------------------------------------------------------
@@ -243,8 +273,8 @@
       ENDIF
       !
-      REWIND( numnatp_ref )              ! Namelist nampiszoo in reference namelist : Pisces microzooplankton
+      REWIND( numnatp_ref )              ! Namelist namp4zzoo in reference namelist : Pisces microzooplankton
       READ  ( numnatp_ref, namp4zzoo, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zzoo in reference namelist' )
-      REWIND( numnatp_cfg )              ! Namelist nampiszoo in configuration namelist : Pisces microzooplankton
+      REWIND( numnatp_cfg )              ! Namelist namp4zzoo in configuration namelist : Pisces microzooplankton
       READ  ( numnatp_cfg, namp4zzoo, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zzoo in configuration namelist' )
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmort.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 12537)
@@ -20,12 +20,12 @@
    PRIVATE
 
-   PUBLIC   p4z_mort    
-   PUBLIC   p4z_mort_init    
-
-   REAL(wp), PUBLIC ::   wchln    !:
-   REAL(wp), PUBLIC ::   wchld    !:
-   REAL(wp), PUBLIC ::   wchldm   !:
-   REAL(wp), PUBLIC ::   mpratn   !:
-   REAL(wp), PUBLIC ::   mpratd   !:
+   PUBLIC   p4z_mort           ! Called from p4zbio.F90 
+   PUBLIC   p4z_mort_init      ! Called from trcini_pisces.F90 
+
+   REAL(wp), PUBLIC ::   wchln    !: Quadratic mortality rate of nanophytoplankton
+   REAL(wp), PUBLIC ::   wchld    !: Quadratic mortality rate of diatoms
+   REAL(wp), PUBLIC ::   wchldm   !: Maximum quadratic mortality rate of diatoms
+   REAL(wp), PUBLIC ::   mpratn   !: Linear mortality rate of nanophytoplankton
+   REAL(wp), PUBLIC ::   mpratd   !: Linear mortality rate of diatoms
 
    !!----------------------------------------------------------------------
@@ -40,5 +40,5 @@
       !!                     ***  ROUTINE p4z_mort  ***
       !!
-      !! ** Purpose :   Calls the different subroutine to initialize and compute
+      !! ** Purpose :   Calls the different subroutine to compute
       !!                the different phytoplankton mortality terms
       !!
@@ -49,5 +49,4 @@
       !
       CALL p4z_nano            ! nanophytoplankton
-      !
       CALL p4z_diat            ! diatoms
       !
@@ -82,17 +81,19 @@
                !     to 0
                zsizerat = MIN(1., MAX( 0., (quotan(ji,jj,jk) - 0.2) / 0.3) ) * trb(ji,jj,jk,jpphy)
-               !     Squared mortality of Phyto similar to a sedimentation term during
-               !     blooms (Doney et al. 1996)
+
+               ! Quadratic mortality of nano due to aggregation during
+               ! blooms (Doney et al. 1996)
+               ! -----------------------------------------------------
                zrespp = wchln * 1.e6 * xstep * xdiss(ji,jj,jk) * zcompaph * zsizerat 
 
-               !     Phytoplankton mortality. This mortality loss is slightly
-               !     increased when nutrients are limiting phytoplankton growth
-               !     as observed for instance in case of iron limitation.
+               ! Phytoplankton linear mortality
+               ! A michaelis-menten like term is introduced to avoid 
+               ! extinction of nanophyto in highly limited areas
+               ! ----------------------------------------------------
                ztortp = mpratn * xstep * zcompaph / ( xkmort + trb(ji,jj,jk,jpphy) ) * zsizerat
 
                zmortp = zrespp + ztortp
 
-               !   Update the arrays TRA which contains the biological sources and sinks
-
+               ! Update the arrays TRA which contains the biological sources and sinks
                zfactfe = trb(ji,jj,jk,jpnfe)/(trb(ji,jj,jk,jpphy)+rtrn)
                zfactch = trb(ji,jj,jk,jpnch)/(trb(ji,jj,jk,jpphy)+rtrn)
@@ -100,8 +101,7 @@
                tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zmortp * zfactch
                tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zmortp * zfactfe
+               ! Production PIC particles due to mortality
                zprcaca = xfracal(ji,jj,jk) * zmortp
-               !
                prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
-               !
                zfracal = 0.5 * xfracal(ji,jj,jk)
                tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprcaca
@@ -146,8 +146,8 @@
       IF( ln_timing )   CALL timing_start('p4z_diat')
       !
-      !    Aggregation term for diatoms is increased in case of nutrient
-      !    stress as observed in reality. The stressed cells become more
-      !    sticky and coagulate to sink quickly out of the euphotic zone
-      !     ------------------------------------------------------------
+      ! Aggregation term for diatoms is increased in case of nutrient
+      ! stress as observed in reality. The stressed cells become more
+      ! sticky and coagulate to sink quickly out of the euphotic zone
+      ! -------------------------------------------------------------
 
       DO jk = 1, jpkm1
@@ -157,21 +157,21 @@
                zcompadi = MAX( ( trb(ji,jj,jk,jpdia) - 1e-9), 0. )
 
-               !    Aggregation term for diatoms is increased in case of nutrient
-               !    stress as observed in reality. The stressed cells become more
-               !    sticky and coagulate to sink quickly out of the euphotic zone
-               !     ------------------------------------------------------------
-               !  Phytoplankton respiration 
-               !     ------------------------
+               ! Aggregation term for diatoms is increased in case of nutrient
+               ! stress as observed in reality. The stressed cells become more
+               ! sticky and coagulate to sink quickly out of the euphotic zone
+               ! ------------------------------------------------------------
                zlim2   = xlimdia(ji,jj,jk) * xlimdia(ji,jj,jk)
                zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) 
                zrespp2 = 1.e6 * xstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
 
-               !     Phytoplankton mortality. 
-               !     ------------------------
+               ! Phytoplankton linear mortality
+               ! A michaelis-menten like term is introduced to avoid 
+               ! extinction of diatoms in highly limited areas
+               !  ---------------------------------------------------
                ztortp2 = mpratd * xstep * trb(ji,jj,jk,jpdia)  / ( xkmort + trb(ji,jj,jk,jpdia) ) * zcompadi 
 
                zmortp2 = zrespp2 + ztortp2
 
-               !   Update the arrays tra which contains the biological sources and sinks
+               ! Update the arrays tra which contains the biological sources and sinks
                !   ---------------------------------------------------------------------
                zfactch = trb(ji,jj,jk,jpdch) / ( trb(ji,jj,jk,jpdia) + rtrn )
@@ -210,8 +210,8 @@
       !! ** Purpose :   Initialization of phytoplankton parameters
       !!
-      !! ** Method  :   Read the nampismort namelist and check the parameters
+      !! ** Method  :   Read the namp4zmort namelist and check the parameters
       !!              called at the first timestep
       !!
-      !! ** input   :   Namelist nampismort
+      !! ** input   :   Namelist namp4zmort
       !!
       !!----------------------------------------------------------------------
@@ -227,8 +227,8 @@
       ENDIF
       !
-      REWIND( numnatp_ref )              ! Namelist nampismort in reference namelist : Pisces phytoplankton
+      REWIND( numnatp_ref )              ! Namelist namp4zmort in reference namelist : Pisces phytoplankton
       READ  ( numnatp_ref, namp4zmort, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zmort in reference namelist' )
-      REWIND( numnatp_cfg )              ! Namelist nampismort in configuration namelist : Pisces phytoplankton
+      REWIND( numnatp_cfg )              ! Namelist namp4zmort in configuration namelist : Pisces phytoplankton
       READ  ( numnatp_cfg, namp4zmort, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zmort in configuration namelist' )
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zpoc.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 12537)
@@ -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
 
 
@@ -49,6 +52,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 ???
@@ -79,5 +89,5 @@
       solgoc = 0.04/ 2.56 * 1./ ( 1.-50**(-0.04) )
 
-      ! Initialisation of temprary arrys
+      ! Initialisation of temporary arrys
       IF( ln_p4z ) THEN
          zremipoc(:,:,:) = xremip
@@ -92,4 +102,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)
@@ -97,8 +109,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(:,:,:)
@@ -108,9 +120,8 @@
               IF (tmask(ji,jj,jk) == 1.) THEN
                 zdep = hmld(ji,jj)
-                !
+
                 ! In the case of GOC, lability is constant in the mixed layer 
                 ! It is computed only below the mixed layer depth
                 ! ------------------------------------------------------------
-                !
                 IF( gdept_n(ji,jj,jk) > zdep ) THEN
                   alphat = 0.
@@ -153,9 +164,8 @@
                     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               &
                     &   * e3t_n(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk)   &
@@ -192,9 +202,12 @@
 
       IF( ln_p4z ) THEN
+         ! The standard PISCES part
          DO jk = 1, jpkm1
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! 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 * trb(ji,jj,jk,jpgoc)
@@ -204,5 +217,5 @@
                   zofer3 = zremig * solgoc * trb(ji,jj,jk,jpbfe)
 
-                  ! -------------------------------------
+                  ! update of the TRA arrays
                   tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zorem3(ji,jj,jk)
                   tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) - zorem2 - zorem3(ji,jj,jk)
@@ -216,9 +229,12 @@
          END DO
       ELSE
+         ! PISCES-QUOTA part
          DO jk = 1, jpkm1
             DO jj = 1, jpj
                DO ji = 1, jpi
-                   ! POC disaggregation by turbulence and bacterial activity. 
-                  ! --------------------------------------------------------
+                  ! POC disaggregation by bacterial activity. It is a function
+                  ! of the mean lability and of temperature. Shrinking of 
+                  ! particles due to the bacterial activity is also included 
+                  ! ----------------------------------------------------------
                   zremig = zremigoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
                   zopoc2 = zremig  * trb(ji,jj,jk,jpgoc)
@@ -229,5 +245,5 @@
                   zofer2 = xremipn / xremipc * zremig * trb(ji,jj,jk,jpbfe)
 
-                  ! -------------------------------------
+                  ! Update of the TRA arrays
                   tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zorem3(ji,jj,jk)
                   tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + solgoc * zopon2 
@@ -254,18 +270,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 jk = 1, jpkm1
         DO jj = 1, jpj
@@ -284,5 +298,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(:,:,:)
@@ -318,5 +333,4 @@
      ENDIF
 
-     ! -----------------------------------------------------------------------
      ! The lability parameterization is used here. The code is here 
      ! almost identical to what is done for big particles. The only difference
@@ -325,5 +339,4 @@
      ! should be determined before.
      ! -----------------------------------------------------------------------
-     !
      DO jk = 2, jpkm1
         DO jj = 1, jpj
@@ -334,13 +347,11 @@
                   alphat = 0.
                   remint = 0.
-                  !
                   ! the scale factors are corrected with temperature
                   zsizek1  = e3t_n(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
                   zsizek = e3t_n(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                  !
+
                   ! Special treatment of the level just below the MXL
                   ! See the comments in the GOC section
                   ! ---------------------------------------------------
-                  !
                   IF ( gdept_n(ji,jj,jk-1) <= zdep ) THEN
                     !
@@ -362,10 +373,9 @@
                     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               &
                     &   * e3t_n(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk)   &
@@ -404,14 +414,17 @@
 
      IF( ln_p4z ) THEN
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF (tmask(ji,jj,jk) == 1.) THEN
-                    ! POC disaggregation by turbulence and bacterial activity. 
-                    ! --------------------------------------------------------
+        ! Standard PISCES part
+        DO jk = 1, jpkm1
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+                 IF (tmask(ji,jj,jk) == 1.) THEN
+                    ! POC disaggregation by 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 * trb(ji,jj,jk,jppoc)
                     zofer           = zremip * trb(ji,jj,jk,jpsfe)
 
+                    ! Update of the TRA arrays
                     tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zorem
                     orem(ji,jj,jk)      = orem(ji,jj,jk) + zorem
@@ -420,31 +433,33 @@
                     tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zofer
                     zfolimi(ji,jj,jk)   = zfolimi(ji,jj,jk) + zofer
-                  ENDIF
-               END DO
-            END DO
-         END DO
+                 ENDIF
+              END DO
+           END DO
+        END DO
      ELSE
-       DO jk = 1, jpkm1
-          DO jj = 1, jpj
-             DO ji = 1, jpi
-                ! POC disaggregation by turbulence and bacterial activity. 
-                ! --------------------------------------------------------
-                zremip = zremipoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
-                zopoc  = zremip * trb(ji,jj,jk,jppoc)
-                orem(ji,jj,jk)  = orem(ji,jj,jk) + zopoc
-                zopon  = xremipn / xremipc * zremip * trb(ji,jj,jk,jppon)
-                zopop  = xremipp / xremipc * zremip * trb(ji,jj,jk,jppop)
-                zofer  = xremipn / xremipc * zremip * trb(ji,jj,jk,jpsfe)
-
-                tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zopoc
-                tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) - zopon
-                tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) - zopop
-                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zofer
-                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zopoc
-                tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zopon 
-                tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + zopop 
-                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer 
-                zfolimi(ji,jj,jk)   = zfolimi(ji,jj,jk) + zofer
-             END DO
+        ! PISCES QUOTA part
+        DO jk = 1, jpkm1
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+                 ! POC disaggregation by 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 * trb(ji,jj,jk,jppoc)
+                 orem(ji,jj,jk)  = orem(ji,jj,jk) + zopoc
+                 zopon  = xremipn / xremipc * zremip * trb(ji,jj,jk,jppon)
+                 zopop  = xremipp / xremipc * zremip * trb(ji,jj,jk,jppop)
+                 zofer  = xremipn / xremipc * zremip * trb(ji,jj,jk,jpsfe)
+
+                 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zopoc
+                 tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) - zopon
+                 tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) - zopop
+                 tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zofer
+                 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zopoc
+                 tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zopon 
+                 tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + zopop 
+                 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer 
+                 zfolimi(ji,jj,jk)   = zfolimi(ji,jj,jk) + zofer
+              END DO
            END DO
         END DO
@@ -454,7 +469,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
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 12537)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                         ***  MODULE p4zprod  ***
-   !! TOP :  Growth Rate of the two phytoplanktons groups 
+   !! TOP :  Growth Rate of the two phytoplankton groups of PISCES 
    !!======================================================================
    !! History :   1.0  !  2004     (O. Aumont) Original code
@@ -24,5 +24,5 @@
    PUBLIC   p4z_prod         ! called in p4zbio.F90
    PUBLIC   p4z_prod_init    ! called in trcsms_pisces.F90
-   PUBLIC   p4z_prod_alloc
+   PUBLIC   p4z_prod_alloc   ! called in trcini_pisces.F90
 
    REAL(wp), PUBLIC ::   pislopen     !:
@@ -40,5 +40,5 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotan   !: proxy of N quota in Nanophyto
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotad   !: proxy of N quota in diatomee
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotad   !: proxy of N quota in diatoms
    
    REAL(wp) ::   r1_rday    ! 1 / rday
@@ -57,8 +57,8 @@
       !!                     ***  ROUTINE p4z_prod  ***
       !!
-      !! ** Purpose :   Compute the phytoplankton production depending on
-      !!              light, temperature and nutrient availability
-      !!
-      !! ** Method  : - ???
+      !! ** Purpose :   Computes phytoplankton production depending on
+      !!                light, temperature and nutrient availability
+      !!                Computes also the uptake of Iron and Si as well 
+      !!                as the chlorophyll content of the cells
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, knt   !
@@ -95,5 +95,6 @@
       zmxl_fac(:,:,:) = 0._wp ; zmxl_chl(:,:,:) = 0._wp 
 
-      ! Computation of the optimal production
+      ! Computation of the maximimum production
+      ! Parameters are taken from Bissinger et al. (2008)
       zprmaxn(:,:,:) = 0.8_wp * r1_rday * tgfunc(:,:,:)
       zprmaxd(:,:,:) = zprmaxn(:,:,:)
@@ -114,4 +115,9 @@
 
       ! Impact of the day duration and light intermittency on phytoplankton growth
+      ! Intermittency is supposed to have a similar effect on production as 
+      ! day length. The correcting factor is zmxl_fac. zmxl_chl is the fractional
+      ! day length and is used to compute the mean PAR during daytime.
+      ! Formulation for the impact of day length on PP is from Thompson (1999)
+      ! ------------------------------------------------------------------------- 
       DO jk = 1, jpkm1
          DO jj = 1 ,jpj
@@ -132,8 +138,11 @@
       zprdia(:,:,:) = zprmaxd(:,:,:) * zmxl_fac(:,:,:)
 
-      ! Maximum light intensity
       WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
 
       ! Computation of the P-I slope for nanos and diatoms
+      ! The formulation proposed by Geider et al. (1997) has been modified 
+      ! to exclude the effect of nutrient limitation and temperature in the PI
+      ! curve following Vichi et al. (2007)
+      ! -----------------------------------------------------------------------
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -145,4 +154,6 @@
                   zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp
                   !
+                  ! The initial slope of the PI curve can be increased for nano
+                  ! to account for photadaptation, for instance in the DCM
                   zpislopeadn(ji,jj,jk) = pislopen * ( 1.+ zadap  * EXP( -0.25 * enano(ji,jj,jk) ) )  &
                   &                   * trb(ji,jj,jk,jpnch) /( trb(ji,jj,jk,jpphy) * 12. + rtrn)
@@ -179,4 +190,5 @@
 
       !  Computation of a proxy of the N/C ratio
+      !  Steady state is assumed
       !  ---------------------------------------
       DO jk = 1, jpkm1
@@ -199,9 +211,10 @@
 
                 IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                   !    Si/C of diatoms
-                   !    ------------------------
-                   !    Si/C increases with iron stress and silicate availability
-                   !    Si/C is arbitrariliy increased for very high Si concentrations
-                   !    to mimic the very high ratios observed in the Southern Ocean (silpot2)
+                   ! Si/C of diatoms
+                   ! ------------------------
+                   ! Si/C increases with iron stress and silicate availability (zsilfac)
+                   ! Si/C is arbitrariliy increased for very high Si concentrations
+                   ! to mimic the very high ratios observed in the Southern Ocean (zsilfac2)
+                   ! -----------------------------------------------------------------------
                   zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
                   zsilim = MIN( zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) )
@@ -219,7 +232,7 @@
       END DO
 
-      !  Mixed-layer effect on production 
-      !  Sea-ice effect on production
-
+      ! Sea-ice effect on production
+      ! No production is assumed below sea ice
+      ! -------------------------------------- 
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -231,14 +244,22 @@
       END DO
 
-      ! Computation of the various production terms 
+      ! Computation of the various production  and nutrient uptake terms
+      ! ---------------------------------------------------------------
       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. (C)
+                  !  production term of nanophyto. (C)
                   zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trb(ji,jj,jk,jpphy) * rfact2
+
+                  !  New production (uptake of NO3)
                   zpronewn(ji,jj,jk)  = zprorcan(ji,jj,jk)* xnanono3(ji,jj,jk) / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn )
                   !
                   zratio = trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * fecnm + rtrn )
+
+                  !  Iron uptake rates of nanophytoplankton. Upregulation  
+                  !  is parameterized at low iron concentrations. Typical 
+                  !  formulation used in quota formulations. Uptake is downregulated
+                  !  when the quota is close to the maximum quota 
                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
                   zprofen(ji,jj,jk) = fecnm * zprmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
@@ -246,8 +267,14 @@
                   &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concnfe(ji,jj,jk) )  &
                   &             * zmax * trb(ji,jj,jk,jpphy) * rfact2
-                  !  production terms for diatoms (C)
+                  !  production terms of diatoms (C)
                   zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trb(ji,jj,jk,jpdia) * rfact2
+
+                  ! New production (uptake of NO3)
                   zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk) / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn )
-                  !
+
+                  !  Iron uptake rates of nanophytoplankton. Upregulation  
+                  !  is parameterized at low iron concentrations. Typical 
+                  !  formulation used in quota formulations. Uptake is downregulated
+                  !  when the quota is close to the maximum quota 
                   zratio = trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * fecdm + rtrn )
                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
@@ -262,19 +289,25 @@
 
       ! Computation of the chlorophyll production terms
+      ! The parameterization is taken from Geider et al. (1997)
+      ! -------------------------------------------------------
       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 )
+                  !  production term for nanophyto. ( chlorophyll )
                   znanotot = enanom(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
                   zprod    = rday * zprorcan(ji,jj,jk) * zprchln(ji,jj,jk) * xlimphy(ji,jj,jk)
                   zprochln = chlcmin * 12. * zprorcan (ji,jj,jk)
+                  ! The maximum reachable Chl quota is modulated by temperature
+                  ! following Geider (1987)
                   chlcnm_n   = MIN ( chlcnm, ( chlcnm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem))) * (1. - 1.14 / 43.4 * 20.))
                   zprochln = zprochln + (chlcnm_n-chlcmin) * 12. * zprod / &
                                         & (  zpislopeadn(ji,jj,jk) * znanotot +rtrn)
-                  !  production terms for diatoms ( chlorophyll )
+                  !  production terms of diatoms ( chlorophyll )
                   zdiattot = ediatm(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
                   zprod    = rday * zprorcad(ji,jj,jk) * zprchld(ji,jj,jk) * xlimdia(ji,jj,jk)
                   zprochld = chlcmin * 12. * zprorcad(ji,jj,jk)
+                  ! The maximum reachable Chl quota is modulated by temperature
+                  ! following Geider (1987)
                   chlcdm_n   = MIN ( chlcdm, ( chlcdm / (1. - 1.14 / 43.4 * tsn(ji,jj,jk,jp_tem))) * (1. - 1.14 / 43.4 * 20.))
                   zprochld = zprochld + (chlcdm_n-chlcmin) * 12. * zprod / &
@@ -319,4 +352,9 @@
      END DO
      !
+     ! Production and uptake of ligands by phytoplankton. This part is activated 
+     ! when ln_ligand is set to .true. in the namelist. Ligand uptake is small 
+     ! and based on the FeL model by Morel et al. (2008) and on the study of
+     ! Shaked and Lis (2012)
+     ! -------------------------------------------------------------------------
      IF( ln_ligand ) THEN
          zpligprod1(:,:,:) = 0._wp    ;    zpligprod2(:,:,:) = 0._wp
@@ -328,7 +366,8 @@
                     zfeup    = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)
                     tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp    &
-                    &       - zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 2.E3 * (1.0 - plig(ji,jj,jk) ) )
+                    &       - zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
                     zpligprod1(ji,jj,jk) = zdocprod * ldocp
-                    zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) * lthet
+                    zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) &
+                    &                      + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
                  ENDIF
               END DO
@@ -373,9 +412,9 @@
           ENDIF
           IF( iom_use( "LPRODP" ) )  THEN
-              zw3d(:,:,:) = zpligprod1(:,:,:) * 1e9 * zfact * tmask(:,:,:)
+              zw3d(:,:,:) = zpligprod1(:,:,:) * 1e9 * zfact * tmask(:,:,:)  ! Ligand production by phytoplankton
               CALL iom_put( "LPRODP"  , zw3d )
           ENDIF
           IF( iom_use( "LDETP" ) )  THEN
-              zw3d(:,:,:) = zpligprod2(:,:,:) * 1e9 * zfact * tmask(:,:,:)
+              zw3d(:,:,:) = zpligprod2(:,:,:) * 1e9 * zfact * tmask(:,:,:)  ! Uptake of ligands by phytoplankton
               CALL iom_put( "LDETP"  , zw3d )
           ENDIF
@@ -392,8 +431,8 @@
           ENDIF
           IF( iom_use( "LNlight" ) .OR. iom_use( "LDlight" ) )  THEN
-              zw3d(:,:,:) = zprbio (:,:,:) / (zprmaxn(:,:,:) + rtrn) * tmask(:,:,:) ! light limitation term
+              zw3d(:,:,:) = zprbio (:,:,:) / (zprmaxn(:,:,:) + rtrn) * tmask(:,:,:) ! light limitation term of nanophytoplankton
               CALL iom_put( "LNlight"  , zw3d )
               !
-              zw3d(:,:,:) = zprdia (:,:,:) / (zprmaxd(:,:,:) + rtrn) * tmask(:,:,:)  ! light limitation term
+              zw3d(:,:,:) = zprdia (:,:,:) / (zprmaxd(:,:,:) + rtrn) * tmask(:,:,:)  ! light limitation term of diatoms
               CALL iom_put( "LDlight"  , zw3d )
           ENDIF
@@ -474,8 +513,8 @@
       !! ** Purpose :   Initialization of phytoplankton production parameters
       !!
-      !! ** Method  :   Read the nampisprod namelist and check the parameters
+      !! ** Method  :   Read the namp4zprod namelist and check the parameters
       !!      called at the first timestep (nittrc000)
       !!
-      !! ** input   :   Namelist nampisprod
+      !! ** input   :   Namelist namp4zprod
       !!----------------------------------------------------------------------
       INTEGER ::   ios   ! Local integer
@@ -491,8 +530,8 @@
       ENDIF
       !
-      REWIND( numnatp_ref )              ! Namelist nampisprod in reference namelist : Pisces phytoplankton production
+      REWIND( numnatp_ref )              ! Namelist namp4zprod in reference namelist : Pisces phytoplankton production
       READ  ( numnatp_ref, namp4zprod, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zprod in reference namelist' )
-      REWIND( numnatp_cfg )              ! Namelist nampisprod in configuration namelist : Pisces phytoplankton production
+      REWIND( numnatp_cfg )              ! Namelist namp4zprod in configuration namelist : Pisces phytoplankton production
       READ  ( numnatp_cfg, namp4zprod, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zprod in configuration namelist' )
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zrem.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zrem.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zrem.F90	(revision 12537)
@@ -3,4 +3,6 @@
    !!                         ***  MODULE p4zrem  ***
    !! TOP :   PISCES Compute remineralization/dissolution of organic compounds
+   !!         except for POC which is treated in p4zpoc.F90
+   !!         This module is common to both PISCES and PISCES-QUOTA
    !!=========================================================================
    !! History :   1.0  !  2004     (O. Aumont) Original code
@@ -17,5 +19,5 @@
    USE p4zche          !  chemical model
    USE p4zprod         !  Growth rate of the 2 phyto groups
-   USE p4zlim
+   USE p4zlim          !  Nutrient limitation terms
    USE prtctl_trc      !  print control for debugging
    USE iom             !  I/O manager
@@ -26,17 +28,18 @@
 
    PUBLIC   p4z_rem         ! called in p4zbio.F90
-   PUBLIC   p4z_rem_init    ! called in trcsms_pisces.F90
-   PUBLIC   p4z_rem_alloc
-
-   REAL(wp), PUBLIC ::   xremikc    !: remineralisation rate of DOC 
-   REAL(wp), PUBLIC ::   xremikn    !: remineralisation rate of DON 
-   REAL(wp), PUBLIC ::   xremikp    !: remineralisation rate of DOP 
-   REAL(wp), PUBLIC ::   xremik     !: remineralisation rate of POC 
+   PUBLIC   p4z_rem_init    ! called in trcini_pisces.F90
+   PUBLIC   p4z_rem_alloc   ! called in trcini_pisces.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC ::   xremikc    !: remineralisation rate of DOC (p5z) 
+   REAL(wp), PUBLIC ::   xremikn    !: remineralisation rate of DON (p5z) 
+   REAL(wp), PUBLIC ::   xremikp    !: remineralisation rate of DOP (p5z) 
+   REAL(wp), PUBLIC ::   xremik     !: remineralisation rate of DOC (p4z) 
    REAL(wp), PUBLIC ::   nitrif     !: NH4 nitrification rate 
-   REAL(wp), PUBLIC ::   xsirem     !: remineralisation rate of POC 
-   REAL(wp), PUBLIC ::   xsiremlab  !: fast remineralisation rate of POC 
+   REAL(wp), PUBLIC ::   xsirem     !: remineralisation rate of biogenic silica
+   REAL(wp), PUBLIC ::   xsiremlab  !: fast remineralisation rate of BSi
    REAL(wp), PUBLIC ::   xsilab     !: fraction of labile biogenic silica 
    REAL(wp), PUBLIC ::   feratb     !: Fe/C quota in bacteria
-   REAL(wp), PUBLIC ::   xkferb     !: Half-saturation constant for bacteria Fe/C
+   REAL(wp), PUBLIC ::   xkferb     !: Half-saturation constant for bacterial Fe/C
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   denitr   !: denitrification array
@@ -53,7 +56,13 @@
       !!                     ***  ROUTINE p4z_rem  ***
       !!
-      !! ** Purpose :   Compute remineralization/scavenging of organic compounds
-      !!
-      !! ** Method  : - ???
+      !! ** Purpose :   Compute remineralization/dissolution of organic compounds
+      !!                Computes also nitrification of ammonium 
+      !!                The solubilization/remineralization of POC is treated 
+      !!                in p4zpoc.F90. The dissolution of calcite is processed
+      !!                in p4zlys.F90. 
+      !!
+      !! ** Method  : - Bacterial biomass is computed implicitely based on a 
+      !!                parameterization developed from an explicit modeling
+      !!                of PISCES in an alternative version 
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, knt ! ocean time step
@@ -73,5 +82,5 @@
       IF( ln_timing )   CALL timing_start('p4z_rem')
       !
-      ! Initialisation of arrys
+      ! Initialisation of arrays
       zdepprod(:,:,:) = 1._wp
       zdepeff (:,:,:) = 0.3_wp
@@ -81,6 +90,5 @@
       zfacsi(:,:,:)   = xsilab
 
-      ! Computation of the mean phytoplankton concentration as
-      ! a crude estimate of the bacterial biomass
+      ! Computation of the mean bacterial concentration
       ! this parameterization has been deduced from a model version
       ! that was modeling explicitely bacteria
@@ -107,17 +115,24 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! DOC ammonification. Depends on depth, phytoplankton biomass
-                  ! and a limitation term which is supposed to be a parameterization of the bacterial activity. 
+                  ! DOC ammonification. Depends on a limitation term of the bacterial activity
+                  ! and on the implicit bacteria concentration
+                  ! --------------------------------------------------------------------------
                   zremik = xremik * xstep / 1.e-6 * xlimbac(ji,jj,jk) * zdepbac(ji,jj,jk) 
                   zremik = MAX( zremik, 2.74e-4 * xstep )
+
                   ! Ammonification in oxic waters with oxygen consumption
                   ! -----------------------------------------------------
                   zolimit = zremik * ( 1.- nitrfac(ji,jj,jk) ) * trb(ji,jj,jk,jpdoc) 
                   zolimi(ji,jj,jk) = MIN( ( trb(ji,jj,jk,jpoxy) - rtrn ) / o2ut, zolimit ) 
+
                   ! Ammonification in suboxic waters with denitrification
-                  ! -------------------------------------------------------
+                  ! -----------------------------------------------------
                   zammonic = zremik * nitrfac(ji,jj,jk) * trb(ji,jj,jk,jpdoc)
                   denitr(ji,jj,jk)  = zammonic * ( 1. - nitrfac2(ji,jj,jk) )
                   denitr(ji,jj,jk)  = MIN( ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenit, denitr(ji,jj,jk) )
+
+                  ! Ammonification in waters depleted in O2 and NO3 based on 
+                  ! other redox processes
+                  ! --------------------------------------------------------
                   zoxyremc          = zammonic - denitr(ji,jj,jk)
                   !
@@ -125,6 +140,5 @@
                   denitr (ji,jj,jk) = MAX( 0.e0, denitr (ji,jj,jk) )
                   zoxyremc          = MAX( 0.e0, zoxyremc )
-
-                  !
+                  ! Update of the TRA arrays
                   tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk) + zoxyremc
                   tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk) + zoxyremc
@@ -142,7 +156,8 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! DOC ammonification. Depends on depth, phytoplankton biomass
-                  ! and a limitation term which is supposed to be a parameterization of the bacterial activity. 
-                  ! -----------------------------------------------------------------
+                  ! DOC ammonification. Depends on a limitation term of the bacterial activity
+                  ! and on the implicit bacteria concentration
+                  ! ---------------------------------------------------------------
+-----------
                   zremik = xstep / 1.e-6 * MAX(0.01, xlimbac(ji,jj,jk)) * zdepbac(ji,jj,jk) 
                   zremik = MAX( zremik, 2.74e-4 * xstep / xremikc )
@@ -170,5 +185,5 @@
                   zoxyremn  = zremikn * zoxyremc * trb(ji,jj,jk,jpdon) / ( trb(ji,jj,jk,jpdoc) + rtrn )
                   zoxyremp  = zremikp * zoxyremc * trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdoc) + rtrn )
-
+                  ! Update of the TRA arrays
                   tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zolimip + zdenitrp + zoxyremp
                   tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zolimin + zdenitrn + zoxyremn
@@ -185,5 +200,4 @@
          !
       ENDIF
-
 
       DO jk = 1, jpkm1
@@ -220,8 +234,12 @@
                ! Bacteries are obliged to take up iron from the water. Some
                ! studies (especially at Papa) have shown this uptake to be significant
-               ! ----------------------------------------------------------
+               ! ---------------------------------------------------------------------
                zbactfer = feratb *  rfact2 * 0.6_wp / rday * tgfunc(ji,jj,jk) * xlimbacl(ji,jj,jk)     &
                   &              * trb(ji,jj,jk,jpfer) / ( xkferb + trb(ji,jj,jk,jpfer) )    &
                   &              * zdepeff(ji,jj,jk) * zdepbac(ji,jj,jk)
+
+               ! Only the transfer of iron from its dissolved form to particles
+               ! is treated here. The GGE of bacteria supposed to be equal to 
+               ! 0.33. This is hard-coded. 
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zbactfer*0.33
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zbactfer*0.25
@@ -242,14 +260,16 @@
       ! of bSi. Set to a constant in the upper ocean
       ! ---------------------------------------------------------------
-
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
+               ! Remineralization rate of BSi depedant on T and saturation
+               ! The parameterization is taken from Ridgwell et al. (2002) 
+               ! ---------------------------------------------------------
                zdep     = MAX( hmld(ji,jj), heup_01(ji,jj) )
                zsatur   = MAX( rtrn, ( sio3eq(ji,jj,jk) - trb(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn ) )
                zsatur2  = ( 1. + tsn(ji,jj,jk,jp_tem) / 400.)**37
                znusil   = 0.225  * ( 1. + tsn(ji,jj,jk,jp_tem) / 15.) * zsatur + 0.775 * zsatur2 * zsatur**9.25
-               ! Remineralization rate of BSi depedant on T and saturation
-               ! ---------------------------------------------------------
+               ! Computation of the vertical evolution of the labile fraction
+               ! of bSi. This is computed assuming steady state.
                IF ( gdept_n(ji,jj,jk) > zdep ) THEN
                   zfacsib(ji,jj,jk) = zfacsib(ji,jj,jk-1) * EXP( -0.5 * ( xsiremlab - xsirem )  &
@@ -261,5 +281,5 @@
                zsiremin = ( xsiremlab * zfacsi(ji,jj,jk) + xsirem * ( 1. - zfacsi(ji,jj,jk) ) ) * xstep * znusil
                zosil    = zsiremin * trb(ji,jj,jk,jpgsi)
-               !
+               ! Update of the TRA arrays
                tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) - zosil
                tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) + zosil
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 12537)
@@ -2,5 +2,8 @@
    !!======================================================================
    !!                         ***  MODULE p4sed  ***
-   !! TOP :   PISCES Compute loss of organic matter in the sediments
+   !! TOP :   PISCES Compute loss of biogenic matter in the sediments
+   !!                Compute gain/loss of tracers from dust, rivers and 
+   !!                sediments 
+   !!                This module is used both by PISCES and PISCES-QUOTA
    !!======================================================================
    !! History :   1.0  !  2004-03 (O. Aumont) Original code
@@ -8,7 +11,9 @@
    !!             3.4  !  2011-06 (C. Ethe) USE of fldread
    !!             3.5  !  2012-07 (O. Aumont) improvment of river input of nutrients 
-   !!----------------------------------------------------------------------
+   !!-----------------------------------------------------------------------
    !!   p4z_sed        :  Compute loss of organic matter in the sediments
-   !!----------------------------------------------------------------------
+   !!                  :  Compute gain/loss of tracers from dust, rivers and 
+   !!                     sediments 
+   !!-----------------------------------------------------------------------
    USE oce_trc         !  shared variables between ocean and passive tracers
    USE trc             !  passive tracers common variables 
@@ -43,7 +48,10 @@
       !!                     ***  ROUTINE p4z_sed  ***
       !!
-      !! ** Purpose :   Compute loss of organic matter in the sediments. This
-      !!              is by no way a sediment model. The loss is simply 
-      !!              computed to balance the inout from rivers and dust
+      !! ** Purpose :   Compute loss of biogenic matter in the sediments. This
+      !!              is by no way a real sediment model. The loss is simply 
+      !!              computed from metamodels.
+      !!              Loss/gain of tracers are also computed here for 
+      !!              dust, rivers, sediments and hydrothermal vents (Fe) 
+      !!              N2 fixation is modeled using an implicit approach
       !!
       !! ** Method  : - ???
@@ -93,5 +101,8 @@
       zsedc   (:,:) = 0.e0
 
-      ! Iron input/uptake due to sea ice : Crude parameterization based on Lancelot et al.
+      ! Iron input/uptake due to sea ice : Crude parameterization based on 
+      ! Lancelot et al. Iron concentration in sea-ice is constant and set 
+      ! in the namelist_pisces (icefeinput). ln_ironice is forced to false
+      ! when nn_ice_tr = 1
       ! ----------------------------------------------------
       IF( ln_ironice ) THEN  
@@ -99,4 +110,5 @@
          ALLOCATE( zironice(jpi,jpj) )
          !                                              
+         ! Compute the iron flux between sea ice and sea water
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -106,5 +118,5 @@
             END DO
          END DO
-         !
+         ! Update of the TRA array
          tra(:,:,1,jpfer) = tra(:,:,1,jpfer) + zironice(:,:) 
          ! 
@@ -121,5 +133,9 @@
          !                                              
          ALLOCATE( zsidep(jpi,jpj), zpdep(jpi,jpj,jpk), zirondep(jpi,jpj,jpk) )
-         !                                              ! Iron and Si deposition at the surface
+         ! Iron, P and Si deposition at the surface
+         ! Iron flux at the surface due to dust deposition. Solubility can be 
+         ! be variable if ln_solub is set to true. In that case, solubility 
+         ! has to be provided in the specific input file (read in p4zsbc)
+         ! ------------------------------------------------------------------
          IF( ln_solub ) THEN
             zirondep(:,:,1) = solub(:,:) * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss 
@@ -127,14 +143,23 @@
             zirondep(:,:,1) = dustsolub  * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss 
          ENDIF
+         ! Si and P flux at the surface due to dust deposition. The content 
+         ! and the solubility are hard coded
+         ! ----------------------------------------------------------------
          zsidep(:,:)   = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 28.1 
          zpdep (:,:,1) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 31. / po4r 
-         !                                              ! Iron solubilization of particles in the water column
-         !                                              ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/j
-         zwdust = 0.03 * rday / ( wdust * 55.85 ) / ( 270. * rday )
+         ! Iron solubilization of particles in the water column
+         ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/d
+         ! Dust are supposed to sink at wdust sinking speed. 3% of the iron 
+         ! in dust is hypothesized to be soluble at a dissolution rate set to 
+         ! 1/(250 days). The vertical distribution of iron in dust is computed 
+         ! from a steady state assumption. Parameters are very uncertain and 
+         ! are estimated from the literature quoted in Raiswell et al. (2011) 
+         ! ------------------------------------------------------------------- 
+         zwdust = 0.03 * rday / ( wdust * 55.85 ) / ( 250. * rday )
          DO jk = 2, jpkm1
-            zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -gdept_n(:,:,jk) / 540. )
+            zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -gdept_n(:,:,jk) / (250. * wdust) )
             zpdep   (:,:,jk) = zirondep(:,:,jk) * 0.023
          END DO
-         !                                              ! Iron solubilization of particles in the water column
+         ! Solubilization of particles in the water column (Si, P, Fe)
          tra(:,:,1,jpsil) = tra(:,:,1,jpsil) + zsidep  (:,:)
          DO jk = 1, jpkm1
@@ -156,5 +181,5 @@
      
       ! Add the external input of nutrients from river
-      ! ----------------------------------------------------------
+      ! ----------------------------------------------
       IF( ln_river ) THEN
          DO jj = 1, jpj
@@ -171,4 +196,8 @@
             ENDDO
          ENDDO
+         ! When prognostic ligands are activated, ligands are supplied 
+         ! to the ocean by rivers. We assume that the amount of ligands
+         ! is equal to that of iron (iron is completely complexed)
+         ! ------------------------------------------------------------
          IF (ln_ligand) THEN
             DO jj = 1, jpj
@@ -180,4 +209,5 @@
             ENDDO
          ENDIF
+         ! PISCES-QUOTA part
          IF( ln_p5z ) THEN
             DO jj = 1, jpj
@@ -200,5 +230,6 @@
 
       ! Add the external input of iron from hydrothermal vents
-      ! ------------------------------------------------------
+      ! Please refer to Tagliabue et al. (2010) for more information
+      ! ------------------------------------------------------------
       IF( ln_hydrofe ) THEN
             tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2
@@ -211,5 +242,7 @@
       ENDIF
 
-      ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments
+      ! OA: Warning, the following part is necessary to avoid CFL problems 
+      ! above the sediments. Vertical sinking speed is limited using the 
+      ! typical CFL criterion
       ! --------------------------------------------------------------------
       DO jj = 1, jpj
@@ -222,4 +255,5 @@
       END DO
       !
+      ! No sediment module activated
       IF( .NOT.lk_sed ) THEN
 !
@@ -233,6 +267,8 @@
          ENDIF
 
-         ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used
-         ! Computation of the fraction of organic matter that is permanently buried from Dunne's model
+         ! Computation of the sediment denitrification proportion: The metamodel 
+         ! from Middleburg (2006) is used
+         ! Computation of the fraction of organic matter that is permanently 
+         ! buried from Dunne's model (2007)
          ! -------------------------------------------------------
          DO jj = 1, jpj
@@ -259,9 +295,11 @@
       ENDIF
 
-      ! This loss is scaled at each bottom grid cell for equilibrating the total budget of silica in the ocean.
-      ! Thus, the amount of silica lost in the sediments equal the supply at the surface (dust+rivers)
-      ! ------------------------------------------------------
+      ! Fraction of dSi that is remineralized in the sediments. This is 
+      ! set so that the burial in sediments equals the total input of Si
+      ! by rivers and dust (sedsilfrac)
+      ! ----------------------------------------------------------------
       IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac
 
+      ! Loss of bSi and CaCO3 to the sediments
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -278,4 +316,9 @@
       !
       IF( .NOT.lk_sed ) THEN
+         ! Dissolution of CaCO3 and bSi in the sediments. This is 
+         ! instantaneous since here sediments are not explicitly 
+         ! modeled. The amount of CaCO3 that dissolves in the sediments
+         ! is computed using a metamodel constructed from Archer (1996)
+         ! ------------------------------------------------------------
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -298,4 +341,5 @@
       ENDIF
       !
+      ! Loss of particulate organic carbon and Fe to the sediments
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -311,4 +355,5 @@
       END DO
       !
+      ! Loss of particulate organic N and P to the sediments (p5z)
       IF( ln_p5z ) THEN
          DO jj = 1, jpj
@@ -327,6 +372,9 @@
 
       IF( .NOT.lk_sed ) THEN
+         ! Degradation of organic matter in the sediments. The metamodel of 
+         ! Middleburg (2006) is used here to mimic the diagenetic reactions. 
          ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after
          ! denitrification in the sediments. Not very clever, but simpliest option.
+         ! ------------------------------------------------------------------------
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -360,19 +408,26 @@
 
 
-      ! Nitrogen fixation process
-      ! Small source iron from particulate inorganic iron
-      !-----------------------------------
+      ! Nitrogen fixation process : light limitation of diazotrophy
+      ! Small source of iron from particulate inorganic iron (photochemistry)
+      !----------------------------------------------------------------------
       DO jk = 1, jpkm1
          zlight (:,:,jk) =  ( 1.- EXP( -etot_ndcy(:,:,jk) / diazolight ) ) * ( 1. - fr_i(:,:) ) 
          zsoufer(:,:,jk) = zlight(:,:,jk) * 2E-11 / ( 2E-11 + biron(:,:,jk) )
       ENDDO
+
+      ! Diazotrophy (nitrogen fixation) is modeled according to an empirical
+      ! formulation. This is described in Aumont et al. (2015). Limitation 
+      ! by P and Fe is computed. Inhibition by high N concentrations is imposed.
+      ! Diazotrophy sensitivity to temperature is parameterized as in 
+      ! Ye et al. (2012)  
+      ! ------------------------------------------------------------------------
       IF( ln_p4z ) THEN
+         ! PISCES part
          DO jk = 1, jpkm1
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  !                      ! Potential nitrogen fixation dependant on temperature and iron
+                  ! Potential nitrogen fixation dependant on temperature and iron
                   ztemp = tsn(ji,jj,jk,jp_tem)
                   zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
-                  !       Potential nitrogen fixation dependant on temperature and iron
                   xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) )
                   xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
@@ -388,11 +443,11 @@
          END DO
       ELSE       ! p5z
+         ! PISCES-QUOTA part
          DO jk = 1, jpkm1
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  !                      ! Potential nitrogen fixation dependant on temperature and iron
+                  ! Potential nitrogen fixation dependant on temperature and iron
                   ztemp = tsn(ji,jj,jk,jp_tem)
                   zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
-                  !       Potential nitrogen fixation dependant on temperature and iron
                   xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) )
                   xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
@@ -410,7 +465,8 @@
       ENDIF
 
-      ! Nitrogen change due to nitrogen fixation
-      ! ----------------------------------------
+      ! Update of the TRA arrays due to nitrogen fixation
+      ! -------------------------------------------------
       IF( ln_p4z ) THEN
+         ! PISCES part
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -428,10 +484,9 @@
                   tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
-                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trb(ji,jj,jk,jppo4) ) &
-                  &                     * 0.001 * trb(ji,jj,jk,jpdoc) * xstep
               END DO
             END DO 
          END DO
       ELSE    ! p5z
+         ! PISCES-QUOTA part
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -475,8 +530,8 @@
                CALL iom_put( "INTNFIX" , zwork ) 
             ENDIF
-            IF( iom_use("SedCal" ) ) CALL iom_put( "SedCal", zsedcal(:,:) * zfact )
-            IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * zfact )
-            IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * zfact )
-            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 )
+            IF( iom_use("SedCal" ) ) CALL iom_put( "SedCal", zsedcal(:,:) * zfact ) ! Permanent burial of CaCO3 in sediments
+            IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * zfact ) ! Permanent burial of bSi in sediments
+            IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * zfact ) ! Permanent burial of OC in sediments
+            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 ) ! Denitrification in the sediments
          ENDIF
       ENDIF
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsink.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 12537)
@@ -2,5 +2,7 @@
    !!======================================================================
    !!                         ***  MODULE p4zsink  ***
-   !! TOP :  PISCES  vertical flux of particulate matter due to gravitational sinking
+   !! TOP :  PISCES  vertical flux of particulate matter due to 
+   !!        gravitational sinking
+   !!        This module is the same for both PISCES and PISCES-QUOTA
    !!======================================================================
    !! History :   1.0  !  2004     (O. Aumont) Original code
@@ -8,4 +10,6 @@
    !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Change aggregation formula
    !!             3.5  !  2012-07  (O. Aumont) Introduce potential time-splitting
+   !!             4.0  !  2019     (O. Aumont) an external subroutine is called
+   !!                                          to compute the impact of sinking
    !!----------------------------------------------------------------------
    !!   p4z_sink       :  Compute vertical flux of particulate matter due to gravitational sinking
@@ -25,11 +29,11 @@
 
    PUBLIC   p4z_sink         ! called in p4zbio.F90
-   PUBLIC   p4z_sink_init    ! called in trcsms_pisces.F90
-   PUBLIC   p4z_sink_alloc
+   PUBLIC   p4z_sink_init    ! called in trcini_pisces.F90
+   PUBLIC   p4z_sink_alloc   ! called in trcini_pisces.F90
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinking, sinking2  !: POC sinking fluxes 
    !                                                          !  (different meanings depending on the parameterization)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkingn, sinking2n  !: POC sinking fluxes 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkingp, sinking2p  !: POC sinking fluxes 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkingn, sinking2n  !: PON sinking fluxes 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkingp, sinking2p  !: POP sinking fluxes 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkcal, sinksil   !: CaCO3 and BSi sinking fluxes
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer            !: Small BFe sinking fluxes
@@ -45,8 +49,4 @@
 CONTAINS
 
-   !!----------------------------------------------------------------------
-   !!   'standard sinking parameterisation'                  ???
-   !!----------------------------------------------------------------------
-
    SUBROUTINE p4z_sink ( kt, knt )
       !!---------------------------------------------------------------------
@@ -54,7 +54,10 @@
       !!
       !! ** Purpose :   Compute vertical flux of particulate matter due to 
-      !!                gravitational sinking
-      !!
-      !! ** Method  : - ???
+      !!                gravitational sinking. 
+      !!
+      !! ** Method  : - An external advection subroutine is called to compute
+      !!                the impact of sinking on the particles. The tracers
+      !!                concentrations are updated in this subroutine which
+      !!                is mandatory to deal with negative concentrations
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) :: kt, knt
@@ -75,8 +78,11 @@
       consgoc(:,:,:) = 0.
 
-      !
-      !    Sinking speeds of detritus is increased with depth as shown
-      !    by data and from the coagulation theory
-      !    -----------------------------------------------------------
+      ! Sinking speeds of big detritus is increased with depth as shown
+      ! by data and from the coagulation theory. This is controled by
+      ! wsbio2max and wsbio2scale. If wsbio2max is set to wsbio2, then
+      ! sinking speed is constant with depth.
+      ! CaCO3 and bSi are supposed to sink at the big particles speed 
+      ! due to their high density
+      ! ---------------------------------------------------------------
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -89,10 +95,9 @@
       END DO
 
-      ! limit the values of the sinking speeds to avoid numerical instabilities  
+      ! Sinking speed of the small particles is constant
       wsbio3(:,:,:) = wsbio
 
-      !
-      !  Initializa to zero all the sinking arrays 
-      !   -----------------------------------------
+      ! Initialize to zero all the sinking arrays 
+      ! -----------------------------------------
       sinking (:,:,:) = 0.e0
       sinking2(:,:,:) = 0.e0
@@ -102,6 +107,6 @@
       sinkfer2(:,:,:) = 0.e0
 
-      !   Compute the sedimentation term using p4zsink2 for all the sinking particles
-      !   -----------------------------------------------------
+      ! Compute the sedimentation term using trc_sink for all the sinking particles
+      ! ---------------------------------------------------------------------------
       CALL trc_sink( kt, wsbio3, sinking , jppoc, rfact2 )
       CALL trc_sink( kt, wsbio3, sinkfer , jpsfe, rfact2 )
@@ -111,4 +116,5 @@
       CALL trc_sink( kt, wsbio4, sinkcal , jpcal, rfact2 )
 
+      ! PISCES-QUOTA part
       IF( ln_p5z ) THEN
          sinkingn (:,:,:) = 0.e0
@@ -117,6 +123,6 @@
          sinking2p(:,:,:) = 0.e0
 
-         !   Compute the sedimentation term using p4zsink2 for all the sinking particles
-         !   -----------------------------------------------------
+         ! Compute the sedimentation term using trc_sink for all the sinking particles
+         ! ---------------------------------------------------------------------------
          CALL trc_sink( kt, wsbio3, sinkingn , jppon, rfact2 )
          CALL trc_sink( kt, wsbio3, sinkingp , jppop, rfact2 )
@@ -186,4 +192,10 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE p4z_sink_init  ***
+      !!
+      !! ** Purpose :   Initialization of sinking parameters
+      !!
+      !! ** Method  :   
+      !!
+      !! ** input   :   
       !!----------------------------------------------------------------------
       INTEGER :: jk
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 12537)
@@ -30,9 +30,9 @@
    PRIVATE
 
-   PUBLIC   p4z_sms_init   ! called in p4zsms.F90
-   PUBLIC   p4z_sms        ! called in p4zsms.F90
+   PUBLIC   p4z_sms_init   ! called in trcini_pisces.F90
+   PUBLIC   p4z_sms        ! called in trcsms_pisces.F90
 
    INTEGER ::    numco2, numnut, numnit      ! logical unit for co2 budget
-   REAL(wp) ::   alkbudget, no3budget, silbudget, ferbudget, po4budget
+   REAL(wp) ::   alkbudget, no3budget, silbudget, ferbudget, po4budget ! total budget of the different conservative elements
    REAL(wp) ::   xfact1, xfact2, xfact3
 
@@ -53,7 +53,8 @@
       !!              routines of PISCES bio-model
       !!
-      !! ** Method  : - at each new day ...
-      !!              - several calls of bio and sed ???
-      !!              - ...
+      !! ** Method  : - calls the various SMS subroutines
+      !!              - calls the sediment module (if ln_sediment)  
+      !!              - several calls of bio and sed (possible time-splitting)
+      !!              - handles the potential negative concentrations (xnegtr)
       !!---------------------------------------------------------------------
       !
@@ -72,5 +73,5 @@
         !
         IF( .NOT. ln_rsttr ) THEN
-            CALL p4z_che                              ! initialize the chemical constants
+            CALL p4z_che            ! initialize the chemical constants
             CALL ahini_for_at(hi)   !  set PH at kt=nit000
             t_oce_co2_flx_cum = 0._wp
@@ -110,9 +111,16 @@
       DO jnt = 1, nrdttrc          ! Potential time splitting if requested
          !
-         CALL p4z_bio( kt, jnt )   ! Biology
-         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation
+         CALL p4z_bio( kt, jnt )   ! Biological SMS
+         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation and dissolution
          CALL p4z_sed( kt, jnt )   ! Surface and Bottom boundary conditions
-         CALL p4z_flx( kt, jnt )   ! Compute surface fluxes
-         !
+         CALL p4z_flx( kt, jnt )   ! Compute surface air-sea fluxes
+         !
+         ! Handling of the negative concentrations
+         ! The biological SMS may generate negative concentrations
+         ! Trends are tested at each grid cell. If a negative concentrations 
+         ! is created at a grid cell, all the sources and sinks at that grid 
+         ! cell are scale to avoid that negative concentration. This approach 
+         ! is quite simplistic but it conserves mass.
+         ! ------------------------------------------------------------------
          xnegtr(:,:,:) = 1.e0
          DO jn = jp_pcs0, jp_pcs1
@@ -130,8 +138,10 @@
          !                                ! where at least 1 tracer concentration becomes negative
          !                                ! 
+         ! Concentrations are updated
          DO jn = jp_pcs0, jp_pcs1
            trb(:,:,:,jn) = trb(:,:,:,jn) + xnegtr(:,:,:) * tra(:,:,:,jn)
          END DO
         !
+        ! Trends are are reset to 0
          DO jn = jp_pcs0, jp_pcs1
             tra(:,:,:,jn) = 0._wp
@@ -153,4 +163,5 @@
 #endif
       !
+      ! If ln_sediment is set to .true. then the sediment module is called
       IF( ln_sediment ) THEN 
          !
@@ -181,8 +192,8 @@
       !!                     ***  p4z_sms_init  ***  
       !!
-      !! ** Purpose :   read PISCES namelist
-      !!
-      !! ** input   :   file 'namelist.trc.s' containing the following
-      !!             namelist: natext, natbio, natsms
+      !! ** Purpose :   read the general PISCES namelist
+      !!
+      !! ** input   :   file 'namelist_pisces' containing the following
+      !!                namelist: nampisbio, nampisdmp, nampismass 
       !!----------------------------------------------------------------------
       INTEGER :: ios                 ! Local integer output status for namelist read
@@ -268,5 +279,5 @@
       !!                   ***  ROUTINE p4z_rst  ***
       !!
-      !!  ** Purpose : Read or write variables in restart file:
+      !!  ** Purpose : Read or write specific PISCES variables in restart file:
       !!
       !!  WRITE(READ) mode:
@@ -280,8 +291,11 @@
       IF( TRIM(cdrw) == 'READ' ) THEN
          !
+         ! Read the specific variable of PISCES
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' p4z_rst : Read specific variables from pisces model '
          IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
          ! 
+         ! Read the pH. If not in the restart file, then it is initialized from
+         ! the initial conditions
          IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN
             CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
@@ -290,4 +304,6 @@
             CALL ahini_for_at(hi)
          ENDIF
+
+         ! Read the Si half saturation constant and the maximum Silica concentration
          CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) )
          IF( iom_varid( numrtr, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN
@@ -296,5 +312,6 @@
             xksimax(:,:) = xksi(:,:)
          ENDIF
-         !
+
+         ! Read the cumulative total flux. If not in the restart file, it is set to 0          
          IF( iom_varid( numrtr, 'tcflxcum', ldstop = .FALSE. ) > 0 ) THEN  ! cumulative total flux of carbon
             CALL iom_get( numrtr, 'tcflxcum' , t_oce_co2_flx_cum  )
@@ -303,5 +320,8 @@
          ENDIF
          !
+         ! PISCES-QUOTA specific part
          IF( ln_p5z ) THEN
+            ! Read the size of the different phytoplankton groups
+            ! If not in the restart file, they are set to 1
             IF( iom_varid( numrtr, 'sized', ldstop = .FALSE. ) > 0 ) THEN
                CALL iom_get( numrtr, jpdom_autoglo, 'sizep' , sizep(:,:,:)  )
@@ -316,4 +336,5 @@
         !
       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
+         ! write the specific variables of PISCES
          IF( kt == nitrst ) THEN
             IF(lwp) WRITE(numout,*)
@@ -321,12 +342,12 @@
             IF(lwp) WRITE(numout,*) '~~~~~~~'
          ENDIF
-         CALL iom_rstput( kt, nitrst, numrtw, 'PH', hi(:,:,:) )
-         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )
-         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) )
-         CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum )
+         CALL iom_rstput( kt, nitrst, numrtw, 'PH', hi(:,:,:) )           ! pH
+         CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )    ! Si 1/2 saturation constant
+         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) ) ! Si max concentration
+         CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum ) ! Cumulative CO2 flux
          IF( ln_p5z ) THEN
-            CALL iom_rstput( kt, nitrst, numrtw, 'sizep', sizep(:,:,:) )
-            CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sizen(:,:,:) )
-            CALL iom_rstput( kt, nitrst, numrtw, 'sized', sized(:,:,:) )
+            CALL iom_rstput( kt, nitrst, numrtw, 'sizep', sizep(:,:,:) )  ! Size of picophytoplankton
+            CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sizen(:,:,:) )  ! Size of nanophytoplankton
+            CALL iom_rstput( kt, nitrst, numrtw, 'sized', sized(:,:,:) )  ! Size of diatoms
          ENDIF
       ENDIF
@@ -339,5 +360,8 @@
       !!                    ***  p4z_dmp  ***
       !!
-      !! ** purpose  : Relaxation of some tracers
+      !! ** purpose  : Relaxation of the total budget of some elements
+      !!               This routine avoids the model to drift far from the 
+      !!               observed content in various elements
+      !!               Elements that may be relaxed : Alk, P, N, Si
       !!----------------------------------------------------------------------
       !
@@ -345,5 +369,5 @@
       !
       REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
-      REAL(wp) ::  po4mean = 2.165     ! mean value of phosphates
+      REAL(wp) ::  po4mean = 2.165     ! mean value of phosphate
       REAL(wp) ::  no3mean = 30.90     ! mean value of nitrate
       REAL(wp) ::  silmean = 91.51     ! mean value of silicate
@@ -368,13 +392,17 @@
             zsilsumn = glob_sum( 'p4zsms', trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
  
+            ! Correct the trn mean content of alkalinity
             IF(lwp) WRITE(numout,*) '       TALKN mean : ', zalksumn
             trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksumn
 
+            ! Correct the trn mean content of PO4
             IF(lwp) WRITE(numout,*) '       PO4N  mean : ', zpo4sumn
             trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sumn
 
+            ! Correct the trn mean content of NO3
             IF(lwp) WRITE(numout,*) '       NO3N  mean : ', zno3sumn
             trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sumn
 
+            ! Correct the trn mean content of SiO3
             IF(lwp) WRITE(numout,*) '       SiO3N mean : ', zsilsumn
             trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsumn )
@@ -388,13 +416,17 @@
  
                IF(lwp) WRITE(numout,*) ' '
+               ! Correct the trb mean content of alkalinity
                IF(lwp) WRITE(numout,*) '       TALKB mean : ', zalksumb
                trb(:,:,:,jptal) = trb(:,:,:,jptal) * alkmean / zalksumb
 
+               ! Correct the trb mean content of PO4
                IF(lwp) WRITE(numout,*) '       PO4B  mean : ', zpo4sumb
                trb(:,:,:,jppo4) = trb(:,:,:,jppo4) * po4mean / zpo4sumb
 
+               ! Correct the trb mean content of NO3
                IF(lwp) WRITE(numout,*) '       NO3B  mean : ', zno3sumb
                trb(:,:,:,jpno3) = trb(:,:,:,jpno3) * no3mean / zno3sumb
 
+               ! Correct the trb mean content of SiO3
                IF(lwp) WRITE(numout,*) '       SiO3B mean : ', zsilsumb
                trb(:,:,:,jpsil) = MIN( 400.e-6,trb(:,:,:,jpsil) * silmean / zsilsumb )
@@ -435,6 +467,6 @@
       ENDIF
 
+      ! Compute the budget of NO3
       IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
-         !   Compute the budget of NO3, ALK, Si, Fer
          IF( ln_p4z ) THEN
             zwork(:,:,:) =    trn(:,:,:,jpno3) + trn(:,:,:,jpnh4)                      &
@@ -454,4 +486,5 @@
       ENDIF
       !
+      ! Compute the budget of PO4
       IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          IF( ln_p4z ) THEN
@@ -472,4 +505,5 @@
       ENDIF
       !
+      ! Compute the budget of SiO3
       IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          zwork(:,:,:) =  trn(:,:,:,jpsil) + trn(:,:,:,jpgsi) + trn(:,:,:,jpdsi) 
@@ -488,4 +522,5 @@
       ENDIF
       !
+      ! Compute the budget of Iron
       IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          zwork(:,:,:) =   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe) + trn(:,:,:,jpdfe)   &
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zlim.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 12537)
@@ -2,5 +2,6 @@
    !!======================================================================
    !!                         ***  MODULE p5zlim  ***
-   !! TOP :   PISCES with variable stoichiometry 
+   !! TOP :   PISCES-QUOTA : Computes the various nutrient limitation terms
+   !!                        of phytoplankton
    !!======================================================================
    !! History :   1.0  !  2004     (O. Aumont) Original code
@@ -14,5 +15,5 @@
    USE oce_trc         ! Shared ocean-passive tracers variables
    USE trc             ! Tracers defined
-   USE p4zlim
+   USE p4zlim          ! Nutrient limitation 
    USE sms_pisces      ! PISCES variables
    USE iom             !  I/O manager
@@ -21,7 +22,7 @@
    PRIVATE
 
-   PUBLIC p5z_lim    
-   PUBLIC p5z_lim_init    
-   PUBLIC p5z_lim_alloc
+   PUBLIC p5z_lim           ! called in p4zbio.F90  
+   PUBLIC p5z_lim_init      ! called in trcsms_pisces.F90 
+   PUBLIC p5z_lim_alloc     ! called in trcini_pisces.F90
 
    !! * Shared module variables
@@ -72,23 +73,23 @@
 
    !!* Phytoplankton nutrient limitation terms
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicono3   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpiconh4   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicopo4   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanodop   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicodop   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatdop   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanofer   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicofer   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatfer   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimpic    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimpics   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimphys   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdias   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimpfe    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvnuptk
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvpuptk
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvduptk
-
-   ! Coefficient for iron limitation
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicono3   !: Limitation of NO3 uptake by picophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpiconh4   !: Limitation of NH4 uptake by picophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicopo4   !: Limitation of PO4 uptake by picophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanodop   !: Limitation of DOP uptake by nanophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicodop   !: Limitation of DOP uptake by picophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatdop   !: Limitation of DOP uptake by diatoms
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanofer   !: Limitation of Fe uptake by nanophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicofer   !: Limitation of Fe uptake by picophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatfer   !: Limitation of Fe uptake by diatoms
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimpic    !: Limitation of picophyto PP by nutrients
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimpics   !: Limitation of picophyto PP by nutrients
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimphys   !: Limitation of nanophyto PP by nutrients
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdias   !: Limitation of diatoms PP by nutrients
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimpfe    !: Limitation of picophyto PP by Fe
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvnuptk    !: Maximum potential uptake rate of nanophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvpuptk    !: Maximum potential uptake rate of picophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvduptk    !: Maximum potential uptake rate of diatoms
+
+   ! Coefficient for iron limitation following Flynn and Hipkin (1999)
    REAL(wp) ::  xcoef1   = 0.00167  / 55.85
    REAL(wp) ::  xcoef2   = 1.21E-5 * 14. / 55.85 / 7.625 * 0.5 * 1.5
@@ -110,6 +111,6 @@
       !!                approach. The quota model is derived from theoretical
       !!                models proposed by Pahlow and Oschlies (2009) and 
-      !!                Flynn (2001). Various adaptations from several publications
-      !!                by these authors have been also adopted. 
+      !!                Flynn (2001). Various adaptations from several 
+      !!                publications by these authors have been also adopted. 
       !!
       !! ** Method  : Quota based approach. The quota model is derived from 
@@ -145,6 +146,7 @@
             DO ji = 1, jpi
                ! 
-               ! Tuning of the iron concentration to a minimum level that is set to the detection limit
-               !-------------------------------------
+               ! Tuning of the iron concentration to a minimum level that
+               ! is set to the detection limit
+               ! --------------------------------------------------------
                zno3    = trb(ji,jj,jk,jpno3) / 40.e-6
                zferlim = MAX( 3e-11 * zno3 * zno3, 5e-12 )
@@ -152,6 +154,6 @@
                trb(ji,jj,jk,jpfer) = MAX( trb(ji,jj,jk,jpfer), zferlim )
 
-               ! Computation of the mean relative size of each community
-               ! -------------------------------------------------------
+               ! Computation of the Chl/C ratio of each phytoplankton group
+               ! ----------------------------------------------------------
                z1_trnphy   = 1. / ( trb(ji,jj,jk,jpphy) + rtrn )
                z1_trnpic   = 1. / ( trb(ji,jj,jk,jppic) + rtrn )
@@ -161,7 +163,10 @@
                zdiatchl = trb(ji,jj,jk,jpdch) * z1_trndia
 
-               ! Computation of a variable Ks for iron on diatoms taking into account
-               ! that increasing biomass is made of generally bigger cells
-               !------------------------------------------------
+               ! Computation of a variable Ks for the different phytoplankton
+               ! group as a function of their relative size. Allometry
+               ! from Edwards et al. (2012)
+               ! ------------------------------------------------------------
+
+               ! diatoms
                zsized            = sized(ji,jj,jk)**0.81
                zconcdfe          = concdfer * zsized
@@ -170,4 +175,5 @@
                zconc0dpo4        = concdpo4 * zsized
 
+               ! picophytoplankton
                zsizep            = sizep(ji,jj,jk)**0.81
                zconcpfe          = concpfer * zsizep
@@ -176,4 +182,5 @@
                zconc0ppo4        = concppo4 * zsizep
 
+               ! nanophytoplankton
                zsizen            = sizen(ji,jj,jk)**0.81
                zconcnfe          = concnfer * zsizen
@@ -193,6 +200,9 @@
 
                ! Computation of the optimal allocation parameters
-               ! Based on the different papers by Pahlow et al., and Smith et al.
-               ! -----------------------------------------------------------------
+               ! Based on the different papers by Pahlow et al., and 
+               ! Smith et al.
+               ! ---------------------------------------------------
+
+               ! Nanophytoplankton
                znutlim = MAX( trb(ji,jj,jk,jpnh4) / zconc0nnh4,    &
                  &         trb(ji,jj,jk,jpno3) / zconc0n)
@@ -202,4 +212,6 @@
                znutlim = biron(ji,jj,jk) / zconcnfe
                fananof = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+
+               ! Picophytoplankton
                znutlim = MAX( trb(ji,jj,jk,jpnh4) / zconc0pnh4,    &
                  &         trb(ji,jj,jk,jpno3) / zconc0p)
@@ -209,4 +221,6 @@
                znutlim = biron(ji,jj,jk) / zconcpfe
                fapicof = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+
+               ! Diatoms
                znutlim = MAX( trb(ji,jj,jk,jpnh4) / zconc1dnh4,    &
                  &         trb(ji,jj,jk,jpno3) / zconc1d )
@@ -216,7 +230,8 @@
                znutlim = biron(ji,jj,jk) / zconcdfe
                fadiatf = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
-               !
-               ! Michaelis-Menten Limitation term for nutrients Small bacteria
-               ! -------------------------------------------------------------
+
+               ! Michaelis-Menten Limitation term by nutrients of
+               ! heterotrophic bacteria
+               ! -------------------------------------------------
                zbactnh4 = trb(ji,jj,jk,jpnh4) / ( concbnh4 + trb(ji,jj,jk,jpnh4) )
                zbactno3 = trb(ji,jj,jk,jpno3) / ( concbno3 + trb(ji,jj,jk,jpno3) ) * (1. - zbactnh4)
@@ -226,9 +241,14 @@
                zlim3    = biron(ji,jj,jk) / ( concbfe + biron(ji,jj,jk) )
                zlim4    = trb(ji,jj,jk,jpdoc) / ( xkdoc   + trb(ji,jj,jk,jpdoc) )
+               ! Xlimbac is used for DOC solubilization whereas xlimbacl
+               ! is used for all the other bacterial-dependent terms
+               ! -------------------------------------------------------
                xlimbacl(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 )
                xlimbac (ji,jj,jk) = xlimbacl(ji,jj,jk) * zlim4
-               !
-               ! Michaelis-Menten Limitation term for nutrients Small flagellates
-               ! -----------------------------------------------
+               
+               ! Michaelis-Menten Limitation term by nutrients: Nanophyto
+               ! --------------------------------------------------------
+               !
+               ! Limitation of N based nutrients uptake (NO3 and NH4)
                zfalim = (1.-fanano) / fanano
                xnanonh4(ji,jj,jk) = (1. - fanano) * trb(ji,jj,jk,jpnh4) / ( zfalim * zconc0nnh4 + trb(ji,jj,jk,jpnh4) )
@@ -236,4 +256,5 @@
                &                    * (1. - xnanonh4(ji,jj,jk))
                !
+               ! Limitation of P based nutrients (PO4 and DOP)
                zfalim = (1.-fananop) / fananop
                xnanopo4(ji,jj,jk) = (1. - fananop) * trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zfalim * zconc0npo4 )
@@ -242,7 +263,11 @@
                xnanodop(ji,jj,jk) = 0.
                !
+               ! Limitation of Fe uptake
                zfalim = (1.-fananof) / fananof
                xnanofer(ji,jj,jk) = (1. - fananof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcnfe )
                !
+               ! The minimum iron quota depends on the size of PSU, respiration
+               ! and the reduction of nitrate following the parameterization 
+               ! proposed by Flynn and Hipkin (1999)
                zratiof   = trb(ji,jj,jk,jpnfe) * z1_trnphy
                zqfemn = xcoef1 * znanochl + xcoef2 + xcoef3 * xnanono3(ji,jj,jk)
@@ -257,16 +282,20 @@
                &          / (xqnnmax(ji,jj,jk) - xqnnmin(ji,jj,jk) ) ) * xqnnmax(ji,jj,jk)  &
                &          / (zration + rtrn)
-               !  The value of the optimal quota in the formulation below
-               !  has been found by solving a non linear equation
+               ! The value of the optimal quota in the formulation below
+               ! has been found by solving a non linear equation
                zlim1f = max(0., ( 1.086 - xqnnmin(ji,jj,jk) )  &
                &          / (xqnnmax(ji,jj,jk) - xqnnmin(ji,jj,jk) ) ) * xqnnmax(ji,jj,jk)
                zlim3  = MAX( 0.,( zratiof - zqfemn ) / qfnopt )
+               ! computation of the various limitation terms of nanophyto
+               ! growth and PP
                xlimnfe (ji,jj,jk) = MIN( 1., zlim3 )
                xlimphy (ji,jj,jk) = MIN( 1., zlim1, zlim3 )
                xlimphys(ji,jj,jk) = MIN( 1., zlim1/( zlim1f + rtrn ), zlim3 )
                xlimnpn (ji,jj,jk) = MIN( 1., zlim1)
-               !
-               ! Michaelis-Menten Limitation term for nutrients picophytoplankton
-               ! ----------------------------------------------------------------
+
+               ! Michaelis-Menten Limitation term by nutrients: Picophyto
+               ! --------------------------------------------------------
+               !
+               ! Limitation of N based nutrients uptake (NO3 and NH4) 
                zfalim = (1.-fapico) / fapico 
                xpiconh4(ji,jj,jk) = (1. - fapico) * trb(ji,jj,jk,jpnh4) / ( zfalim * zconc0pnh4 + trb(ji,jj,jk,jpnh4) )
@@ -274,4 +303,5 @@
                &                    * (1. - xpiconh4(ji,jj,jk))
                !
+               ! Limitation of P based nutrients uptake (PO4 and DOP)
                zfalim = (1.-fapicop) / fapicop 
                xpicopo4(ji,jj,jk) = (1. - fapicop) * trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zfalim * zconc0ppo4 )
@@ -283,4 +313,7 @@
                xpicofer(ji,jj,jk) = (1. - fapicof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcpfe )
                !
+               ! The minimum iron quota depends on the size of PSU, respiration
+               ! and the reduction of nitrate following the parameterization 
+               ! proposed by Flynn and Hipkin (1999)
                zratiof   = trb(ji,jj,jk,jppfe) * z1_trnpic
                zqfemp = xcoef1 * zpicochl + xcoef2 + xcoef3 * xpicono3(ji,jj,jk)
@@ -295,16 +328,21 @@
                &          / (xqnpmax(ji,jj,jk) - xqnpmin(ji,jj,jk) ) ) * xqnpmax(ji,jj,jk)  &
                &          / (zration + rtrn)
-               !  The value of the optimal quota in the formulation below
-               !  has been found by solving a non linear equation
+               ! The value of the optimal quota in the formulation below
+               ! has been found by solving a non linear equation
                zlim1f   = max(0., (1.367 - xqnpmin(ji,jj,jk) )  &
                &          / (xqnpmax(ji,jj,jk) - xqnpmin(ji,jj,jk) ) ) * xqnpmax(ji,jj,jk)
                zlim3    = MAX( 0.,( zratiof - zqfemp ) / qfpopt )
+
+               ! computation of the various limitation terms of picophyto
+               ! growth and PP
                xlimpfe (ji,jj,jk) = MIN( 1., zlim3 )
                xlimpic (ji,jj,jk) = MIN( 1., zlim1, zlim3 )
                xlimnpp (ji,jj,jk) = MIN( 1., zlim1 )
                xlimpics(ji,jj,jk) = MIN( 1., zlim1/( zlim1f + rtrn ), zlim3 )
-               !
-               !   Michaelis-Menten Limitation term for nutrients Diatoms
-               !   ------------------------------------------------------
+
+               ! Michaelis-Menten Limitation term by nutrients : Diatoms
+               ! -------------------------------------------------------
+               !
+               ! Limitation of N based nutrients uptake (NO3 and NH4)
                zfalim = (1.-fadiat) / fadiat 
                xdiatnh4(ji,jj,jk) = (1. - fadiat) * trb(ji,jj,jk,jpnh4) / ( zfalim * zconc1dnh4 + trb(ji,jj,jk,jpnh4) )
@@ -312,4 +350,5 @@
                &                    * (1. - xdiatnh4(ji,jj,jk))
                !
+               ! Limitation of P based nutrients uptake (PO4 and DOP)
                zfalim = (1.-fadiatp) / fadiatp
                xdiatpo4(ji,jj,jk) = (1. - fadiatp) * trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zfalim * zconc0dpo4 )
@@ -318,7 +357,11 @@
                xdiatdop(ji,jj,jk) = 0.
                !
+               ! Limitation of Fe uptake
                zfalim = (1.-fadiatf) / fadiatf
                xdiatfer(ji,jj,jk) = (1. - fadiatf) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcdfe )
                !
+               ! The minimum iron quota depends on the size of PSU, respiration
+               ! and the reduction of nitrate following the parameterization 
+               ! proposed by Flynn and Hipkin (1999)
                zratiof   = trb(ji,jj,jk,jpdfe) * z1_trndia
                zqfemd = xcoef1 * zdiatchl + xcoef2 + xcoef3 * xdiatno3(ji,jj,jk)
@@ -333,6 +376,6 @@
                &          / (xqndmax(ji,jj,jk) - xqndmin(ji,jj,jk) ) )   &
                &          * xqndmax(ji,jj,jk) / (zration + rtrn)
-               !  The value of the optimal quota in the formulation below
-               !  has been found by solving a non linear equation
+               ! The value of the optimal quota in the formulation below
+               ! has been found by solving a non linear equation
                zlim1f   = max(0., (1.077 - xqndmin(ji,jj,jk) )    &
                &          / (xqndmax(ji,jj,jk) - xqndmin(ji,jj,jk) ) )   &
@@ -340,4 +383,6 @@
                zlim3    = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi(ji,jj) )
                zlim4    = MAX( 0., ( zratiof - zqfemd ) / qfdopt )
+               ! computation of the various limitation terms of diatoms
+               ! growth and PP
                xlimdfe(ji,jj,jk) = MIN( 1., zlim4 )
                xlimdia(ji,jj,jk) = MIN( 1., zlim1, zlim3, zlim4 )
@@ -357,6 +402,7 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ! Size estimation of nanophytoplankton
-               ! ------------------------------------
+               ! Size estimation of nanophytoplankton based on total biomass
+               ! Assumes that larger biomass implies addition of larger cells
+               ! ------------------------------------------------------------
                zcoef = trb(ji,jj,jk,jpphy) - MIN(xsizephy, trb(ji,jj,jk,jpphy) )
                sizena(ji,jj,jk) = 1. + ( xsizern -1.0 ) * zcoef / ( xsizephy + zcoef )
@@ -373,6 +419,7 @@
 
 
-               ! Size estimation of picophytoplankton
-               ! ------------------------------------
+               ! Size estimation of picophytoplankton based on total biomass
+               ! Assumes that larger biomass implies addition of larger cells
+               ! ------------------------------------------------------------
                zcoef = trb(ji,jj,jk,jppic) - MIN(xsizepic, trb(ji,jj,jk,jppic) )
                sizepa(ji,jj,jk) = 1. + ( xsizerp -1.0 ) * zcoef / ( xsizepic + zcoef )
@@ -389,6 +436,7 @@
                xqppmax(ji,jj,jk) = MIN( qppmax, xqppmax(ji,jj,jk) )
 
-               ! Size estimation of diatoms
-               ! --------------------------
+               ! Size estimation of diatoms based on total biomass
+               ! Assumes that larger biomass implies addition of larger cells
+               ! ------------------------------------------------------------
                zcoef = trb(ji,jj,jk,jpdia) - MIN(xsizedia, trb(ji,jj,jk,jpdia) )
                sized(ji,jj,jk) = 1. + ( xsizerd - 1.0 ) * zcoef / ( xsizedia + zcoef )
@@ -412,4 +460,6 @@
 
       ! Compute the fraction of nanophytoplankton that is made of calcifiers
+      ! This is a purely adhoc formulation described in Aumont et al. (2015)
+      ! This fraction depends on nutrient limitation, light, temperature
       ! --------------------------------------------------------------------
       DO jk = 1, jpkm1
@@ -441,4 +491,9 @@
                   &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  )
                nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
+               !
+               ! redox factor computed from NO3 levels
+               nitrfac2(ji,jj,jk) = MAX( 0.e0,       ( 1.E-6 - trb(ji,jj,jk,jpno3) )  &
+                  &                                / ( 1.E-6 + trb(ji,jj,jk,jpno3) ) )
+               nitrfac2(ji,jj,jk) = MIN( 1., nitrfac2(ji,jj,jk) )
             END DO
          END DO
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmeso.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 12537)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                         ***  MODULE p5zmeso  ***
-   !! TOP :   PISCES Compute the sources/sinks for mesozooplankton
+   !! TOP :   PISCES-QUOTA Compute the sources/sinks for mesozooplankton
    !!======================================================================
    !! History :   1.0  !  2002     (O. Aumont) Original code
@@ -9,6 +9,7 @@
    !!             3.6  !  2015-05  (O. Aumont) PISCES quota
    !!----------------------------------------------------------------------
-   !!   p5z_meso       :   Compute the sources/sinks for mesozooplankton
-   !!   p5z_meso_init  :   Initialization of the parameters for mesozooplankton
+   !!   p5z_meso       : Compute the sources/sinks for mesozooplankton
+   !!   p5z_meso_init  : Initialization of the parameters for mesozooplankton
+   !!   p5z_meso_alloc : Allocate variables for mesozooplankton 
    !!----------------------------------------------------------------------
    USE oce_trc         !  shared variables between ocean and passive tracers
@@ -23,5 +24,5 @@
    PUBLIC   p5z_meso              ! called in p5zbio.F90
    PUBLIC   p5z_meso_init         ! called in trcsms_pisces.F90
-   PUBLIC   p5z_meso_alloc
+   PUBLIC   p5z_meso_alloc        ! called in trcini_pisces.F90
 
    !! * Shared module variables
@@ -53,7 +54,6 @@
    LOGICAL,  PUBLIC ::  bmetexc2     !: Use of excess carbon for respiration
    LOGICAL , PUBLIC ::  ln_dvm_meso  !: Boolean to activate DVM of mesozooplankton
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: depmig
-   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: kmig
-
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: depmig  !: DVM of mesozooplankton : migration depth
+   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: kmig    !: Vertical indice of the the migration depth
 
    !!----------------------------------------------------------------------
@@ -70,8 +70,14 @@
       !!
       !! ** Purpose :   Compute the sources/sinks for mesozooplankton
-      !!
-      !! ** Method  : - ???
+      !!                This includes ingestion and assimilation, flux feeding
+      !!                and mortality. We use an active prey switching  
+      !!                parameterization Morozov and Petrovskii (2013). 
+      !!                All living compartments and mesozooplankton
+      !!                are potential preys of mesozooplankton as well as small
+      !!                sinking particles 
+      !!
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, knt ! ocean time step
+      !
       INTEGER  :: ji, jj, jk, jkt
       REAL(wp) :: zcompadi, zcompaph, zcompapoc, zcompaz, zcompam, zcompames
@@ -115,7 +121,9 @@
       !
       ! Diurnal vertical migration of mesozooplankton
+      ! Computation of the migration depth
       ! ---------------------------------------------
       IF (ln_dvm_meso) CALL p5z_meso_depmig
 
+      ! Use of excess carbon for metabolism
       zmetexcess = 0.0
       IF ( bmetexc2 ) zmetexcess = 1.0
@@ -127,12 +135,15 @@
                zfact     = xstep * tgfunc2(ji,jj,jk) * zcompam
 
-               !   Michaelis-Menten mortality rates of mesozooplankton
-               !   ---------------------------------------------------
+               !  linear mortality of mesozooplankton
+               !  A michaelis menten modulation term is used to avoid extinction of 
+               !  mesozooplankton at very low food concentrations
+               !  -----------------------------------------------------------------
                zrespz   = resrat2 * zfact * ( trb(ji,jj,jk,jpmes) / ( xkmort + trb(ji,jj,jk,jpmes) )  &
                &          + 3. * nitrfac(ji,jj,jk) )
 
-               !   Zooplankton mortality. A square function has been selected with
-               !   no real reason except that it seems to be more stable and may mimic predation
-               !   ---------------------------------------------------------------
+               !  Zooplankton quadratic mortality. A square function has been selected with
+               !  to mimic predation and disease (density dependent mortality). It also tends
+               !  to stabilise the model
+               !  -------------------------------------------------------------------------
                ztortz   = mzrat2 * 1.e6 * zfact * trb(ji,jj,jk,jpmes) * (1. - nitrfac(ji,jj,jk))
 
@@ -147,5 +158,10 @@
 
                !   Mesozooplankton grazing
-               !   ------------------------
+               ! The total amount of food is the sum of all preys accessible to mesozooplankton 
+               ! multiplied by their food preference
+               ! A threshold can be specified in the namelist (xthresh2). However, when food 
+               ! concentration is close to this threshold, it is decreased to avoid the 
+               ! accumulation of food in the mesozoopelagic domain
+               ! -------------------------------------------------------------------------------
                zfood     = xpref2d * zcompadi + xpref2z * zcompaz + xpref2n * zcompaph + xpref2c * zcompapoc   &
                &           + xpref2m * zcompames 
@@ -154,14 +170,21 @@
                zgraze2   = grazrat2 * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpmes) * (1. - nitrfac(ji,jj,jk)) 
 
-               !   An active switching parameterization is used here.
-               !   We don't use the KTW parameterization proposed by 
-               !   Vallina et al. because it tends to produce too steady biomass
-               !   composition and the variance of Chl is too low as it grazes
-               !   too strongly on winning organisms. We use a generalized
-               !   switching parameterization proposed by Morozov and 
-               !   Petrovskii (2013)
-               !   ------------------------------------------------------------  
+               ! An active switching parameterization is used here.
+               ! We don't use the KTW parameterization proposed by 
+               ! Vallina et al. because it tends to produce too steady biomass
+               ! composition and the variance of Chl is too low as it grazes
+               ! too strongly on winning organisms. We use a generalized
+               ! switching parameterization proposed by Morozov and 
+               ! Petrovskii (2013)
+               ! ------------------------------------------------------------  
+               ! The width of the selection window is increased when preys
+               ! have low abundance, .i.e. zooplankton become less specific 
+               ! to avoid starvation.
+               ! ----------------------------------------------------------
                zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
                zsigma = 0.5 + 1.0 * zsigma
+               ! Nanophytoplankton and diatoms are the only preys considered
+               ! to be close enough to have potential interference
+               ! -----------------------------------------------------------
                zdiffdn = exp( -ABS(log(3.0 * sizen(ji,jj,jk) / (5.0 * sized(ji,jj,jk) + rtrn )) )**2 / zsigma**2 )
                ztmp1 = xpref2n * zcompaph * ( zcompaph + zdiffdn * zcompadi )
@@ -194,6 +217,7 @@
                zgrazpof  = zgrazpoc * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn)
 
-               !   Mesozooplankton flux feeding on GOC
-               !   ----------------------------------
+               !  Mesozooplankton flux feeding on GOC and POC. The feeding pressure
+               ! is proportional to the flux
+               !  ------------------------------------------------------------------
                zgrazffeg = grazflux  * xstep * wsbio4(ji,jj,jk)      &
                &           * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpgoc) * trb(ji,jj,jk,jpmes)  &
@@ -211,6 +235,6 @@
                zgraztotc  = zgrazdc + zgrazz + zgraznc + zgrazm + zgrazpoc + zgrazffep + zgrazffeg
 
-               !   Compute the proportion of filter feeders
-               !   ----------------------------------------  
+               ! Compute the proportion of filter feeders. It is assumed steady state.
+               ! ---------------------------------------------------------------------  
                zproport  = (zgrazffep + zgrazffeg)/(rtrn + zgraztotc)
 
@@ -228,4 +252,5 @@
                zfracp    = zfracc * trb(ji,jj,jk,jpgop) / (trb(ji,jj,jk,jpgoc) + rtrn)
 
+               ! Flux feeding is multiplied by the fractional biomass of flux feeders
                zgrazffep = zproport * zgrazffep   ;   zgrazffeg = zproport * zgrazffeg
                zgrazfffp = zproport * zgrazfffp   ;   zgrazfffg = zproport * zgrazfffg
@@ -241,5 +266,4 @@
                &            + zgrazffpp + zgrazffpg
 
-
                ! Total grazing ( grazing by microzoo is already computed in p5zmicro )
                zgrazing(ji,jj,jk) = zgraztotc
@@ -251,7 +275,12 @@
                zgrasratp  =  (zgraztotp + rtrn) / ( zgraztotc + rtrn )
 
-               !   Growth efficiency is made a function of the quality 
-               !   and the quantity of the preys
-               !   ---------------------------------------------------
+               ! Mesozooplankton efficiency. 
+               ! We adopt a formulation proposed by Mitra et al. (2007)
+               ! The gross growth efficiency is controled by the most limiting nutrient.
+               ! Growth is also further decreased when the food quality is poor. This is currently
+               ! hard coded : it can be decreased by up to 50% (zepsherq)
+               ! GGE can also be decreased when food quantity is high, zepsherf (Montagnes and 
+               ! Fulton, 2012)
+               ! -----------------------------------------------------------------------------------
                zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / ferat3)
                zbeta     = MAX(0., (epsher2 - epsher2min) )
@@ -261,5 +290,6 @@
                !   Respiration of mesozooplankton
                !   Excess carbon in the food is used preferentially
-               !   ----------------  ------------------------------
+               !   when bmetexc2 is set to .true.
+               !   -----------------------------------------------
                zexcess  = zgraztotc * zepsherf * (1.0 - zepshert) * zmetexcess 
                zbasresb = MAX(0., zrespz - zexcess)
@@ -344,6 +374,8 @@
                zfracal = trb(ji,jj,jk,jpcal) / ( trb(ji,jj,jk,jpgoc) + rtrn )
                zgrazcal = zgrazffeg * (1. - part2) * zfracal
-               !  calcite production
-               !  ------------------
+               ! Calcite production
+               ! Calcite remineralization due to zooplankton activity
+               ! part2 of the ingested calcite is dissolving in the acidic gut
+               ! -------------------------------------------------------------
                zprcaca = xfracal(ji,jj,jk) * zgraznc
                prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
@@ -356,9 +388,12 @@
       END DO
 
+      ! Computation of the effect of DVM by mesozooplankton
+      ! This part is only activated if ln_dvm_meso is set to true
+      ! The parameterization has been published in Gorgues et al. (2019).
+      ! -----------------------------------------------------------------
       IF (ln_dvm_meso) THEN
          ALLOCATE( zgramigrem(jpi,jpj), zgramigref(jpi,jpj), zgramigpoc(jpi,jpj), zgramigpof(jpi,jpj) )
          ALLOCATE( zgramigrep(jpi,jpj), zgramigren(jpi,jpj), zgramigpop(jpi,jpj), zgramigpon(jpi,jpj) )
          ALLOCATE( zgramigdoc(jpi,jpj), zgramigdon(jpi,jpj), zgramigdop(jpi,jpj) )
-
          ALLOCATE( zstrn(jpi,jpj) )
          zgramigrem(:,:)  = 0.0   ;   zgramigref(:,:) = 0.0
@@ -385,9 +420,11 @@
 
 
+        ! Compute the amount of materials that will go into vertical migration
+        ! This fraction is sumed over the euphotic zone and is removed from 
+        ! the fluxes driven by mesozooplankton in the euphotic zone.
+        ! --------------------------------------------------------------------
          DO jk = 1, jpk
             DO jj = 1, jpj
                DO ji = 1, jpi
-
-                  !   Compute the amount of materials that will go into vertical migration
                   zmigreltime = (1. - zstrn(ji,jj))
                   IF ( gdept_n(ji,jj,jk) <= heup(ji,jj) ) THEN
@@ -431,4 +468,7 @@
          END DO
 
+         ! The inorganic and organic fluxes induced by migrating organisms are added at the 
+         ! the migration depth (corresponding indice is set by kmig)
+         ! --------------------------------------------------------------------------------
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -456,11 +496,14 @@
          DEALLOCATE( zgramigdoc, zgramigdon, zgramigdop )
          DEALLOCATE( zstrn )
-
+      ! End of the ln_dvm_meso part
       ENDIF
 
+      !   Update the arrays TRA which contain the biological sources and sinks
+      !   This only concerns the variables which are affected by DVM (inorganic 
+      !   nutrients, DOC agands, and particulate organic carbon). 
+      !   ---------------------------------------------------------------------
       DO jk = 1, jpk
          DO jj = 1, jpj
             DO ji = 1, jpi
-               !   Update the arrays TRA which contain the biological sources and sinks
                tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zgrarep(ji,jj,jk) 
                tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zgraren(ji,jj,jk)
@@ -526,11 +569,11 @@
       !! ** Purpose :   Initialization of mesozooplankton parameters
       !!
-      !! ** Method  :   Read the nampismes namelist and check the parameters
+      !! ** Method  :   Read the namp5zmes namelist and check the parameters
       !!      called at the first timestep (nittrc000)
       !!
-      !! ** input   :   Namelist nampismes
+      !! ** input   :   Namelist namp5zmes
       !!
       !!----------------------------------------------------------------------
-      INTEGER :: ios                 ! Local integer output status for namelist read
+      INTEGER :: ios    ! Local integer output status for namelist read
       !!
       NAMELIST/namp5zmes/part2, bmetexc2, grazrat2, resrat2, mzrat2, xpref2c, xpref2n, xpref2z, &
@@ -540,11 +583,11 @@
       !!----------------------------------------------------------------------
       !
-      REWIND( numnatp_ref )              ! Namelist nampismes in reference namelist : Pisces mesozooplankton
+      REWIND( numnatp_ref )              ! Namelist namp5zmes in reference namelist : Pisces mesozooplankton
       READ  ( numnatp_ref, namp5zmes, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampismes in reference namelist' )
-      !
-      REWIND( numnatp_cfg )              ! Namelist nampismes in configuration namelist : Pisces mesozooplankton
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp5zmes in reference namelist' )
+      !
+      REWIND( numnatp_cfg )              ! Namelist namp5zmes in configuration namelist : Pisces mesozooplankton
       READ  ( numnatp_cfg, namp5zmes, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nampismes in configuration namelist' )
+902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp5zmes in configuration namelist' )
       IF(lwm) WRITE ( numonp, namp5zmes )
       !
@@ -628,4 +671,7 @@
       END DO
 
+      ! Compute the difference between surface values and the mean values in the mesopelagic
+      ! domain
+      ! ------------------------------------------------------------------------------------
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -634,7 +680,8 @@
          END DO
       END DO
-      !
+
       ! Computation of the migration depth based on the parameterization of 
-      !
+      ! Bianchi et al. (2013)
+      ! -------------------------------------------------------------------
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -648,5 +695,4 @@
       ! Computation of the corresponding jk indice 
       ! ------------------------------------------
-      ! 
       DO jk = 1, jpk-1
          DO jj = 1, jpj
@@ -665,5 +711,4 @@
       ! that it falls above the OMZ
       ! -----------------------------------------------------------------------
-      !
       DO ji =1, jpi
          DO jj = 1, jpj
@@ -685,5 +730,5 @@
    INTEGER FUNCTION p5z_meso_alloc()
       !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE p4z_meso_alloc  ***
+      !!                     ***  ROUTINE p5z_meso_alloc  ***
       !!----------------------------------------------------------------------
       !
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmicro.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmicro.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmicro.F90	(revision 12537)
@@ -97,4 +97,5 @@
       ENDIF
       !
+      ! Use of excess carbon for metabolism
       zmetexcess = 0.0
       IF ( bmetexc ) zmetexcess = 1.0
@@ -109,17 +110,23 @@
                zproport  = min(1.0, exp(-1.1 * MAX(0., ( sized(ji,jj,jk) - 1.8 ))**0.8 ))
                zproport2 = sizen(ji,jj,jk)**(-0.54)
-               !   Michaelis-Menten mortality rates of microzooplankton
-               !   -----------------------------------------------------
+               !  linear mortality of mesozooplankton
+               !  A michaelis menten modulation term is used to avoid extinction of 
+               !  microzooplankton at very low food concentrations. Mortality is 
+               !  enhanced in low O2 waters
+               !  -----------------------------------------------------------------
                zrespz = resrat * zfact * ( trb(ji,jj,jk,jpzoo) / ( xkmort + trb(ji,jj,jk,jpzoo) )  &
                &        + 3. * nitrfac(ji,jj,jk) )
 
-               !   Zooplankton mortality. A square function has been selected with
-               !   no real reason except that it seems to be more stable and may mimic predation.
-               !   ------------------------------------------------------------------------------
+               !  Zooplankton quadratic mortality. A square function has been selected with
+               !  to mimic predation and disease (density dependent mortality). It also tends
+               !  to stabilise the model
+               !  -------------------------------------------------------------------------
                ztortz = mzrat * 1.e6 * zfact * trb(ji,jj,jk,jpzoo) * (1. - nitrfac(ji,jj,jk))
 
                !   Computation of the abundance of the preys
                !   A threshold can be specified in the namelist
-               !   --------------------------------------------
+               !   Nanophyto and diatoms have a specific treatment with 
+               !   teir preference decreasing with size.
+               !   --------------------------------------------------------
                zcompadi  = zproport * MAX( ( trb(ji,jj,jk,jpdia) - xthreshdia ), 0.e0 )
                zcompaph  = zproport2 * MAX( ( trb(ji,jj,jk,jpphy) - xthreshphy ), 0.e0 )
@@ -128,6 +135,11 @@
                zcompapoc = MAX( ( trb(ji,jj,jk,jppoc) - xthreshpoc ), 0.e0 )
                
-               !   Microzooplankton grazing
-               !   ------------------------
+               ! Microzooplankton grazing
+               ! The total amount of food is the sum of all preys accessible to mesozooplankton 
+               ! multiplied by their food preference
+               ! A threshold can be specified in the namelist (xthresh). However, when food 
+               ! concentration is close to this threshold, it is decreased to avoid the 
+               ! accumulation of food in the mesozoopelagic domain
+               ! -------------------------------------------------------------------------------
                zfood     = xprefn * zcompaph + xprefc * zcompapoc + xprefd * zcompadi   &
                &           + xprefz * zcompaz + xprefp * zcompapi
@@ -136,12 +148,16 @@
                zgraze    = grazrat * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpzoo) * (1. - nitrfac(ji,jj,jk)) 
 
-               !   An active switching parameterization is used here.
-               !   We don't use the KTW parameterization proposed by 
-               !   Vallina et al. because it tends to produce to steady biomass
-               !   composition and the variance of Chl is too low as it grazes
-               !   too strongly on winning organisms. Thus, instead of a square
-               !   a 1.5 power value is used which decreases the pressure on the
-               !   most abundant species
-               !   ------------------------------------------------------------  
+               ! An active switching parameterization is used here.
+               ! We don't use the KTW parameterization proposed by 
+               ! Vallina et al. because it tends to produce too steady biomass
+               ! composition and the variance of Chl is too low as it grazes
+               ! too strongly on winning organisms. We use a generalized
+               ! switching parameterization proposed by Morozov and 
+               ! Petrovskii (2013)
+               ! ------------------------------------------------------------  
+               ! The width of the selection window is increased when preys
+               ! have low abundance, .i.e. zooplankton become less specific 
+               ! to avoid starvation.
+               ! ----------------------------------------------------------
                zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
                zsigma = 0.5 + 1.0*zsigma
@@ -189,13 +205,18 @@
                zgrazing(ji,jj,jk) = zgraztotc
 
-               !   Stoichiometruc ratios of the food ingested by zooplanton 
-               !   --------------------------------------------------------
+               ! Stoichiometruc ratios of the food ingested by zooplanton 
+               ! --------------------------------------------------------
                zgrasratf =  (zgraztotf + rtrn) / ( zgraztotc + rtrn )
                zgrasratn =  (zgraztotn + rtrn) / ( zgraztotc + rtrn )
                zgrasratp =  (zgraztotp + rtrn) / ( zgraztotc + rtrn )
 
-               !   Growth efficiency is made a function of the quality 
-               !   and the quantity of the preys
-               !   ---------------------------------------------------
+               ! Mesozooplankton efficiency. 
+               ! We adopt a formulation proposed by Mitra et al. (2007)
+               ! The gross growth efficiency is controled by the most limiting nutrient.
+               ! Growth is also further decreased when the food quality is poor. This is currently
+               ! hard coded : it can be decreased by up to 50% (zepsherq)
+               ! GGE can also be decreased when food quantity is high, zepsherf (Montagnes and 
+               ! Fulton, 2012)
+               ! -----------------------------------------------------------------------------------
                zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / ferat3)
                zbeta     = MAX( 0., (epsher - epshermin) )
@@ -204,7 +225,7 @@
                zepsherv  = zepsherf * zepshert * zepsherq
 
-               !   Respiration of microzooplankton
-               !   Excess carbon in the food is used preferentially
-               !   ------------------------------------------------
+               ! Respiration of microzooplankton
+               ! Excess carbon in the food is used preferentially
+               ! ------------------------------------------------
                zexcess  = zgraztotc * zepsherf * (1.0 - zepshert) * zmetexcess
                zbasresb = MAX(0., zrespz - zexcess)
@@ -212,7 +233,7 @@
                zrespirc = srespir * zepsherv * zgraztotc + zbasresb
                
-               !   When excess carbon is used, the other elements in excess
-               !   are also used proportionally to their abundance
-               !   --------------------------------------------------------
+               ! When excess carbon is used, the other elements in excess
+               ! are also used proportionally to their abundance
+               ! --------------------------------------------------------
                zexcess  = ( zgrasratn/ no3rat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
                zbasresn = zbasresi * zexcess * zgrasratn 
@@ -222,6 +243,6 @@
                zbasresf = zbasresi * zexcess * zgrasratf
 
-               !   Voiding of the excessive elements as DOM
-               !   ----------------------------------------
+               ! Voiding of the excessive elements as DOM
+               ! ----------------------------------------
                zgradoct   = (1. - unassc - zepsherv) * zgraztotc - zbasresi  
                zgradont   = (1. - unassn) * zgraztotn - zepsherv * no3rat3 * zgraztotc - zbasresn
@@ -229,8 +250,8 @@
                zgrareft   = (1. - unassc) * zgraztotf - zepsherv * ferat3 * zgraztotc - zbasresf
 
-               !  Since only semilabile DOM is represented in PISCES
-               !  part of DOM is in fact labile and is then released
-               !  as dissolved inorganic compounds (ssigma)
-               !  --------------------------------------------------
+               ! Since only semilabile DOM is represented in PISCES
+               ! part of DOM is in fact labile and is then released
+               ! as dissolved inorganic compounds (ssigma)
+               ! --------------------------------------------------
                zgradoc =  zgradoct * ssigma
                zgradon =  zgradont * ssigma
@@ -241,6 +262,6 @@
                zgraref = zgrareft
 
-               !   Defecation as a result of non assimilated products
-               !   --------------------------------------------------
+               ! Defecation as a result of non assimilated products
+               ! --------------------------------------------------
                zgrapoc   = zgraztotc * unassc
                zgrapon   = zgraztotn * unassn
@@ -248,6 +269,6 @@
                zgrapof   = zgraztotf * unassc
 
-               !  Addition of respiration to the release of inorganic nutrients
-               !  -------------------------------------------------------------
+               ! Addition of respiration to the release of inorganic nutrients
+               ! -------------------------------------------------------------
                zgrarem = zgrarem + zbasresi + zrespirc
                zgraren = zgraren + zbasresn + zrespirc * no3rat3
@@ -255,6 +276,6 @@
                zgraref = zgraref + zbasresf + zrespirc * ferat3
 
-               !   Update of the TRA arrays
-               !   ------------------------
+               ! Update of the TRA arrays
+               ! ------------------------
                tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zgrarep
                tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zgraren
@@ -296,5 +317,8 @@
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * ztortz  + zgrapof - zgrazpof
                !
-               ! calcite production
+               ! Calcite production
+               ! Calcite remineralization due to zooplankton activity
+               ! part of the ingested calcite is dissolving in the acidic gut
+               ! -------------------------------------------------------------
                zprcaca = xfracal(ji,jj,jk) * zgraznc
                prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
@@ -345,8 +369,8 @@
       !! ** Purpose :   Initialization of microzooplankton parameters
       !!
-      !! ** Method  :   Read the nampiszoo namelist and check the parameters
+      !! ** Method  :   Read the namp5zzoo namelist and check the parameters
       !!                called at the first timestep (nittrc000)
       !!
-      !! ** input   :   Namelist nampiszoo
+      !! ** input   :   Namelist namp5zzoo
       !!
       !!----------------------------------------------------------------------
@@ -359,9 +383,9 @@
       !!----------------------------------------------------------------------
       !
-      REWIND( numnatp_ref )              ! Namelist nampiszoo in reference namelist : Pisces microzooplankton
+      REWIND( numnatp_ref )              ! Namelist namp5zzoo in reference namelist : Pisces microzooplankton
       READ  ( numnatp_ref, namp5zzoo, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp5zzoo in reference namelist' )
       !
-      REWIND( numnatp_cfg )              ! Namelist nampiszoo in configuration namelist : Pisces microzooplankton
+      REWIND( numnatp_cfg )              ! Namelist namp5zzoo in configuration namelist : Pisces microzooplankton
       READ  ( numnatp_cfg, namp5zzoo, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp5zzoo in configuration namelist' )
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmort.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmort.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmort.F90	(revision 12537)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                         ***  MODULE p5zmort  ***
-   !! TOP :   PISCES Compute the mortality terms for phytoplankton
+   !! TOP :   PISCES-QUOTA Compute the mortality terms for phytoplankton
    !!======================================================================
    !! History :   1.0  !  2002     (O. Aumont)  Original code
@@ -14,6 +14,6 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p4zlim
-   USE p5zlim          !  Phytoplankton limitation terms
+   USE p4zlim          !  Phytoplankton limitation terms (p4z)
+   USE p5zlim          !  Phytoplankton limitation terms (p5z)
    USE prtctl_trc      !  print control for debugging
 
@@ -21,15 +21,15 @@
    PRIVATE
 
-   PUBLIC   p5z_mort    
-   PUBLIC   p5z_mort_init    
+   PUBLIC   p5z_mort           ! Called from p4zbio.F90 
+   PUBLIC   p5z_mort_init      ! Called from trcini_pisces.F90 
 
    !! * Shared module variables
-   REAL(wp), PUBLIC :: wchln    !:
-   REAL(wp), PUBLIC :: wchlp   !:
-   REAL(wp), PUBLIC :: wchld   !:
-   REAL(wp), PUBLIC :: wchldm  !:
-   REAL(wp), PUBLIC :: mpratn   !:
-   REAL(wp), PUBLIC :: mpratp  !:
-   REAL(wp), PUBLIC :: mpratd  !:
+   REAL(wp), PUBLIC :: wchln   !! Quadratic mortality rate of nanophytoplankton
+   REAL(wp), PUBLIC :: wchlp   !: Quadratic mortality rate of picophytoplankton
+   REAL(wp), PUBLIC :: wchld   !: Quadratic mortality rate of diatoms
+   REAL(wp), PUBLIC :: wchldm  !: Maximum quadratic mortality rate of diatoms
+   REAL(wp), PUBLIC :: mpratn  !: Linear mortality rate of nanophytoplankton
+   REAL(wp), PUBLIC :: mpratp  !: Linear mortality rate of picophytoplankton
+   REAL(wp), PUBLIC :: mpratd  !: Linear mortality rate of diatoms
 
    !!----------------------------------------------------------------------
@@ -45,5 +45,5 @@
       !!                     ***  ROUTINE p5z_mort  ***
       !!
-      !! ** Purpose :   Calls the different subroutine to initialize and compute
+      !! ** Purpose :   Calls the different subroutine to compute
       !!                the different phytoplankton mortality terms
       !!
@@ -82,15 +82,17 @@
             DO ji = 1, jpi
                zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e-9 ), 0.e0 )
-               !   Squared mortality of Phyto similar to a sedimentation term during
-               !   blooms (Doney et al. 1996)
-               !   -----------------------------------------------------------------
+               ! Quadratic mortality of nano due to aggregation during
+               ! blooms (Doney et al. 1996)
+               ! -----------------------------------------------------
                zrespp = wchln * 1.e6 * xstep * xdiss(ji,jj,jk) * zcompaph * trb(ji,jj,jk,jpphy)
 
-               !   Phytoplankton linear mortality
-               !   ------------------------------
-               ztortp = mpratn * xstep  * zcompaph
+               ! Phytoplankton linear mortality
+               ! A michaelis-menten like term is introduced to avoid 
+               ! extinction of nanophyto in highly limited areas
+               ! ----------------------------------------------------
+               ztortp = mpratn * xstep  * zcompaphi * trb(ji,jj,jk,jpphy) /  ( xkmort + trb(ji,jj,jk,jpphy) )
                zmortp = zrespp + ztortp
 
-               !   Update the arrays TRA which contains the biological sources and sinks
+               ! Update the arrays TRA which contains the biological sources and sinks
 
                zfactn  = trb(ji,jj,jk,jpnph)/(trb(ji,jj,jk,jpphy)+rtrn)
@@ -103,6 +105,6 @@
                tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zmortp * zfactch
                tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zmortp * zfactfe
+               ! Production PIC particles due to mortality
                zprcaca = xfracal(ji,jj,jk) * zmortp
-               !
                prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
                !
@@ -151,15 +153,18 @@
             DO ji = 1, jpi
                zcompaph = MAX( ( trb(ji,jj,jk,jppic) - 1e-9 ), 0.e0 )
-               !  Squared mortality of Phyto similar to a sedimentation term during
-               !  blooms (Doney et al. 1996)
-               !  -----------------------------------------------------------------
+
+               ! Quadratic mortality of pico due to aggregation during
+               ! blooms (Doney et al. 1996)
+               ! -----------------------------------------------------
                zrespp = wchlp * 1.e6 * xstep * xdiss(ji,jj,jk) * zcompaph * trb(ji,jj,jk,jppic)
 
-               !     Phytoplankton mortality 
-               ztortp = mpratp * xstep  * zcompaph
+               ! Phytoplankton linear mortality
+               ! A michaelis-menten like term is introduced to avoid 
+               ! extinction of picophyto in highly limited areas
+               ! ----------------------------------------------------
+               ztortp = mpratp * xstep  * zcompaph * trb(ji,jj,jk,jppic) /  ( xkmort + trb(ji,jj,jk,jppic) )
                zmortp = zrespp + ztortp
 
                !   Update the arrays TRA which contains the biological sources and sinks
-
                zfactn = trb(ji,jj,jk,jpnpi)/(trb(ji,jj,jk,jppic)+rtrn)
                zfactp = trb(ji,jj,jk,jpppi)/(trb(ji,jj,jk,jppic)+rtrn)
@@ -225,7 +230,9 @@
                zrespp2 = 1.e6 * xstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
 
-               !  Phytoplankton linear mortality 
-               !  ------------------------------
-               ztortp2 = mpratd * xstep  * zcompadi
+               ! Phytoplankton linear mortality
+               ! A michaelis-menten like term is introduced to avoid 
+               ! extinction of diatoms in highly limited areas
+               !  ---------------------------------------------------
+               ztortp2 = mpratd * xstep  * zcompadi * trb(ji,jj,jk,jpdia) /  ( xkmort + trb(ji,jj,jk,jpdia) )
                zmortp2 = zrespp2 + ztortp2
 
@@ -273,22 +280,22 @@
       !!                  ***  ROUTINE p5z_mort_init  ***
       !!
-      !! ** Purpose :   Initialization of phytoplankton parameters
-      !!
-      !! ** Method  :   Read the nampismort namelist and check the parameters
+      !! ** Purpose :   Initialization of phytoplankton mortality parameters
+      !!
+      !! ** Method  :   Read the namp5zmort namelist and check the parameters
       !!      called at the first timestep
       !!
-      !! ** input   :   Namelist nampismort
+      !! ** input   :   Namelist namp5zmort
       !!
       !!----------------------------------------------------------------------
-      INTEGER :: ios                 ! Local integer output status for namelist read
+      INTEGER :: ios   ! Local integer output status for namelist read
       !!
       NAMELIST/namp5zmort/ wchln, wchlp, wchld, wchldm, mpratn, mpratp, mpratd
       !!----------------------------------------------------------------------
 
-      REWIND( numnatp_ref )              ! Namelist nampismort in reference namelist : Pisces phytoplankton
+      REWIND( numnatp_ref )              ! Namelist namp5zmort in reference namelist : Pisces phytoplankton
       READ  ( numnatp_ref, namp5zmort, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp5zmort in reference namelist' )
 
-      REWIND( numnatp_cfg )              ! Namelist nampismort in configuration namelist : Pisces phytoplankton
+      REWIND( numnatp_cfg )              ! Namelist namp5zmort in configuration namelist : Pisces phytoplankton
       READ  ( numnatp_cfg, namp5zmort, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp5zmort in configuration namelist' )
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zprod.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zprod.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zprod.F90	(revision 12537)
@@ -238,5 +238,11 @@
       END DO
 
-      ! Computation of the various production terms of nanophytoplankton 
+      ! Computation of the various production and uptake terms of nanophytoplankton 
+      ! Interactions between N and P are modeled according to the Chain Model 
+      ! of Pahlow et al. (2009). Iron uptake is modeled following traditional
+      ! Droop kinetics. When the quota is approaching the maximum achievable
+      ! quota, uptake is downregulated according to a sigmoidal function 
+      ! (power 2), as proposed by Flynn (2003)
+      ! ---------------------------------------------------------------------------
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -245,9 +251,15 @@
                   !  production terms for nanophyto.
                   zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trb(ji,jj,jk,jpphy) * rfact2
-                  ! size computation
+
+                  ! Size computation
+                  ! Size is made a function of the limitation of of phytoplankton growth
+                  ! Strongly limited cells are supposed to be smaller. sizena is the 
+                  ! size at time step t+1 and is thus updated at the end of the 
+                  ! current time step
+                  ! --------------------------------------------------------------------
                   zlimfac = xlimphys(ji,jj,jk) * zprchln(ji,jj,jk) / ( zprmaxn(ji,jj,jk) + rtrn )
                   zsizetmp = 1.0 + 1.3 * ( xsizern - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3)
                   sizena(ji,jj,jk) = min(xsizern, max( sizena(ji,jj,jk), zsizetmp ) )
-                  ! Nutrient uptake terms
+                  ! Maximum potential uptake rate
                   zration = trb(ji,jj,jk,jpnph) / ( trb(ji,jj,jk,jpphy) + rtrn )
                   zratiop = trb(ji,jj,jk,jppph) / ( trb(ji,jj,jk,jpphy) + rtrn )
@@ -261,5 +273,5 @@
                   zpronewn(ji,jj,jk) = zpronmax * xnanono3(ji,jj,jk)
                   zproregn(ji,jj,jk) = zpronmax * xnanonh4(ji,jj,jk)
-                  ! Uptake of phosphorus
+                  ! Uptake of phosphorus and DOP
                   zrat = 1.0 - MIN( 1., zratiop / (xqpnmax(ji,jj,jk) + rtrn) )
                   zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
@@ -279,5 +291,11 @@
       END DO
 
-      ! Computation of the various production terms of picophytoplankton 
+      ! Computation of the various production and uptake terms of picophytoplankton 
+      ! Interactions between N and P are modeled according to the Chain Model 
+      ! of Pahlow et al. (2009). Iron uptake is modeled following traditional
+      ! Droop kinetics. When the quota is approaching the maximum achievable
+      ! quota, uptake is downregulated according to a sigmoidal function 
+      ! (power 2), as proposed by Flynn (2003)
+      ! ---------------------------------------------------------------------------
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -287,8 +305,13 @@
                   zprorcap(ji,jj,jk) = zprpic(ji,jj,jk)  * xlimpic(ji,jj,jk) * trb(ji,jj,jk,jppic) * rfact2
                   ! Size computation
+                  ! Size is made a function of the limitation of of phytoplankton growth
+                  ! Strongly limited cells are supposed to be smaller. sizepa is
+                  ! size at time step t+1 and is thus updated at the end of the 
+                  ! current time step
+                  ! --------------------------------------------------------------------
                   zlimfac = zprchlp(ji,jj,jk)  * xlimpics(ji,jj,jk) / ( zprmaxp(ji,jj,jk) + rtrn )
                   zsizetmp = 1.0 + 1.3 * ( xsizerp - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3)
                   sizepa(ji,jj,jk) = min(xsizerp, max( sizepa(ji,jj,jk), zsizetmp ) )
-                  ! Nutrient uptake
+                  ! Maximum potential uptake rate of nutrients
                   zration = trb(ji,jj,jk,jpnpi) / ( trb(ji,jj,jk,jppic) + rtrn )
                   zratiop = trb(ji,jj,jk,jpppi) / ( trb(ji,jj,jk,jppic) + rtrn )
@@ -320,5 +343,11 @@
       END DO
 
-      ! Computation of the various production terms of diatoms
+      ! Computation of the various production and uptake terms of diatoms
+      ! Interactions between N and P are modeled according to the Chain Model 
+      ! of Pahlow et al. (2009). Iron uptake is modeled following traditional
+      ! Droop kinetics. When the quota is approaching the maximum achievable
+      ! quota, uptake is downregulated according to a sigmoidal function 
+      ! (power 2), as proposed by Flynn (2003)
+      ! ---------------------------------------------------------------------------
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -327,9 +356,14 @@
                   !  production terms for diatomees
                   zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trb(ji,jj,jk,jpdia) * rfact2
-                  ! size computation
+                  ! Size computation
+                  ! Size is made a function of the limitation of of phytoplankton growth
+                  ! Strongly limited cells are supposed to be smaller. sizeda is
+                  ! size at time step t+1 and is thus updated at the end of the 
+                  ! current time step. 
+                  ! --------------------------------------------------------------------
                   zlimfac = zprchld(ji,jj,jk) * xlimdias(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn )
                   zsizetmp = 1.0 + 1.3 * ( xsizerd - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3)
                   sizeda(ji,jj,jk) = min(xsizerd, max( sizeda(ji,jj,jk), zsizetmp ) )
-                  ! Nutrient uptake
+                  ! Maximum potential uptake rate of nutrients
                   zration = trb(ji,jj,jk,jpndi) / ( trb(ji,jj,jk,jpdia) + rtrn )
                   zratiop = trb(ji,jj,jk,jppdi) / ( trb(ji,jj,jk,jpdia) + rtrn )
@@ -361,4 +395,7 @@
       END DO
 
+      ! Production of Chlorophyll. The formulation proposed by Geider et al.
+      ! is adopted here.
+      ! --------------------------------------------------------------------
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -375,5 +412,5 @@
                   zprochlp = thetanpm * zprod / ( zpislopeadp(ji,jj,jk) * zpicotot + rtrn )
                   zprochlp = MAX(zprochlp, chlcmin * 12. * zprorcap(ji,jj,jk) )
-                  !  production terms for diatomees ( chlorophyll )
+                  !  production terms for diatoms ( chlorophyll )
                   zdiattot = ediatm(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
                   zprod = rday * (zpronewd(ji,jj,jk) + zproregd(ji,jj,jk)) * zprchld(ji,jj,jk) * xlimdia(ji,jj,jk)
@@ -454,5 +491,10 @@
         END DO
      END DO
-     !
+     
+     ! Production and uptake of ligands by phytoplankton. This part is activated 
+     ! when ln_ligand is set to .true. in the namelist. Ligand uptake is small 
+     ! and based on the FeL model by Morel et al. (2008) and on the study of
+     ! Shaked and Lis (2012)
+     ! -------------------------------------------------------------------------
      IF( ln_ligand ) THEN
          zpligprod1(:,:,:) = 0._wp    ;    zpligprod2(:,:,:) = 0._wp
@@ -463,14 +505,12 @@
                  zfeup    = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) + texcretp * zprofep(ji,jj,jk)
                  tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp    &
-                 &       - zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 2.E3 * (1.0 - plig(ji,jj,jk) ) )
+                 &       - zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
                  zpligprod1(ji,jj,jk) = zdocprod * ldocp
-                 zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) * lthet
+                 zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) &
+                 &                      + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
               END DO
            END DO
         END DO
      ENDIF
-
-
-     ! Total primary production per year
 
     ! Total primary production per year
@@ -490,5 +530,5 @@
               CALL iom_put( "PPPHYP"  , zw3d )
               !
-              zw3d(:,:,:) = zprorcad(:,:,:) * zfact * tmask(:,:,:)  ! primary production by diatomes
+              zw3d(:,:,:) = zprorcad(:,:,:) * zfact * tmask(:,:,:)  ! primary production by diatoms
               CALL iom_put( "PPPHYD"  , zw3d )
           ENDIF
@@ -500,5 +540,5 @@
               CALL iom_put( "PPNEWP"  , zw3d )
               !
-              zw3d(:,:,:) = zpronewd(:,:,:) * zfact * tmask(:,:,:)  ! new primary production by diatomes
+              zw3d(:,:,:) = zpronewd(:,:,:) * zfact * tmask(:,:,:)  ! new primary production by diatoms
               CALL iom_put( "PPNEWD"  , zw3d )
           ENDIF
@@ -508,20 +548,20 @@
           ENDIF
           IF( iom_use( "PFeN" ) .OR. iom_use( "PFeD" ) .OR. iom_use( "PFeP" ) )  THEN
-              zw3d(:,:,:) = zprofen(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron production by nanophyto
+              zw3d(:,:,:) = zprofen(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron uptake by nanophyto
               CALL iom_put( "PFeN"  , zw3d )
               !
-              zw3d(:,:,:) = zprofep(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron production by picophyto
+              zw3d(:,:,:) = zprofep(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron uptake by picophyto
               CALL iom_put( "PFeP"  , zw3d )
               !
-              zw3d(:,:,:) = zprofed(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron production by  diatomes
+              zw3d(:,:,:) = zprofed(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron uptake by  diatoms
               CALL iom_put( "PFeD"  , zw3d )
           ENDIF
           IF( iom_use( "LPRODP" ) )  THEN
               zw3d(:,:,:) = zpligprod1(:,:,:) * 1e9 * zfact * tmask(:,:,:)
-              CALL iom_put( "LPRODP"  , zw3d )
+              CALL iom_put( "LPRODP"  , zw3d )  ! Ligand production by phytoplankton
           ENDIF
           IF( iom_use( "LDETP" ) )  THEN
               zw3d(:,:,:) = zpligprod2(:,:,:) * 1e9 * zfact * tmask(:,:,:)
-              CALL iom_put( "LDETP"  , zw3d )
+              CALL iom_put( "LDETP"  , zw3d )  ! Uptake of ligands by phytoplankton
           ENDIF
           IF( iom_use( "Mumax" ) )  THEN
@@ -540,11 +580,11 @@
           ENDIF
           IF( iom_use( "LNlight" ) .OR. iom_use( "LDlight" ) .OR. iom_use( "LPlight" ) )  THEN
-              zw3d(:,:,:) = zprbio (:,:,:) / (zprmaxn(:,:,:) + rtrn) * tmask(:,:,:) ! light limitation term
+              zw3d(:,:,:) = zprbio (:,:,:) / (zprmaxn(:,:,:) + rtrn) * tmask(:,:,:) ! light limitation term of nanophytoplankton
               CALL iom_put( "LNlight"  , zw3d )
               !
-              zw3d(:,:,:) = zprpic (:,:,:) / (zprmaxp(:,:,:) + rtrn) * tmask(:,:,:) ! light limitation term
+              zw3d(:,:,:) = zprpic (:,:,:) / (zprmaxp(:,:,:) + rtrn) * tmask(:,:,:) ! light limitation term of picophytoplankton
               CALL iom_put( "LPlight"  , zw3d )
               !
-              zw3d(:,:,:) =  zprdia (:,:,:) / (zprmaxd(:,:,:) + rtrn) * tmask(:,:,:)  ! light limitation term
+              zw3d(:,:,:) =  zprdia (:,:,:) / (zprmaxd(:,:,:) + rtrn) * tmask(:,:,:)  ! light limitation term of diatoms
               CALL iom_put( "LDlight"  , zw3d )
           ENDIF
@@ -556,5 +596,5 @@
               CALL iom_put( "MunetP"  , zw3d )
               !
-              zw3d(:,:,:) = zcroissd(:,:,:) * tmask(:,:,:) ! ! Realized growth rate for diatomes
+              zw3d(:,:,:) = zcroissd(:,:,:) * tmask(:,:,:) ! ! Realized growth rate for diatoms
               CALL iom_put( "MunetD"  , zw3d )
               !
@@ -584,10 +624,10 @@
       !! ** Purpose :   Initialization of phytoplankton production parameters
       !!
-      !! ** Method  :   Read the nampisprod namelist and check the parameters
+      !! ** Method  :   Read the namp5zprod namelist and check the parameters
       !!      called at the first timestep (nittrc000)
       !!
-      !! ** input   :   Namelist nampisprod
+      !! ** input   :   Namelist namp5zprod
       !!----------------------------------------------------------------------
-      INTEGER :: ios                 ! Local integer output status for namelist read
+      INTEGER :: ios    ! Local integer output status for namelist read
       !!
       NAMELIST/namp5zprod/ pislopen, pislopep, pisloped, excretn, excretp, excretd,     &
@@ -595,9 +635,9 @@
       !!----------------------------------------------------------------------
 
-      REWIND( numnatp_ref )              ! Namelist nampisprod in reference namelist : Pisces phytoplankton production
+      REWIND( numnatp_ref )              ! Namelist namp5zprod in reference namelist : Pisces phytoplankton production
       READ  ( numnatp_ref, namp5zprod, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp5zprod in reference namelist' )
 
-      REWIND( numnatp_cfg )              ! Namelist nampisprod in configuration namelist : Pisces phytoplankton production
+      REWIND( numnatp_cfg )              ! Namelist namp5zprod in configuration namelist : Pisces phytoplankton production
       READ  ( numnatp_cfg, namp5zprod, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp5zprod in configuration namelist' )
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/par_pisces.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/par_pisces.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/par_pisces.F90	(revision 12537)
@@ -10,6 +10,6 @@
 
    ! productive layer depth
-   INTEGER, PUBLIC ::   jpkb       !: first vertical layers where biology is active
-   INTEGER, PUBLIC ::   jpkbm1     !: first vertical layers where biology is active
+   INTEGER, PUBLIC ::   jpkb      !: first vertical layers where biology is active
+   INTEGER, PUBLIC ::   jpkbm1    !: first vertical layers where biology is active
 
    ! assign an index in trc arrays for each LOBSTER prognostic variables
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/sms_pisces.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/sms_pisces.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/sms_pisces.F90	(revision 12537)
@@ -2,5 +2,5 @@
    !!----------------------------------------------------------------------
    !!                     ***  sms_pisces.F90  ***  
-   !! TOP :   PISCES Source Minus Sink variables
+   !! TOP :   PISCES Source Minus Sink variables are declared and allocated
    !!----------------------------------------------------------------------
    !! History :   1.0  !  2000-02 (O. Aumont) original code
@@ -13,9 +13,7 @@
    PUBLIC
 
-   INTEGER ::   numnatp_ref = -1           !! Logical units for namelist pisces
-   INTEGER ::   numnatp_cfg = -1           !! Logical units for namelist pisces
-   INTEGER ::   numonp      = -1           !! Logical unit for namelist pisces output
-
-   !                                                       !:  PISCES  : silicon dependant half saturation
+   INTEGER ::   numnatp_ref = -1  !: Logical units for namelist pisces
+   INTEGER ::   numnatp_cfg = -1  !: Logical units for namelist pisces
+   INTEGER ::   numonp      = -1  !: Logical unit for namelist pisces output
 
    !!* Model used
@@ -27,45 +25,43 @@
 
    !!*  Time variables
-   INTEGER  ::   nrdttrc           !: ???
-   REAL(wp) ::   rfact , rfactr    !: ???
-   REAL(wp) ::   rfact2, rfact2r   !: ???
-   REAL(wp) ::   xstep             !: Time step duration for biology
-   REAL(wp) ::   ryyss             !: number of seconds per year 
-   REAL(wp) ::   r1_ryyss          !: inverse number of seconds per year 
-
+   INTEGER  ::   nrdttrc          !: ???
+   REAL(wp) ::   rfact , rfactr   !: time step duration (in seconds)
+   REAL(wp) ::   rfact2, rfact2r  !: time step duration (in seconds) when timesplitting is activated for PISCES
+   REAL(wp) ::   xstep            !: Time step duration for biology
+   REAL(wp) ::   ryyss            !: number of seconds per year 
+   REAL(wp) ::   r1_ryyss         !: inverse number of seconds per year 
 
    !!*  Biological parameters 
-   REAL(wp) ::   rno3              !: ???
-   REAL(wp) ::   o2ut              !: ???
-   REAL(wp) ::   po4r              !: ???
-   REAL(wp) ::   rdenit            !: ???
-   REAL(wp) ::   rdenita           !: ???
-   REAL(wp) ::   o2nit             !: ???
-   REAL(wp) ::   wsbio, wsbio2     !: ???
-   REAL(wp) ::   wsbio2max         !: ???
-   REAL(wp) ::   wsbio2scale       !: ???
-   REAL(wp) ::   xkmort            !: ???
-   REAL(wp) ::   ferat3            !: ???
-   REAL(wp) ::   ldocp             !: ???
-   REAL(wp) ::   ldocz             !: ???
-   REAL(wp) ::   lthet             !: ???
-   REAL(wp) ::   no3rat3           !: ???
-   REAL(wp) ::   po4rat3           !: ???
-
+   REAL(wp) ::   rno3             !: C/N stoichiometric ratio
+   REAL(wp) ::   o2ut             !: O2/N stoichiometric ratio for ammonification
+   REAL(wp) ::   po4r             !: C/P stoichiometric ratio
+   REAL(wp) ::   rdenit           !: C/N ratio for denitrification
+   REAL(wp) ::   rdenita          !: C/N ratio for denitrification
+   REAL(wp) ::   o2nit            !: O2/N ratio for nitrification
+   REAL(wp) ::   wsbio, wsbio2    !: Sinking speeds of particles
+   REAL(wp) ::   wsbio2max        !: Maximum sinking speed of the largest particles
+   REAL(wp) ::   wsbio2scale      !: Length scale for the variations of wsbio2
+   REAL(wp) ::   xkmort           !: Mortality half-saturation constant
+   REAL(wp) ::   ferat3           !: Fe/C in heterotrophic organisms
+   REAL(wp) ::   ldocp            !: Ligand production ratio during PP
+   REAL(wp) ::   ldocz            !: Ligand production ratio by grazing
+   REAL(wp) ::   lthet            !: Uptake of ligand by phytoplankton
+   REAL(wp) ::   no3rat3          !: C/N ratio of zooplankton
+   REAL(wp) ::   po4rat3          !: C/P ratio of zooplankton
 
    !!*  diagnostic parameters 
-   REAL(wp) ::  tpp                !: total primary production
-   REAL(wp) ::  t_oce_co2_exp      !: total carbon export
-   REAL(wp) ::  t_oce_co2_flx      !: Total ocean carbon flux
-   REAL(wp) ::  t_oce_co2_flx_cum  !: Cumulative Total ocean carbon flux
-   REAL(wp) ::  t_atm_co2_flx      !: global mean of atmospheric pco2
+   REAL(wp) ::  tpp               !: total primary production
+   REAL(wp) ::  t_oce_co2_exp     !: total carbon export
+   REAL(wp) ::  t_oce_co2_flx     !: Total ocean carbon flux
+   REAL(wp) ::  t_oce_co2_flx_cum !: Cumulative Total ocean carbon flux
+   REAL(wp) ::  t_atm_co2_flx     !: global mean of atmospheric pco2
 
    !!* restoring
-   LOGICAL  ::  ln_pisdmp          !: restoring or not of nutrients to a mean value
-   INTEGER  ::  nn_pisdmp          !: frequency of relaxation or not of nutrients to a mean value
+   LOGICAL  ::  ln_pisdmp         !: restoring or not of nutrients to a mean value
+   INTEGER  ::  nn_pisdmp         !: frequency of relaxation or not of nutrients to a mean value
 
    !!* Mass conservation
-   LOGICAL  ::  ln_check_mass      !: Flag to check mass conservation
-   LOGICAL , PUBLIC ::   ln_ironice   !: boolean for Fe input from sea ice
+   LOGICAL  ::  ln_check_mass     !: Flag to check mass conservation
+   LOGICAL, PUBLIC ::   ln_ironice   !: boolean for Fe input from sea ice
 
    !!*  Biological fluxes for light : variables shared by pisces & lobster
@@ -75,13 +71,13 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  etot_ndcy      !: PAR over 24h in case of diurnal cycle
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  enano, ediat   !: PAR for phyto, nano and diat 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  enanom, ediatm !: PAR for phyto, nano and diat 
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  enanom, ediatm !: mean PAR for phyto, nano and diat 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  epico          !: PAR for pico
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  epicom         !: PAR for pico
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  epicom         !: mean PAR for pico
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  emoy           !: averaged PAR in the mixed layer
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  heup_01 !: Absolute euphotic layer depth
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  xksi  !:  LOBSTER : zooplakton closure
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  xksi  !:  Half-saturation con,stant for diatoms
 
    !!*  Biological fluxes for primary production
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   xksimax    !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   xksimax    !: Maximum half-saturation constant over the year (Si)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   biron      !: bioavailable fraction of iron
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   plig       !: proportion of iron organically complexed
@@ -92,14 +88,14 @@
 
    !!*  SMS for the organic matter
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xfracal    !: ??
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   nitrfac    !: ??
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   nitrfac2   !: ??
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   orem       !: ??
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xdiss      !: ??
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xfracal    !: Fraction of nanophytoplankton that are calcifying organisms
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   nitrfac    !: OMZ 
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   nitrfac2   !: N depleted indice
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   orem       !: oxic remineralisation
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xdiss      !: Shear rate
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodcal    !: Calcite production
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodpoc    !: Calcite production
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   conspoc    !: Calcite production
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodgoc    !: Calcite production
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   consgoc    !: Calcite production
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodpoc    !: POC production
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   conspoc    !: POC consumption
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodgoc    !: GOC production
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   consgoc    !: GOC consumption
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   blim       !: bacterial production factor
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizen      !: size of nanophyto
@@ -111,9 +107,9 @@
 
    !!* Variable for chemistry of the CO2 cycle
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak13       !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak23       !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aksp       !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hi         !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   excess     !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak13       !: Carbonate chemistry constant
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak23       !: Carbonate chemistry constant
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aksp       !: Solubility product of CaCO3
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hi         !: Proton concentration
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   excess     !: CO3 saturation
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aphscale   !: 
 
@@ -139,5 +135,4 @@
       !!----------------------------------------------------------------------
       ierr(:) = 0
-      !*  Biological fluxes for light : shared variables for pisces & lobster
       ALLOCATE( etot(jpi,jpj,jpk), neln(jpi,jpj), heup(jpi,jpj),    &
         &       heup_01(jpi,jpj) , xksi(jpi,jpj)               ,  STAT=ierr(1) )
@@ -145,15 +140,16 @@
   
       IF( ln_p4z .OR. ln_p5z ) THEN
-         !*  Biological fluxes for light 
+
+         !* Optics
          ALLOCATE(  enano(jpi,jpj,jpk)    , ediat(jpi,jpj,jpk) ,   &
            &        enanom(jpi,jpj,jpk)   , ediatm(jpi,jpj,jpk),   &
            &        etot_ndcy(jpi,jpj,jpk), emoy(jpi,jpj,jpk)  ,  STAT=ierr(2) ) 
 
-         !*  Biological fluxes for primary production
+         !* Biological SMS
          ALLOCATE( xksimax(jpi,jpj)  , biron(jpi,jpj,jpk)      ,  STAT=ierr(3) )
-         !
-         !*  SMS for the organic matter
+
+         ! Biological SMS
          ALLOCATE( xfracal (jpi,jpj,jpk), orem(jpi,jpj,jpk)    ,    &
-            &      nitrfac(jpi,jpj,jpk), nitrfac2(jpi,jpj,jpk) ,    &
+            &      nitrfac(jpi,jpj,jpk) , nitrfac2(jpi,jpj,jpk),    &
             &      prodcal(jpi,jpj,jpk) , xdiss   (jpi,jpj,jpk),    &
             &      prodpoc(jpi,jpj,jpk) , conspoc(jpi,jpj,jpk) ,    &
@@ -161,18 +157,19 @@
             &      blim   (jpi,jpj,jpk) ,                         STAT=ierr(4) )
 
-         !* Variable for chemistry of the CO2 cycle
-         ALLOCATE( ak13  (jpi,jpj,jpk) ,                            &
+         !* Carbonate chemistry
+         ALLOCATE( ak13  (jpi,jpj,jpk)  ,                           &
             &      ak23(jpi,jpj,jpk)    , aksp  (jpi,jpj,jpk) ,     &
             &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,     &
             &      aphscale(jpi,jpj,jpk),                         STAT=ierr(5) )
          !
-         !* Temperature dependancy of SMS terms
+         !* Temperature dependency of SMS terms
          ALLOCATE( tgfunc (jpi,jpj,jpk) , tgfunc2(jpi,jpj,jpk),     &
             &      tgfunc3(jpi,jpj,jpk) , STAT=ierr(6) )
          !
-         !* Sinkong speed
+         !* Sinking speed
          ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4 (jpi,jpj,jpk),     &
             &                             STAT=ierr(7) )   
          ! 
+         !* Prognostic ligand
          IF( ln_ligand ) THEN
            ALLOCATE( plig(jpi,jpj,jpk)  ,                         STAT=ierr(8) )
@@ -181,5 +178,5 @@
       !
       IF( ln_p5z ) THEN
-         !       
+         ! PISCES-QUOTA specific part      
          ALLOCATE( epico(jpi,jpj,jpk)   , epicom(jpi,jpj,jpk) ,   STAT=ierr(9) ) 
 
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/trcini_pisces.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/trcini_pisces.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/trcini_pisces.F90	(revision 12537)
@@ -3,4 +3,5 @@
    !!                         ***  MODULE trcini_pisces  ***
    !! TOP :   initialisation of the PISCES biochemical model
+   !!         This module is for LOBSTER, PISCES and PISCES-QUOTA
    !!======================================================================
    !! History :    -   !  1988-07  (E. Maier-Reiner) Original code
@@ -37,6 +38,8 @@
       !!
       !! ** Purpose :   Initialisation of the PISCES biochemical model
-      !!----------------------------------------------------------------------
-      !
+      !!                Allocation of the dynamic arrays
+      !!----------------------------------------------------------------------
+      !
+      ! Read the PISCES namelist
       CALL trc_nam_pisces
       !
@@ -54,5 +57,5 @@
       !! ** Purpose :   Initialisation of the PISCES biochemical model
       !!----------------------------------------------------------------------
-      USE p4zsms          ! Main P4Z routine
+      USE p4zsms          !  Main P4Z routine
       USE p4zche          !  Chemical model
       USE p4zsink         !  vertical flux of particulate matter due to sinking
@@ -111,7 +114,9 @@
       ierr = ierr +  p4z_lim_alloc()
       IF( ln_p4z ) THEN
+         ! PISCES part
          ierr = ierr +  p4z_prod_alloc()
          ierr = ierr +  p4z_meso_alloc()
       ELSE
+         ! PISCES-QUOTA part
          ierr = ierr +  p5z_lim_alloc()
          ierr = ierr +  p5z_prod_alloc()
@@ -128,4 +133,5 @@
 
       ! assign an index in trc arrays for each prognostic variables
+      ! This is based on the information read in the namelist_top
       DO jn = 1, jptra
         cltra = ctrcnm(jn) 
@@ -177,13 +183,13 @@
       ! Set biological ratios
       ! ---------------------
-      rno3    =  16._wp / 122._wp
-      po4r    =   1._wp / 122._wp
-      o2nit   =  32._wp / 122._wp
-      o2ut    = 133._wp / 122._wp
-      rdenit  =  ( ( o2ut + o2nit ) * 0.80 - rno3 - rno3 * 0.60 ) / rno3
-      rdenita =   3._wp /  5._wp
+      rno3    =  16._wp / 122._wp   ! C/N
+      po4r    =   1._wp / 122._wp   ! C/P
+      o2nit   =  32._wp / 122._wp   ! O2/C for nitrification
+      o2ut    = 133._wp / 122._wp   ! O2/C for ammonification
+      rdenit  =  ( ( o2ut + o2nit ) * 0.80 - rno3 - rno3 * 0.60 ) / rno3  ! Denitrification
+      rdenita =   3._wp /  5._wp    ! Denitrification
       IF( ln_p5z ) THEN
-         no3rat3 = no3rat3 / rno3
-         po4rat3 = po4rat3 / po4r
+         no3rat3 = no3rat3 / rno3   ! C/N ratio in zooplankton
+         po4rat3 = po4rat3 / po4r   ! C/P ratio in zooplankton
       ENDIF
 
@@ -250,7 +256,9 @@
       CALL p4z_opt_init          !  Optic: PAR in the water column
       IF( ln_p4z ) THEN
+         ! PISCES part
          CALL p4z_lim_init       !  co-limitations by the various nutrients
          CALL p4z_prod_init      !  phytoplankton growth rate over the global ocean.
       ELSE
+         ! PISCES-QUOTA part
          CALL p5z_lim_init       !  co-limitations by the various nutrients
          CALL p5z_prod_init      !  phytoplankton growth rate over the global ocean.
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/trcnam_pisces.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/trcnam_pisces.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/trcnam_pisces.F90	(revision 12537)
@@ -34,8 +34,10 @@
       !!                     ***  trc_nam_pisces  ***  
       !!
-      !! ** Purpose :   read PISCES namelist
+      !! ** Purpose :  read PISCES namelist. The parameters read here activate
+      !!               different functionalities of PISCES : PISCES, PISCES-QUOTA
+      !!               LOBSTER, variable ligands, sediment module
       !!
-      !! ** input   :   file 'namelist.trc.sms' containing the following
-      !!             namelist: natext, natbio, natsms
+      !! ** input   :  file 'namelist_pisces' containing the following
+      !!             namelist: nampismod
       !!----------------------------------------------------------------------
       INTEGER :: jl, jn
@@ -55,8 +57,8 @@
       IF(lwm) CALL ctl_opn( numonp     , 'output.namelist.pis' , 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
       !
-      REWIND( numnatp_ref )              ! Namelist nampisbio in reference namelist : Pisces variables
+      REWIND( numnatp_ref )              ! Namelist nampismod in reference namelist : Pisces variables
       READ  ( numnatp_ref, nampismod, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampismod in reference namelist' )
-      REWIND( numnatp_cfg )              ! Namelist nampisbio in configuration namelist : Pisces variables
+      REWIND( numnatp_cfg )              ! Namelist nampismod in configuration namelist : Pisces variables
       READ  ( numnatp_cfg, nampismod, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampismod in configuration namelist' )
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/trcsms_pisces.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/trcsms_pisces.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/trcsms_pisces.F90	(revision 12537)
@@ -3,4 +3,5 @@
    !!                         ***  MODULE trcsms_pisces  ***
    !! TOP :   PISCES Source Minus Sink manager
+   !!         This module is for LOBSTER, PISCES and PISCES-QUOTA
    !!======================================================================
    !! History :   1.0  !  2004-03 (O. Aumont) Original code
@@ -30,5 +31,5 @@
       !!
       !! ** Purpose :   Managment of the call to Biological sources and sinks 
-      !!                routines of PISCES or LOBSTER bio-model
+      !!                routines of PISCES/PISCES-QUOTA or LOBSTER bio-model
       !!
       !!---------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcbc.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcbc.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcbc.F90	(revision 12537)
@@ -153,5 +153,5 @@
                IF(trcdta_bdy(jn,ib)%cn_obc == 'frs' .AND. nn_trcdmp_bdy(ib) /= 0 )  &
                    & CALL ctl_stop( 'trc_bc_ini: Use FRS OR relaxation' )
-               IF(  .NOT.( 0 < nn_trcdmp_bdy(ib)  .AND.  nn_trcdmp_bdy(ib) <= 2 )  )   &
+               IF( .NOT.( 0 <= nn_trcdmp_bdy(ib) .AND. nn_trcdmp_bdy(ib) <= 2 ) ) &
                    & CALL ctl_stop( 'trc_bc_ini: Not a valid option for nn_trcdmp_bdy. Allowed: 0,1,2.' )
             END DO
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcini.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcini.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcini.F90	(revision 12537)
@@ -235,5 +235,5 @@
       IF( ln_trcdta )   CALL trc_dta_ini( jptra )      ! set initial tracers values
       !
-      IF( ln_my_trc )   CALL trc_bc_ini ( jptra )      ! set tracers Boundary Conditions
+      CALL trc_bc_ini ( jptra )      ! set tracers Boundary Conditions
       !
       !
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcsms.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcsms.F90	(revision 12536)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcsms.F90	(revision 12537)
@@ -21,4 +21,6 @@
    USE trcsms_my_trc      ! MY_TRC  tracers
    USE prtctl_trc         ! Print control for debbuging
+   USE trcbc  , only : trc_bc
+   USE bdy_oce, only : ln_bdy 
 
    IMPLICIT NONE
@@ -55,4 +57,6 @@
       IF( ln_my_trc  )   CALL trc_sms_my_trc ( kt )    ! MY_TRC  tracers
 
+      IF( ln_bdy     )   CALL trc_bc ( kt )
+
       IF(ln_ctl) THEN      ! print mean trends (used for debugging)
          WRITE(charout, FMT="('sms ')")
