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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zagg.F90	(revision 13233)
@@ -100,4 +100,5 @@
                   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
+                  ! The temperature dependency has been omitted.
                   zaggdoc3 =  114. * 0.3 * trb(ji,jj,jk,jpdoc) *xstep * 0.3 * trb(ji,jj,jk,jpdoc)
 
@@ -123,5 +124,5 @@
                   !
                   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
                   zaggtmp = 25.9  * zfact * trb(ji,jj,jk,jppoc)
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 13233)
@@ -66,5 +66,6 @@
 
       ! ASSIGN THE SHEAR RATE THAT IS USED FOR AGGREGATION
-      ! OF PHYTOPLANKTON AND DETRITUS
+      ! OF PHYTOPLANKTON AND DETRITUS. Shear rate is supposed to equal 1
+      ! in the mixed layer and 0.1 below the mixed layer.
       xdiss(:,:,:) = 1.
       DO jk = 2, jpkm1
@@ -80,18 +81,20 @@
       CALL p4z_fechem  ( kt, knt )     ! Iron chemistry/scavenging
       !
-      IF( ln_p4z ) THEN
+      IF( ln_p4z ) THEN  ! PISCES standard
+         ! Phytoplankton only sources/sinks terms
          CALL p4z_lim  ( kt, knt )     ! co-limitations by the various nutrients
          CALL p4z_prod ( kt, knt )     ! phytoplankton growth rate over the global ocean. 
          !                             ! (for each element : C, Si, Fe, Chl )
          CALL p4z_mort ( kt      )     ! phytoplankton mortality
-         !                             ! zooplankton sources/sinks routines 
+         ! zooplankton sources/sinks routines 
          CALL p4z_micro( kt, knt )     ! microzooplankton
          CALL p4z_meso ( kt, knt )     ! mesozooplankton
-      ELSE
+      ELSE  ! PISCES-QUOTA
+         ! Phytoplankton only sources/sinks terms
          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, N, P, Si, Fe, Chl )
          CALL p5z_mort ( kt      )     ! phytoplankton mortality
-         !                             ! zooplankton sources/sinks routines 
+         ! zooplankton sources/sinks routines 
          CALL p5z_micro( kt, knt )     ! microzooplankton
          CALL p5z_meso ( kt, knt )     ! mesozooplankton
@@ -99,5 +102,5 @@
       !
       CALL p4z_agg     ( kt, knt )     ! Aggregation of particles
-      CALL p4z_rem     ( kt, knt )     ! remineralization terms of organic matter+scavenging of Fe
+      CALL p4z_rem     ( kt, knt )     ! remineralization terms of organic matter
       CALL p4z_poc     ( kt, knt )     ! Remineralization of organic particles
       !
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlim.F90	(revision 13233)
@@ -84,4 +84,6 @@
       !!
       !! ** Method  : - Limitation follows the Liebieg law of the minimum
+      !!              - Monod approach for N, P and Si. Quota approach 
+      !!                for Iron
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in)  :: kt, knt
@@ -113,4 +115,5 @@
                ! Computation of a variable Ks of diatoms taking into account
                ! that increasing biomass is made of generally bigger cells
+               ! The allometric relationship is classical.
                !------------------------------------------------------------
                z1_trbphy   = 1. / ( trb(ji,jj,jk,jpphy) + rtrn )
@@ -156,6 +159,10 @@
 
                ! Michaelis-Menten Limitation term by nutrients: Nanophyto
+               ! Optimal parameterization by Smith and Pahlow series of 
+               ! papers is used. Optimal allocation is supposed independant
+               ! for all nutrients. 
                ! --------------------------------------------------------
-               ! Limitation of Fe uptake
+
+               ! Limitation of Fe uptake (Quota formalism)
                zfalim = (1.-fananof) / fananof
                xnanofer(ji,jj,jk) = (1. - fananof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * concnfe(ji,jj,jk) )
@@ -181,5 +188,5 @@
                !   Michaelis-Menten Limitation term by nutrients : Diatoms
                !   -------------------------------------------------------
-               ! Limitation of Fe uptake
+               ! Limitation of Fe uptake (Quota formalism)
                zfalim = (1.-fadiatf) / fadiatf
                xdiatfer(ji,jj,jk) = (1. - fadiatf) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * concdfe(ji,jj,jk) )
@@ -256,4 +263,6 @@
             DO ji = 1, jpi
                ! denitrification factor computed from O2 levels
+               ! This factor diagnoses below which level of O2 denitrification
+               ! is active
                nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    &
                   &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  )
@@ -261,4 +270,6 @@
                !
                ! redox factor computed from NO3 levels
+               ! This factor diagnoses below which level of NO3 additional redox
+               ! reactions are taking place.
                nitrfac2(ji,jj,jk) = MAX( 0.e0,       ( 1.E-6 - trb(ji,jj,jk,jpno3) )  &
                   &                                / ( 1.E-6 + trb(ji,jj,jk,jpno3) ) )
@@ -296,5 +307,6 @@
       !!----------------------------------------------------------------------
       INTEGER ::   ios   ! Local integer
-      !
+
+      ! Namelist block
       NAMELIST/namp4zlim/ concnno3, concdno3, concnnh4, concdnh4, concnfer, concdfer, concbfe,   &
          &                concbno3, concbnh4, xsizedia, xsizephy, xsizern, xsizerd,          & 
@@ -348,4 +360,6 @@
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE p5z_lim_alloc  ***
+      !! 
+      !            Allocation of the arrays used in this module
       !!----------------------------------------------------------------------
       USE lib_mpp , ONLY: ctl_stop
Index: NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlys.F90
===================================================================
--- NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlys.F90	(revision 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlys.F90	(revision 13233)
@@ -51,5 +51,6 @@
       !!                OF CACO3 TO THE CACO3 SEDIMENT POOL.
       !!
-      !! ** Method  : - ???
+      !! ** Method  : - pH is computed using the MOCSY module and dissolution
+      !!                is a function of the saturation state of CO3
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, knt   ! ocean time step and ???
@@ -71,5 +72,6 @@
       !     COMPUTE [CO3--] and [H+] CONCENTRATIONS
       !     -------------------------------------------
-
+      !     Call the carbonate chemistry subroutine to compute pH
+      !     This subroutine is in the p4zche module
       CALL solve_at_general( zhinit, zhi )
 
@@ -85,7 +87,7 @@
 
       !     ---------------------------------------------------------
-      !        CALCULATE DEGREE OF CACO3 SATURATION AND CORRESPONDING
-      !        DISSOLOUTION AND PRECIPITATION OF CACO3 (BE AWARE OF
-      !        MGCO3)
+      !     CALCULATE DEGREE OF CACO3 SATURATION AND CORRESPONDING
+      !     DISSOLOUTION AND PRECIPITATION OF CACO3 (BE AWARE OF
+      !     MGCO3)
       !     ---------------------------------------------------------
 
@@ -106,20 +108,18 @@
                zexcess  = zexcess0**nca
 
-               ! AMOUNT CACO3 (12C) THAT RE-ENTERS SOLUTION
-               !       (ACCORDING TO THIS FORMULATION ALSO SOME PARTICULATE
-               !       CACO3 GETS DISSOLVED EVEN IN THE CASE OF OVERSATURATION)
+               ! AMOUNT CACO3 THAT RE-ENTERS SOLUTION
                zdispot = kdca * zexcess * trb(ji,jj,jk,jpcal)
-              !  CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
-              !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
-              zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss ! calcite dissolution
-              !
-              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk)
-              tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) -      zcaldiss(ji,jj,jk)
-              tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +      zcaldiss(ji,jj,jk)
+               ! CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
+               ! AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
+               zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss ! calcite dissolution
+               !
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk)
+               tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) -      zcaldiss(ji,jj,jk)
+               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +      zcaldiss(ji,jj,jk)
             END DO
          END DO
       END DO
-      !
 
+      ! output of diagnostics
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
          IF( iom_use( "PH"     ) ) CALL iom_put( "PH"    , -1. * LOG10( MAX( hi(:,:,:), rtrn ) ) * tmask(:,:,:) )
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 13233)
@@ -45,4 +45,6 @@
    REAL(wp), PUBLIC ::  epsher2      !: growth efficiency
    REAL(wp), PUBLIC ::  epsher2min   !: minimum growth efficiency at high food for grazing 2
+   REAL(wp), PUBLIC ::  xsigma2      !: Width of the predation window
+   REAL(wp), PUBLIC ::  xsigma2del   !: Maximum width of the predation window at low food density
    REAL(wp), PUBLIC ::  grazflux     !: mesozoo flux feeding rate
    REAL(wp), PUBLIC ::  xfracmig     !: Fractional biomass of meso that performs DVM
@@ -169,5 +171,5 @@
                ! ----------------------------------------------------------
                zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
-               zsigma = 0.5 + 1.0 * zsigma
+               zsigma = xsigma2 + xsigma2del * zsigma
                ! Nanophytoplankton and diatoms are the only preys considered
                ! to be close enough to have potential interference
@@ -186,9 +188,10 @@
                !   Mesozooplankton regular grazing on the different preys
                !   ------------------------------------------------------
-               zgrazdc   = zgraze2  * ztmp3 * zdenom
-               zgraznc   = zgraze2  * ztmp1 * zdenom
-               zgrazpoc  = zgraze2  * ztmp2 * zdenom
-               zgrazz    = zgraze2  * ztmp4 * zdenom
-
+               zgrazdc   = zgraze2  * ztmp3 * zdenom  ! diatoms
+               zgraznc   = zgraze2  * ztmp1 * zdenom  ! nanophytoplankton
+               zgrazpoc  = zgraze2  * ztmp2 * zdenom  ! small POC
+               zgrazz    = zgraze2  * ztmp4 * zdenom  ! microzooplankton
+
+               ! Ingestion rates of the Fe content of the different preys
                zgraznf   = zgraznc  * trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) + rtrn)
                zgrazdf   = zgrazdc  * trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) + rtrn)
@@ -228,4 +231,6 @@
                zgrazfffp = zproport * zgrazfffp
                zgrazfffg = zproport * zgrazfffg
+
+               ! Total ingestion rates in C, N, Fe
                zgraztotc = zgrazdc + zgrazz + zgraznc + zgrazpoc + zgrazffep + zgrazffeg
                zgraztotn = zgrazdc * quotad(ji,jj,jk) + zgrazz + zgraznc * quotan(ji,jj,jk)   &
@@ -248,6 +253,9 @@
                zepshert  = MIN( 1., zgrasratn, zgrasratf / ferat3)
                zbeta     = MAX(0., (epsher2 - epsher2min) )
+               ! Food quantity deprivation of GGE
                zepsherf  = epsher2min + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta ) 
+               ! Food quality deprivation of GGE
                zepsherq  = 0.5 + (1.0 - 0.5) * zepshert * ( 1.0 + 1.0 ) / ( zepshert + 1.0 )
+               ! Actual GGE
                zepsherv  = zepsherf * zepshert * zepsherq
                ! 
@@ -258,4 +266,6 @@
                ! according to a infinite chain of predators (ANderson et al., 2013)
                zmortzgoc = unass2 / ( 1. - epsher2 ) * ztortz + zrespz
+
+               ! Update of the trends
                tra(ji,jj,jk,jpmes) = tra(ji,jj,jk,jpmes) - zmortz + zepsherv * zgraztotc
                tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazdc
@@ -502,4 +512,6 @@
          WRITE(numout,*) '      Fraction of mesozoo excretion as DOM           sigma2       =', sigma2
          WRITE(numout,*) '      half sturation constant for grazing 2          xkgraz2      =', xkgraz2
+         WRITE(numout,*) '      Width of the grazing window                     xsigma2     =', xsigma2
+         WRITE(numout,*) '      Maximum additional width of the grazing window  xsigma2del  =', xsigma2del
          WRITE(numout,*) '      Diurnal vertical migration of mesozoo.         ln_dvm_meso  =', ln_dvm_meso
          WRITE(numout,*) '      Fractional biomass of meso  that performs DVM  xfracmig     =', xfracmig
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 13233)
@@ -42,4 +42,7 @@
    REAL(wp), PUBLIC ::   epsher      !: growth efficiency for grazing 1 
    REAL(wp), PUBLIC ::   epshermin   !: minimum growth efficiency for grazing 1
+   REAL(wp), PUBLIC ::   xsigma      !: Width of the grazing window
+   REAL(wp), PUBLIC ::   xsigmadel   !: Maximum additional width of the grazing window at low food density 
+
 
    !!----------------------------------------------------------------------
@@ -59,5 +62,5 @@
       !!                parameterization.
       !!                All living compartments smaller than microzooplankton
-      !!                are potential preys of mesozooplankton
+      !!                are potential preys of microzooplankton
       !!
       !! ** Method  : - ???
@@ -146,5 +149,5 @@
                ! ----------------------------------------------------------
                zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
-               zsigma = 0.5 + 1.0*zsigma
+               zsigma = xsigma + xsigmadel * zsigma
                zdiffdn = exp( -ABS(log(1.67 * sizen(ji,jj,jk) / (5.0 * sized(ji,jj,jk) + rtrn )) )**2 / zsigma**2)
                ztmp1 = xprefn * zcompaph * ( zcompaph + zdiffdn * zcompadi ) / ( 1.0 + zdiffdn )
@@ -156,12 +159,15 @@
                ztmp3 = ztmp3 / ztmptot
 
-               zgraznc   = zgraze   * ztmp1 * zdenom
-               zgrazdc   = zgraze   * ztmp2 * zdenom
-               zgrazpoc  = zgraze   * ztmp3 * zdenom
-
+               ! Ingestion terms on the different preys of microzooplankton
+               zgraznc   = zgraze   * ztmp1 * zdenom  ! Nanophytoplankton
+               zgrazdc   = zgraze   * ztmp2 * zdenom  ! Diatoms
+               zgrazpoc  = zgraze   * ztmp3 * zdenom  ! POC
+
+               ! Ingestion terms on the iron content of the different preys
                zgraznf   = zgraznc  * trb(ji,jj,jk,jpnfe) / (trb(ji,jj,jk,jpphy) + rtrn)
                zgrazpof  = zgrazpoc * trb(ji,jj,jk,jpsfe) / (trb(ji,jj,jk,jppoc) + rtrn)
                zgrazdf   = zgrazdc  * trb(ji,jj,jk,jpdfe) / (trb(ji,jj,jk,jpdia) + rtrn)
                !
+               ! Total ingestion rate in C, Fe, N units
                zgraztotc = zgraznc + zgrazpoc + zgrazdc
                zgraztotf = zgraznf + zgrazdf  + zgrazpof 
@@ -183,13 +189,20 @@
                zepshert  =  MIN( 1., zgrasratn, zgrasratf / ferat3)
                zbeta     = MAX(0., (epsher - epshermin) )
+               ! Food quantity deprivation of the GGE
                zepsherf  = epshermin + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
+               ! Food quality deprivation of the GGE
                zepsherq  = 0.5 + (1.0 - 0.5) * zepshert * ( 1.0 + 1.0 ) / ( zepshert + 1.0 )
+               ! Actual GGE of microzooplankton
                zepsherv  = zepsherf * zepshert * zepsherq
+               ! Excretion of Fe
                zgrafer   = zgraztotc * MAX( 0. , ( 1. - unass ) * zgrasratf - ferat3 * zepsherv ) 
+               ! Excretion of C, N, P
                zgrarem   = zgraztotc * ( 1. - zepsherv - unass )
+               ! Egestion of C, N, P
                zgrapoc   = zgraztotc * unass
 
                !  Update of the TRA arrays
                !  ------------------------
+               ! Fraction of excretion as inorganic nutrients and DIC
                zgrarsig  = zgrarem * sigma1
                tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zgrarsig
@@ -287,5 +300,6 @@
       NAMELIST/namp4zzoo/ part, grazrat, resrat, mzrat, xprefn, xprefc, &
          &                xprefd,  xthreshdia,  xthreshphy,  xthreshpoc, &
-         &                xthresh, xkgraz, epsher, epshermin, sigma1, unass
+         &                xthresh, xkgraz, epsher, epshermin, sigma1, unass  &
+         &                xsigma, xsigmadel
       !!----------------------------------------------------------------------
       !
@@ -321,5 +335,8 @@
          WRITE(numout,*) '      Minimum efficicency of microzoo growth          epshermin   =', epshermin
          WRITE(numout,*) '      Fraction of microzoo excretion as DOM           sigma1      =', sigma1
-         WRITE(numout,*) '      half sturation constant for grazing 1           xkgraz      =', xkgraz
+         WRITE(numout,*) '      half saturation constant for grazing 1          xkgraz      =', xkgraz
+         WRITE(numout,*) '      Width of the grazing window                     xsigma      =', xsigma
+         WRITE(numout,*) '      Maximum additional width of the grazing window  xsigmadel   =', xsigmadel
+
       ENDIF
       !
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 13233)
@@ -48,17 +48,17 @@
       !!---------------------------------------------------------------------
       !
-      CALL p4z_nano            ! nanophytoplankton
-      CALL p4z_diat            ! diatoms
+      CALL p4z_mort_nano            ! nanophytoplankton
+      CALL p4z_mort_diat            ! diatoms
       !
    END SUBROUTINE p4z_mort
 
 
-   SUBROUTINE p4z_nano
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE p4z_nano  ***
+   SUBROUTINE p4z_mort_nano
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_mort_nano  ***
       !!
       !! ** Purpose :   Compute the mortality terms for nanophytoplankton
       !!
-      !! ** Method  : - ???
+      !! ** Method  :   Both quadratic and simili linear mortality terms
       !!---------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk
@@ -69,5 +69,5 @@
       !!---------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('p4z_nano')
+      IF( ln_timing )   CALL timing_start('p4z_mort_nano')
       !
       prodcal(:,:,:) = 0._wp   ! calcite production variable set to zero
@@ -101,7 +101,11 @@
                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)
+
+               ! POC associated with the shell is supposed to be routed to 
+               ! big particles because of the ballasting effect
                zfracal = 0.5 * xfracal(ji,jj,jk)
                tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprcaca
@@ -112,4 +116,6 @@
                prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + ( 1. - zfracal ) * zmortp
                prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zfracal * zmortp
+
+               ! Update the arrays TRA which contains the biological sources and sinks
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ( 1. - zfracal ) * zmortp * zfactfe
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zfracal * zmortp * zfactfe
@@ -124,16 +130,16 @@
        ENDIF
       !
-      IF( ln_timing )   CALL timing_stop('p4z_nano')
-      !
-   END SUBROUTINE p4z_nano
-
-
-   SUBROUTINE p4z_diat
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE p4z_diat  ***
+      IF( ln_timing )   CALL timing_stop('p4z_mort_nano')
+      !
+   END SUBROUTINE p4z_mort_nano
+
+
+   SUBROUTINE p4z_mort_diat
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_mort_diat  ***
       !!
       !! ** Purpose :   Compute the mortality terms for diatoms
       !!
-      !! ** Method  : - ???
+      !! ** Method  : - Both quadratic and simili linear mortality terms
       !!---------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk
@@ -144,9 +150,10 @@
       !!---------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('p4z_diat')
+      IF( ln_timing )   CALL timing_start('p4z_mort_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
+      ! This is due to the production of EPS by stressed cells
       ! -------------------------------------------------------------
 
@@ -183,4 +190,7 @@
                tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) - zmortp2 * zfactsi
                tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) + zmortp2 * zfactsi
+
+               ! Half of the linear mortality term is routed to big particles
+               ! becaue of the ballasting effect
                tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zrespp2 + 0.5 * ztortp2
                tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + 0.5 * ztortp2
@@ -199,7 +209,7 @@
       ENDIF
       !
-      IF( ln_timing )   CALL timing_stop('p4z_diat')
-      !
-   END SUBROUTINE p4z_diat
+      IF( ln_timing )   CALL timing_stop('p4z_mort_diat')
+      !
+   END SUBROUTINE p4z_mort_diat
 
 
Index: NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zopt.F90
===================================================================
--- NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 13233)
@@ -80,10 +80,15 @@
       ze2(:,:,:) = 0._wp
       ze3(:,:,:) = 0._wp
-      !
-      !                                        !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
-      !                                        !  --------------------------------------------------------
+
+      ! Attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
+      ! Thus the light penetration scheme is based on a decomposition of PAR
+      ! into three wave length domains. This was first officially published
+      ! in Lengaigne et al. (2007).
+      ! --------------------------------------------------------
                      zchl3d(:,:,:) = trb(:,:,:,jpnch) + trb(:,:,:,jpdch)
       IF( ln_p5z )   zchl3d(:,:,:) = zchl3d(:,:,:)    + trb(:,:,:,jppch)
       !
+      ! Computation of the light attenuation parameters based on a 
+      ! look-up table
       DO jk = 1, jpkm1   
          DO jj = 1, jpj
@@ -99,12 +104,26 @@
          END DO
       END DO
-      !                                        !* Photosynthetically Available Radiation (PAR)
-      !                                        !  --------------------------------------
+
+      ! Photosynthetically Available Radiation (PAR)
+      ! Two cases are considered in the following : 
+      ! (1) An explicit diunal cycle is activated. In that case, mean 
+      ! QSR is used as PISCES in its current state has not been parameterized
+      ! for an explicit diurnal cycle
+      ! (2) no diurnal cycle of SW is active and in that case, QSR is used.
+      ! --------------------------------------------
       IF( l_trcdm2dc ) THEN                     !  diurnal cycle
          !
+         ! SW over the ice free zone of the grid cell. This assumes that
+         ! SW is zero below sea ice which is a very crude assumption that is 
+         ! not fully correct with LIM3 and SI3 but no information is 
+         ! currently available to do a better job. SHould be improved in the 
+         ! (near) future.
          zqsr_corr(:,:) = qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
          !
          CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 ) 
          !
+         ! Used PAR is computed for each phytoplankton species
+         ! etot_ndcy is PAR at level jk averaged over 24h.
+         ! Due to their size, they have different light absorption characteristics
          DO jk = 1, nksrp      
             etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
@@ -117,27 +136,40 @@
             END DO
          ENDIF
-         !
+
+         ! SW over the ice free zone of the grid cell. This assumes that
+         ! SW is zero below sea ice which is a very crude assumption that is 
+         ! not fully correct with LIM3 and SI3 but no information is 
+         ! currently available to do a better job. SHould be improved in the 
+         ! (near) future.
          zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
          !
          CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) 
-         !
+
+         ! Total PAR computation at level jk that includes the diurnal cycle
          DO jk = 1, nksrp      
             etot(:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
          END DO
          !
-      ELSE
-         !
+      ELSE ! no diurnal cycle
+         !
+         ! SW over the ice free zone of the grid cell. This assumes that
+         ! SW is zero below sea ice which is a very crude assumption that is 
+         ! not fully correct with LIM3 and SI3 but no information is 
+         ! currently available to do a better job. SHould be improved in the 
+         ! (near) future.
          zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
-         !
+
          CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100  ) 
-         !
+
+         ! Used PAR is computed for each phytoplankton species
+         ! Due to their size, they have different light absorption characteristics
          DO jk = 1, nksrp      
-            etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
-            enano(:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
-            ediat(:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)
+            etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk) ! Total PAR
+            enano(:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk) ! Nanophytoplankton
+            ediat(:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk) ! Diatoms
          END DO
          IF( ln_p5z ) THEN
             DO jk = 1, nksrp      
-              epico(:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
+              epico(:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk) ! Picophytoplankton (PISCES-QUOTA)
             END DO
          ENDIF
@@ -146,4 +178,5 @@
 
 
+      ! Biophysical feedback part (computation of vertical penetration of SW)
       IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics)
          !                                     !  ------------------------
@@ -154,8 +187,12 @@
             etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
          END DO
-         !                                     !  ------------------------
-      ENDIF
-      !                                        !* Euphotic depth and level
-      neln   (:,:) = 1                            !  ------------------------
+
+      ENDIF
+
+      ! Euphotic depth and level
+      ! Two definitions of the euphotic zone are used here. 
+      ! (1) The classical definition based on the relative threshold value
+      ! (2) An alternative definition based on a absolute threshold value.
+      ! -------------------------------------------------------------------
       heup   (:,:) = gdepw_n(:,:,2)
       heup_01(:,:) = gdepw_n(:,:,2)
@@ -176,8 +213,11 @@
       END DO
       !
+      ! The euphotic depth can not exceed 300 meters.
       heup   (:,:) = MIN( 300., heup   (:,:) )
       heup_01(:,:) = MIN( 300., heup_01(:,:) )
-      !                                        !* mean light over the mixed layer
-      zdepmoy(:,:)   = 0.e0                    !  -------------------------------
+
+      ! Mean PAR over the mixed layer
+      ! -----------------------------
+      zdepmoy(:,:)   = 0.e0
       zetmp1 (:,:)   = 0.e0
       zetmp2 (:,:)   = 0.e0
@@ -187,6 +227,6 @@
             DO ji = 1, jpi
                IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t_n(ji,jj,jk) ! remineralisation
-                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t_n(ji,jj,jk) ! production
+                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t_n(ji,jj,jk) ! Actual PAR
+                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t_n(ji,jj,jk) ! Par averaged over 24h
                   zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t_n(ji,jj,jk)
                ENDIF
@@ -195,5 +235,5 @@
       END DO
       !
-      emoy(:,:,:) = etot(:,:,:)       ! remineralisation
+      emoy(:,:,:) = etot(:,:,:)       ! PAR
       zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle 
       !
@@ -209,5 +249,7 @@
          END DO
       END DO
-      !
+
+      ! Computation of the mean usable light for the different phytoplankton
+      ! groups based on their absorption characteristics.
       zdepmoy(:,:)   = 0.e0
       zetmp3 (:,:)   = 0.e0
@@ -218,6 +260,6 @@
             DO ji = 1, jpi
                IF( gdepw_n(ji,jj,jk+1) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
-                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t_n(ji,jj,jk) ! production
-                  zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * e3t_n(ji,jj,jk) ! production
+                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t_n(ji,jj,jk) ! Nanophytoplankton
+                  zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * e3t_n(ji,jj,jk) ! Diatoms 
                   zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t_n(ji,jj,jk)
                ENDIF
@@ -241,4 +283,5 @@
       !
       IF( ln_p5z ) THEN
+         ! Picophytoplankton when using PISCES-QUOTA
          zetmp5 (:,:) = 0.e0
          DO jk = 1, nksrp
@@ -246,5 +289,5 @@
                DO ji = 1, jpi
                   IF( gdepw_n(ji,jj,jk+1) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
-                     zetmp5(ji,jj)  = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t_n(ji,jj,jk) ! production
+                     zetmp5(ji,jj)  = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t_n(ji,jj,jk)
                   ENDIF
                END DO
@@ -265,4 +308,6 @@
          END DO
       ENDIF
+
+      ! Output of the diagnostics
       IF( lk_iomput ) THEN
         IF( knt == nrdttrc ) THEN
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 13233)
@@ -89,5 +89,5 @@
       solgoc = 0.04/ 2.56 * 1./ ( 1.-50**(-0.04) )
 
-      ! Initialisation of temporary arrys
+      ! Initialisation of temporary arrays
       IF( ln_p4z ) THEN
          zremipoc(:,:,:) = xremip
@@ -419,5 +419,5 @@
               DO ji = 1, jpi
                  IF (tmask(ji,jj,jk) == 1.) THEN
-                    ! POC disaggregation by bacterial activity. It is a function
+                    ! POC solubilisation by bacterial activity. It is a function
                     ! of the mean lability and of temperature 
                     ! ----------------------------------------------------------
@@ -538,5 +538,5 @@
       ALLOCATE( alphan(jcpoc) , reminp(jcpoc) , alphap(jpi,jpj,jpk,jcpoc) )
       !
-      IF (jcpoc > 1) THEN
+      IF (jcpoc > 1) THEN  ! Case when more than one lability class is used
          !
          remindelta = LOG(4. * 1000. ) / REAL(jcpoc-1, wp)
@@ -562,5 +562,5 @@
          reminp(jcpoc) = reminp(jcpoc) * xremip / alphan(jcpoc)
 
-      ELSE
+      ELSE  ! Only one lability class is used
          alphan(jcpoc) = 1.
          reminp(jcpoc) = xremip
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 13233)
@@ -26,16 +26,16 @@
    PUBLIC   p4z_prod_alloc   ! called in trcini_pisces.F90
 
-   REAL(wp), PUBLIC ::   pislopen     !:
-   REAL(wp), PUBLIC ::   pisloped     !:
-   REAL(wp), PUBLIC ::   xadap        !:
-   REAL(wp), PUBLIC ::   excretn      !:
-   REAL(wp), PUBLIC ::   excretd      !:
-   REAL(wp), PUBLIC ::   bresp        !:
-   REAL(wp), PUBLIC ::   chlcnm       !:
-   REAL(wp), PUBLIC ::   chlcdm       !:
-   REAL(wp), PUBLIC ::   chlcmin      !:
-   REAL(wp), PUBLIC ::   fecnm        !:
-   REAL(wp), PUBLIC ::   fecdm        !:
-   REAL(wp), PUBLIC ::   grosip       !:
+   REAL(wp), PUBLIC ::   pislopen     !:  P-I slope of nanophytoplankton
+   REAL(wp), PUBLIC ::   pisloped     !:  P-I slope of diatoms
+   REAL(wp), PUBLIC ::   xadap        !:  Adaptation factor to low light 
+   REAL(wp), PUBLIC ::   excretn      !:  Excretion ratio of nanophyto
+   REAL(wp), PUBLIC ::   excretd      !:  Excretion ratio of diatoms
+   REAL(wp), PUBLIC ::   bresp        !:  Basal respiration rate
+   REAL(wp), PUBLIC ::   chlcnm       !:  Maximum Chl/C ratio of nano
+   REAL(wp), PUBLIC ::   chlcdm       !:  Maximum Chl/C ratio of diatoms
+   REAL(wp), PUBLIC ::   chlcmin      !:  Minimum Chl/C ratio of phytoplankton
+   REAL(wp), PUBLIC ::   fecnm        !:  Maximum Fe/C ratio of nano
+   REAL(wp), PUBLIC ::   fecdm        !:  Maximum Fe/C ratio of diatoms
+   REAL(wp), PUBLIC ::   grosip       !:  Mean Si/C ratio of diatoms
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotan   !: proxy of N quota in Nanophyto
@@ -61,4 +61,5 @@
       !!                Computes also the uptake of Iron and Si as well 
       !!                as the chlorophyll content of the cells
+      !!                PISCES relies on a mixed Monod-Quota formalism 
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, knt   !
@@ -94,5 +95,6 @@
       zmxl_fac(:,:,:) = 0._wp ; zmxl_chl(:,:,:) = 0._wp 
 
-      ! Computation of the maximimum production
+      ! Computation of the maximimum production. Based on a Q10 description
+      ! of the thermal dependency
       ! Parameters are taken from Bissinger et al. (2008)
       zprmaxn(:,:,:) = 0.8_wp * r1_rday * tgfunc(:,:,:)
@@ -100,4 +102,5 @@
 
       ! compute the day length depending on latitude and the day
+      ! Astronomical parameterization taken from HAMOCC3
       zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
       zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
@@ -115,7 +118,8 @@
       ! 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)
+      ! day length (Shatwell et al., 2012). The correcting factor is zmxl_fac. 
+      ! zmxl_chl is the fractional day length and is used to compute the mean
+      ! PAR during daytime. The effect of mixing is computed using the 
+      ! absolute light level definition of the euphotic zone
       ! ------------------------------------------------------------------------- 
       DO jk = 1, jpkm1
@@ -137,6 +141,4 @@
       zprdia(:,:,:) = zprmaxd(:,:,:) * zmxl_fac(:,:,:)
 
-      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 
@@ -155,7 +157,12 @@
                   ! The initial slope of the PI curve can be increased for nano
                   ! to account for photadaptation, for instance in the DCM
+                  ! This parameterization is adhoc and should be either 
+                  ! improved or removed in future versions of the model
+
+                  ! Nanophytoplankton
                   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)
-                  !
+
+                  ! Diatoms
                   zpislopeadd(ji,jj,jk) = (pislopen * zconctemp2 + pisloped * zconctemp) / ( trb(ji,jj,jk,jpdia) + rtrn )   &
                   &                   * trb(ji,jj,jk,jpdch) /( trb(ji,jj,jk,jpdia) * 12. + rtrn)
@@ -170,5 +177,6 @@
                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                    ! Computation of production function for Carbon
-                   !  ---------------------------------------------
+                   ! Actual light levels are used here 
+                   ! ----------------------------------------------
                    zpislopen = zpislopeadn(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
                    &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
@@ -177,6 +185,10 @@
                    zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
                    zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
+
                    !  Computation of production function for Chlorophyll
-                   !--------------------------------------------------
+                   !  Mean light level in the mixed layer (when appropriate)
+                   !  is used here (acclimation is in general slower than 
+                   !  the characteristic time scales of vertical mixing)
+                   !  ------------------------------------------------------
                    zpislopen = zpislopeadn(ji,jj,jk) / ( zprmaxn(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
                    zpisloped = zpislopeadd(ji,jj,jk) / ( zprmaxd(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
@@ -188,7 +200,7 @@
       END DO
 
-      !  Computation of a proxy of the N/C ratio
-      !  Steady state is assumed
-      !  ---------------------------------------
+      !  Computation of a proxy of the N/C quota from nutrient limitation 
+      !  and light limitation. Steady state is assumed to allow the computation
+      !  ----------------------------------------------------------------------
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -210,4 +222,5 @@
 
                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+
                    ! Si/C of diatoms
                    ! ------------------------
@@ -215,4 +228,7 @@
                    ! Si/C is arbitrariliy increased for very high Si concentrations
                    ! to mimic the very high ratios observed in the Southern Ocean (zsilfac2)
+                   ! A parameterization derived from Flynn (2003) is used for the control
+                   ! when Si is not limiting which is similar to the parameterisation
+                   ! proposed by Gurney and Davidson (1999).
                    ! -----------------------------------------------------------------------
                  zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
@@ -271,8 +287,8 @@
                   sizena(ji,jj,jk) = min(xsizern, max( sizena(ji,jj,jk), zsizetmp ) )
 
-                  !  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 
+                  ! Iron uptake rates of nanophytoplankton. Upregulation is  
+                  ! not parameterized at low iron concentrations as observations
+                  ! do not suggest it for accimated cells. Uptake is
+                  ! downregulated when the quota is close to the maximum quota
                   zratio = 1.0 - MIN(1.0,trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * fecnm + rtrn ) )
                   zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) ) 
@@ -281,5 +297,5 @@
                   &          + xnanonh4(ji,jj,jk) ) * (1. - xnanofer(ji,jj,jk) ) )   &
                   &          * xnanofer(ji,jj,jk) * zmax * trb(ji,jj,jk,jpphy) * rfact2
-                  !  production terms of 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
 
@@ -297,8 +313,8 @@
                   sizeda(ji,jj,jk) = min(xsizerd, max( sizeda(ji,jj,jk), zsizetmp ) )
 
-                  !  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 
+                  ! Iron uptake rates of diatoms. Upregulation is  
+                  ! not parameterized at low iron concentrations as observations
+                  ! do not suggest it for accimated cells. Uptake is
+                  ! downregulated when the quota is close to the maximum quota
                   zratio = 1.0 - MIN(1.0, trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * fecdm + rtrn ) )
                   zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) ) 
@@ -323,4 +339,5 @@
                   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)
@@ -328,8 +345,10 @@
                   zprochln = zprochln + (chlcnm_n-chlcmin) * 12. * zprod / &
                                         & (  zpislopeadn(ji,jj,jk) * znanotot +rtrn)
-                  !  production terms of 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)
@@ -337,5 +356,6 @@
                   zprochld = zprochld + (chlcdm_n-chlcmin) * 12. * zprod / &
                                         & ( zpislopeadd(ji,jj,jk) * zdiattot +rtrn )
-                  !   Update the arrays TRA which contain the Chla sources and sinks
+
+                  ! Update the arrays TRA which contain the Chla sources and sinks
                   tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) + zprochln * texcretn
                   tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) + zprochld * texcretd
@@ -377,5 +397,5 @@
         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 
@@ -403,4 +423,5 @@
 
 
+    ! Output of the diagnostics
     ! Total primary production per year
     IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
@@ -423,5 +444,5 @@
               CALL iom_put( "PPNEWN"  , 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
@@ -434,5 +455,5 @@
               CALL iom_put( "PFeN"  , zw3d )
               !
-              zw3d(:,:,:) = zprofed(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron production by  diatomes
+              zw3d(:,:,:) = zprofed(:,:,:) * zfact * tmask(:,:,:)  ! biogenic iron production by  diatoms
               CALL iom_put( "PFeD"  , zw3d )
           ENDIF
@@ -546,4 +567,5 @@
       INTEGER ::   ios   ! Local integer
       !
+      ! Namelist block
       NAMELIST/namp4zprod/ pislopen, pisloped, xadap, bresp, excretn, excretd,  &
          &                 chlcnm, chlcdm, chlcmin, fecnm, fecdm, grosip
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zrem.F90	(revision 13233)
@@ -92,6 +92,8 @@
       ! Computation of the mean bacterial concentration
       ! this parameterization has been deduced from a model version
-      ! that was modeling explicitely bacteria
-      ! -------------------------------------------------------
+      ! that was modeling explicitely bacteria. This is a very old param 
+      ! that will be very soon updated based on results from a much more
+      ! recent version of PISCES with bacteria.
+      ! ----------------------------------------------------------------
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -111,5 +113,5 @@
       END DO
 
-      IF( ln_p4z ) THEN
+      IF( ln_p4z ) THEN ! Standard PISCES code
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -152,5 +154,5 @@
             END DO
          END DO
-      ELSE
+      ELSE ! PISCES-QUOTA
          DO jk = 1, jpkm1
             DO jj = 1, jpj
@@ -231,5 +233,5 @@
 
                ! Bacterial uptake of iron. No iron is available in DOC. So
-               ! Bacteries are obliged to take up iron from the water. Some
+               ! Bacteria are obliged to take up iron from the water. Some
                ! studies (especially at Papa) have shown this uptake to be significant
                ! ---------------------------------------------------------------------
@@ -262,5 +264,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ! Remineralization rate of BSi depedant on T and saturation
+
+               ! Remineralization rate of BSi dependent on T and saturation
                ! The parameterization is taken from Ridgwell et al. (2002) 
                ! ---------------------------------------------------------
@@ -269,6 +272,12 @@
                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
+
+               ! Two fractions of bSi are considered : a labile one and a more
+               ! refractory one based on the commonly observed two step 
+               ! dissolution of bSi (initial rapid dissolution followed by 
+               ! more slowly dissolution).
                ! 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 )  &
@@ -280,4 +289,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
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 13233)
@@ -35,5 +35,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: sdenit     !: Nitrate reduction in the sediments
    REAL(wp) :: r1_rday                  !: inverse of rday
-   LOGICAL, SAVE :: lk_sed
+   LOGICAL, SAVE :: lk_sed              !: Explicit sediment module
 
    !!----------------------------------------------------------------------
@@ -48,5 +48,5 @@
       !!                     ***  ROUTINE p4z_sed  ***
       !!
-      !! ** Purpose :   Compute loss of biogenic matter in the sediments. This
+      !! ** Purpose : Compute the loss of biogenic matter in the sediments. This
       !!              is by no way a real sediment model. The loss is simply 
       !!              computed from metamodels.
@@ -55,5 +55,6 @@
       !!              N2 fixation is modeled using an implicit approach
       !!
-      !! ** Method  : - ???
+      !! ** Method  : - Fluxes with the sediments are mainly modeled using
+      !!                statiscal metamodels.
       !!---------------------------------------------------------------------
       !
@@ -81,4 +82,5 @@
       IF( kt == nittrc000 .AND. knt == 1 )   THEN
           r1_rday  = 1. / rday
+          ! Configuration with an active two-way sediment module 
           IF (ln_sediment .AND. ln_sed_2way) THEN
              lk_sed = .TRUE.
@@ -109,6 +111,9 @@
          !                                              
          ALLOCATE( zironice(jpi,jpj) )
-         !                                              
+
          ! Compute the iron flux between sea ice and sea water
+         ! Simple parameterization assuming a fixed constant concentration in
+         ! sea-ice (icefeinput)
+         ! ------------------------------------------------------------------
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -133,8 +138,10 @@
          !                                              
          ALLOCATE( zsidep(jpi,jpj), zpdep(jpi,jpj,jpk), zirondep(jpi,jpj,jpk) )
+
          ! 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)
+         ! mfrac is the mean iron relative weight content of dust
          ! ------------------------------------------------------------------
          IF( ln_solub ) THEN
@@ -143,4 +150,5 @@
             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
@@ -148,4 +156,5 @@
          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/d
@@ -159,7 +168,7 @@
          DO jk = 2, jpkm1
             zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -gdept_n(:,:,jk) / (250. * wdust) )
-!            zpdep   (:,:,jk) = zirondep(:,:,jk) * 0.023
             zpdep   (:,:,jk) = zirondep(:,:,jk) * 0.38 / po4r
          END DO
+
          ! Solubilization of particles in the water column (Si, P, Fe)
          tra(:,:,1,jpsil) = tra(:,:,1,jpsil) + zsidep  (:,:)
@@ -197,4 +206,5 @@
             ENDDO
          ENDDO
+
          ! When prognostic ligands are activated, ligands are supplied 
          ! to the ocean by rivers. We assume that the amount of ligands
@@ -262,5 +272,5 @@
          ! ------------------------------------------------------
          IF( ln_ironsed ) THEN
-                            tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2
+            tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2
             !
             IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironsed" ) )   &
@@ -434,4 +444,6 @@
                   xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
                   zlim = ( 1.- xdiano3 - xdianh4 )
+                  ! Nitrogen fixation is almost fully halted when the N 
+                  ! limitation term (xdiano3+xdianh4) is > 0.9
                   IF( zlim <= 0.1 )   zlim = 0.01
                   zfact = zlim * rfact2
@@ -454,4 +466,7 @@
                   xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
                   zlim = ( 1.- xdiano3 - xdianh4 )
+
+                  ! Nitrogen fixation is almost fully halted when the N 
+                  ! limitation term (xdiano3+xdianh4) is > 0.9
                   IF( zlim <= 0.1 )   zlim = 0.01
                   zfact = zlim * rfact2
@@ -474,4 +489,7 @@
                DO ji = 1, jpi
                   zfact = nitrpot(ji,jj,jk) * nitrfix
+                  ! 1/3 of the diazotrophs growth is supposed to be excreted
+                  ! as NH4. 1/3 as DOC and the rest is routed POC and GOC as 
+                  ! a result of mortality by predation. Completely adhoc param 
                   tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0
                   tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0
@@ -481,4 +499,5 @@
                   tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0
                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
+                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C
                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0
                   tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
@@ -494,6 +513,10 @@
                DO ji = 1, jpi
                   zfact = nitrpot(ji,jj,jk) * nitrfix
+                  ! 1/3 of the diazotrophs growth is supposed to be excreted
+                  ! as NH4. 1/3 as DOC and the rest is routed POC and GOC as 
+                  ! a result of mortality by predation. Completely adhoc param 
                   tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0
                   tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0
+                  ! N/P ratio of diazotrophs is supposed to be 46
                   tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) &
                   &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
@@ -510,4 +533,5 @@
                   tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0
                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
+                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C
                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0 
                   tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 13233)
@@ -95,5 +95,5 @@
       END DO
 
-      ! Sinking speed of the small particles is constant
+      ! Sinking speed of the small particles is always constant
       wsbio3(:,:,:) = wsbio
 
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 13233)
@@ -74,5 +74,5 @@
         IF( .NOT. ln_rsttr ) THEN
             CALL p4z_che            ! initialize the chemical constants
-            CALL ahini_for_at(hi)   !  set PH at kt=nit000
+            CALL ahini_for_at(hi)   ! set PH at kt=nit000
             t_oce_co2_flx_cum = 0._wp
         ELSE
@@ -84,11 +84,11 @@
       IF( ln_pisdmp .AND. MOD( kt - nn_dttrc, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt )      ! Relaxation of some tracers
       !
-      rfact = r2dttrc
+      rfact = r2dttrc  ! time step of PISCES
       !
       IF( ( ln_top_euler .AND. kt == nittrc000 )  .OR. ( .NOT.ln_top_euler .AND. kt <= nittrc000 + nn_dttrc ) ) THEN
-         rfactr  = 1. / rfact
-         rfact2  = rfact / REAL( nrdttrc, wp )
-         rfact2r = 1. / rfact2
-         xstep = rfact2 / rday         ! Time step duration for biology
+         rfactr  = 1. / rfact  ! inverse of the time step
+         rfact2  = rfact / REAL( nrdttrc, wp )  ! time step of the biological SMS
+         rfact2r = 1. / rfact2  ! Inverse of the biological time step
+         xstep = rfact2 / rday         ! Time step duration for biology relative to a day
          IF(lwp) WRITE(numout,*) 
          IF(lwp) WRITE(numout,*) '    Passive Tracer  time step    rfact  = ', rfact, ' rdt = ', rdt
@@ -147,5 +147,5 @@
             tra(:,:,:,jn) = 0._wp
          END DO
-         !
+         ! Euler-forward temporal scheme
          IF( ln_top_euler ) THEN
             DO jn = jp_pcs0, jp_pcs1
@@ -167,5 +167,5 @@
          !
          CALL sed_model( kt )     !  Main program of Sediment model
-         !
+         ! Eulor forward temporal scheme
          IF( ln_top_euler ) THEN
             DO jn = jp_pcs0, jp_pcs1
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 13233)
@@ -59,16 +59,16 @@
 
    !!*  Allometric variations of the quotas
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnnmin    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnnmax    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpnmin    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpnmax    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnpmin    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnpmax    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqppmin    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqppmax    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqndmin    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqndmax    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpdmin    !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpdmax    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnnmin    !: Minimum N quota of nanophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnnmax    !: Maximum N quota of nanophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpnmin    !: Minimum P quota of nanophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpnmax    !: Maximum P quota of picophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnpmin    !: Minimum N quota of picophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnpmax    !: Maximum N quota of picophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqppmin    !: Minimum P quota of picophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqppmax    !: Maximum P quota of picophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqndmin    !: Minimum N quota of diatoms
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqndmax    !: Maximum N quota of diatoms
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpdmin    !: Minimum P quota of diatoms
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpdmax    !: Maximum P quota of diatoms
 
    !!* Phytoplankton nutrient limitation terms
@@ -521,8 +521,8 @@
       !! ** Purpose :   Initialization of nutrient limitation parameters
       !!
-      !! ** Method  :   Read the nampislim and nampisquota namelists and check
+      !! ** Method  :   Read the namp5zlim and nampisquota namelists and check
       !!      the parameters called at the first timestep (nittrc000)
       !!
-      !! ** input   :   Namelist nampislim
+      !! ** input   :   Namelist namp5zlim
       !!
       !!----------------------------------------------------------------------
@@ -540,11 +540,11 @@
       !!----------------------------------------------------------------------
       !
-      REWIND( numnatp_ref )              ! Namelist nampislim in reference namelist : Pisces nutrient limitation parameters
+      REWIND( numnatp_ref )              ! Namelist namp5zlim in reference namelist : Pisces nutrient limitation parameters
       READ  ( numnatp_ref, namp5zlim, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampislim in reference namelist' )
-      !
-      REWIND( numnatp_cfg )              ! Namelist nampislim in configuration namelist : Pisces nutrient limitation parameters 
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp5zlim in reference namelist' )
+      !
+      REWIND( numnatp_cfg )              ! Namelist namp5zlim in configuration namelist : Pisces nutrient limitation parameters 
       READ  ( numnatp_cfg, namp5zlim, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nampislim in configuration namelist' )
+902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp5zlim in configuration namelist' )
       IF(lwm) WRITE ( numonp, namp5zlim )
       !
@@ -614,4 +614,5 @@
       ENDIF
       !
+      ! Metabolic cost of nitrate and ammonium utilisation
       zpsino3  = 2.3 * rno3
       zpsinh4  = 1.8 * rno3
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 13233)
@@ -52,4 +52,6 @@
    REAL(wp), PUBLIC ::  grazflux     !: mesozoo flux feeding rate
    REAL(wp), PUBLIC ::  xfracmig     !: Fractional biomass of meso that performs DVM
+   REAL(wp), PUBLIC ::  xsigma2      !: Width of the predation window
+   REAL(wp), PUBLIC ::  xsigma2del   !: Maximum width of the predation window at low food density
    LOGICAL,  PUBLIC ::  bmetexc2     !: Use of excess carbon for respiration
    LOGICAL , PUBLIC ::  ln_dvm_meso  !: Boolean to activate DVM of mesozooplankton
@@ -189,5 +191,5 @@
                ! ----------------------------------------------------------
                zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
-               zsigma = 0.5 + 1.0 * zsigma
+               zsigma = xsigma2 + xsigma2del * zsigma
                ! Nanophytoplankton and diatoms are the only preys considered
                ! to be close enough to have potential interference
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmicro.F90	(revision 13233)
@@ -15,5 +15,5 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p4zlim
+   USE p4zlim          !  PISCES nutrient limitation term of PISCES std
    USE p5zlim          !  Phytoplankton limitation terms
    USE iom             !  I/O manager
@@ -50,4 +50,6 @@
    REAL(wp), PUBLIC ::  srespir     !: half sturation constant for grazing 1 
    REAL(wp), PUBLIC ::  ssigma      !: Fraction excreted as semi-labile DOM
+   REAL(wp), PUBLIC ::  xsigma      !: Width of the grazing window
+   REAL(wp), PUBLIC ::  xsigmadel   !: Maximum additional width of the grazing window at low food density
    LOGICAL,  PUBLIC ::  bmetexc     !: Use of excess carbon for respiration
 
@@ -161,5 +163,5 @@
                ! ----------------------------------------------------------
                zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
-               zsigma = 0.5 + 1.0*zsigma
+               zsigma = xsigma + xsigmadel * zsigma
                zdiffpn = exp( -ABS(log(0.5 * sizep(ji,jj,jk) / (3.0 * sizen(ji,jj,jk) + rtrn )) )**2 / zsigma**2 )
                zdiffdn = exp( -ABS(log(3.0 * sizen(ji,jj,jk) / (5.0 * sized(ji,jj,jk) + rtrn )) )**2 / zsigma**2)
@@ -179,17 +181,25 @@
                !   Microzooplankton regular grazing on the different preys
                !   -------------------------------------------------------
-               zgraznc   = zgraze  * ztmp1  * zdenom
+               !   Nanophytoplankton
+               zgraznc   = zgraze  * ztmp1  * zdenom  
                zgraznn   = zgraznc * trb(ji,jj,jk,jpnph) / (trb(ji,jj,jk,jpphy) + rtrn)
                zgraznp   = zgraznc * trb(ji,jj,jk,jppph) / (trb(ji,jj,jk,jpphy) + rtrn)
                zgraznf   = zgraznc * trb(ji,jj,jk,jpnfe) / (trb(ji,jj,jk,jpphy) + rtrn)
+
+               ! Picophytoplankton
                zgrazpc   = zgraze  * ztmp2  * zdenom
                zgrazpn   = zgrazpc * trb(ji,jj,jk,jpnpi) / (trb(ji,jj,jk,jppic) + rtrn)
                zgrazpp   = zgrazpc * trb(ji,jj,jk,jpppi) / (trb(ji,jj,jk,jppic) + rtrn)
                zgrazpf   = zgrazpc * trb(ji,jj,jk,jppfe) / (trb(ji,jj,jk,jppic) + rtrn)
+               ! Microzooplankton
                zgrazz    = zgraze  * ztmp5   * zdenom
+
+               ! small POC
                zgrazpoc  = zgraze  * ztmp3   * zdenom
                zgrazpon  = zgrazpoc * trb(ji,jj,jk,jppon) / ( trb(ji,jj,jk,jppoc) + rtrn )
                zgrazpop  = zgrazpoc * trb(ji,jj,jk,jppop) / ( trb(ji,jj,jk,jppoc) + rtrn )
                zgrazpof  = zgrazpoc* trb(ji,jj,jk,jpsfe) / (trb(ji,jj,jk,jppoc) + rtrn)
+
+               ! Diatoms
                zgrazdc   = zgraze  * ztmp4  * zdenom
                zgrazdn   = zgrazdc * trb(ji,jj,jk,jpndi) / (trb(ji,jj,jk,jpdia) + rtrn)
@@ -197,4 +207,5 @@
                zgrazdf   = zgrazdc * trb(ji,jj,jk,jpdfe) / (trb(ji,jj,jk,jpdia) + rtrn)
                !
+               ! Total ingestion rates in C, P, Fe, N
                zgraztotc = zgraznc + zgrazpoc + zgrazdc + zgrazz + zgrazpc
                zgraztotn = zgraznn + zgrazpn + zgrazpon + zgrazdn + zgrazz * no3rat3
@@ -221,10 +232,14 @@
                zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / ferat3)
                zbeta     = MAX( 0., (epsher - epshermin) )
+               ! Food density deprivation of GGE
                zepsherf  = epshermin + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
+               ! Food quality deprivation of GGE
                zepsherq  = 0.5 + (1.0 - 0.5) * zepshert * ( 1.0 + 1.0 ) / ( zepshert + 1.0 )
+               ! Actual GGE
                zepsherv  = zepsherf * zepshert * zepsherq
 
                ! Respiration of microzooplankton
                ! Excess carbon in the food is used preferentially
+               ! when activated by zmetexcess
                ! ------------------------------------------------
                zexcess  = zgraztotc * zepsherf * (1.0 - zepshert) * zmetexcess
@@ -380,5 +395,6 @@
          &                xprefp, xprefd, xprefz, xthreshdia, xthreshphy, &
          &                xthreshpic, xthreshpoc, xthreshzoo, xthresh, xkgraz, &
-         &                epsher, epshermin, ssigma, srespir, unassc, unassn, unassp
+         &                epsher, epshermin, ssigma, srespir, unassc, unassn, unassp,   &
+         &                xsigma, xsigmadel   
       !!----------------------------------------------------------------------
       !
@@ -420,4 +436,6 @@
          WRITE(numout,*) '    half sturation constant for grazing 1           xkgraz      =', xkgraz
          WRITE(numout,*) '    Use of excess carbon for respiration            bmetexc     =', bmetexc
+         WRITE(numout,*) '      Width of the grazing window                     xsigma      =', xsigma
+         WRITE(numout,*) '      Maximum additional width of the grazing window  xsigmadel   =', xsigmadel
       ENDIF
       !
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmort.F90	(revision 13233)
@@ -53,18 +53,18 @@
       !!---------------------------------------------------------------------
 
-      CALL p5z_nano            ! nanophytoplankton
-      CALL p5z_pico            ! picophytoplankton
-      CALL p5z_diat            ! diatoms
+      CALL p5z_mort_nano            ! nanophytoplankton
+      CALL p5z_mort_pico            ! picophytoplankton
+      CALL p5z_mort_diat            ! diatoms
 
    END SUBROUTINE p5z_mort
 
 
-   SUBROUTINE p5z_nano
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE p5z_nano  ***
+   SUBROUTINE p5z_mort_nano
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p5z_mort_nano  ***
       !!
       !! ** Purpose :   Compute the mortality terms for nanophytoplankton
       !!
-      !! ** Method  : - ???
+      !! ** Method  : - Both quadratic and simili linear mortality terms
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jk
@@ -75,5 +75,5 @@
       !!---------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('p5z_nano')
+      IF( ln_timing )   CALL timing_start('p5z_mort_nano')
       !
       prodcal(:,:,:) = 0.  !: calcite production variable set to zero
@@ -82,4 +82,5 @@
             DO ji = 1, jpi
                zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e-9 ), 0.e0 )
+
                ! Quadratic mortality of nano due to aggregation during
                ! blooms (Doney et al. 1996)
@@ -127,16 +128,16 @@
        ENDIF
       !
-      IF( ln_timing )   CALL timing_stop('p5z_nano')
-      !
-   END SUBROUTINE p5z_nano
-
-
-   SUBROUTINE p5z_pico
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE p5z_pico  ***
+      IF( ln_timing )   CALL timing_stop('p5z_mort_nano')
+      !
+   END SUBROUTINE p5z_mort_nano
+
+
+   SUBROUTINE p5z_mort_pico
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p5z_mort_pico  ***
       !!
       !! ** Purpose :   Compute the mortality terms for picophytoplankton
       !!
-      !! ** Method  : - ???
+      !! ** Method  : - Both quadratic and semilininear terms are used
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jk
@@ -147,5 +148,5 @@
       !!---------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('p5z_pico')
+      IF( ln_timing )   CALL timing_start('p5z_mort_pico')
       !
       DO jk = 1, jpkm1
@@ -191,12 +192,12 @@
        ENDIF
       !
-      IF( ln_timing )   CALL timing_stop('p5z_pico')
-      !
-   END SUBROUTINE p5z_pico
-
-
-   SUBROUTINE p5z_diat
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE p5z_diat  ***
+      IF( ln_timing )   CALL timing_stop('p5z_mort_pico')
+      !
+   END SUBROUTINE p5z_mort_pico
+
+
+   SUBROUTINE p5z_mort_diat
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p5z_mort_diat  ***
       !!
       !! ** Purpose :   Compute the mortality terms for diatoms
@@ -211,5 +212,5 @@
       !!---------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('p5z_diat')
+      IF( ln_timing )   CALL timing_start('p5z_mort_diat')
       !
 
@@ -271,7 +272,7 @@
       ENDIF
       !
-      IF( ln_timing )   CALL timing_stop('p5z_diat')
-      !
-   END SUBROUTINE p5z_diat
+      IF( ln_timing )   CALL timing_stop('p5z_mort_diat')
+      !
+   END SUBROUTINE p5z_mort_diat
 
 
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zprod.F90	(revision 13233)
@@ -2,5 +2,6 @@
    !!======================================================================
    !!                         ***  MODULE p5zprod  ***
-   !! TOP :  Growth Rate of the two phytoplanktons groups 
+   !! TOP :  Growth Rate of the three phytoplanktons groups 
+   !!        PISCES-QUOTA version of the module
    !!======================================================================
    !! History :   1.0  !  2004     (O. Aumont) Original code
@@ -29,24 +30,24 @@
 
    !! * Shared module variables
-   REAL(wp), PUBLIC ::  pislopen        !:
-   REAL(wp), PUBLIC ::  pislopep        !:
-   REAL(wp), PUBLIC ::  pisloped        !:
-   REAL(wp), PUBLIC ::  xadap           !:
-   REAL(wp), PUBLIC ::  excretn         !:
-   REAL(wp), PUBLIC ::  excretp         !:
-   REAL(wp), PUBLIC ::  excretd         !:
-   REAL(wp), PUBLIC ::  bresp           !:
-   REAL(wp), PUBLIC ::  thetanpm        !:
-   REAL(wp), PUBLIC ::  thetannm        !:
-   REAL(wp), PUBLIC ::  thetandm        !:
-   REAL(wp), PUBLIC ::  chlcmin         !:
-   REAL(wp), PUBLIC ::  grosip          !:
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   zdaylen
+   REAL(wp), PUBLIC ::  pislopen        !: P-I slope of nanophytoplankton
+   REAL(wp), PUBLIC ::  pislopep        !: P-I slope of picophytoplankton
+   REAL(wp), PUBLIC ::  pisloped        !: P-I slope of diatoms
+   REAL(wp), PUBLIC ::  xadap           !: Adaptation factor to low light
+   REAL(wp), PUBLIC ::  excretn         !: Excretion ratio of nanophyto
+   REAL(wp), PUBLIC ::  excretp         !: Excretion ratio of picophyto
+   REAL(wp), PUBLIC ::  excretd         !: Excretion ratio of diatoms
+   REAL(wp), PUBLIC ::  bresp           !: Basal respiration rate
+   REAL(wp), PUBLIC ::  thetanpm        !: Maximum Chl/N ratio of picophyto
+   REAL(wp), PUBLIC ::  thetannm        !: Maximum Chl/N ratio of nanophyto
+   REAL(wp), PUBLIC ::  thetandm        !: Maximum Chl/N ratio of diatoms
+   REAL(wp), PUBLIC ::  chlcmin         !: Minimum Chl/C ratio of phytoplankton
+   REAL(wp), PUBLIC ::  grosip          !: Mean Si/C ratio of diatoms
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   zdaylen ! day length
    
    REAL(wp) :: r1_rday                !: 1 / rday
-   REAL(wp) :: texcretn               !: 1 - excret 
+   REAL(wp) :: texcretn               !: 1 - excretn 
    REAL(wp) :: texcretp               !: 1 - excretp 
-   REAL(wp) :: texcretd               !: 1 - excret2        
+   REAL(wp) :: texcretd               !: 1 - excretd        
 
    !!----------------------------------------------------------------------
@@ -63,6 +64,6 @@
       !! ** Purpose :   Compute the phytoplankton production depending on
       !!              light, temperature and nutrient availability
-      !!
-      !! ** Method  : - ???
+      !!              Computes also the uptake of nutrients. PISCES-quota
+      !!              relies on a full quota formalism
       !!---------------------------------------------------------------------
       !
@@ -99,5 +100,6 @@
       !
       IF( ln_timing )   CALL timing_start('p5z_prod')
-      !
+
+      ! Initialize the local arrays
       zprorcan(:,:,:) = 0._wp ; zprorcap(:,:,:) = 0._wp ; zprorcad(:,:,:) = 0._wp
       zprofed (:,:,:) = 0._wp ; zprofep (:,:,:) = 0._wp ; zprofen (:,:,:) = 0._wp
@@ -110,5 +112,6 @@
       zrespn  (:,:,:) = 0._wp ; zrespp  (:,:,:) = 0._wp ; zrespd  (:,:,:) = 0._wp 
 
-      ! Computation of the optimal production
+      ! Computation of the optimal production rates and nutrient uptake
+      ! rates. Based on a Q10 description of the thermal dependency.
       zprnut (:,:,:) = 0.6_wp * (1.0 + zpsino3 * qnnmax ) * r1_rday * tgfunc(:,:,:)
       zprnutp(:,:,:) =  0.6_wp * (1. + zpsino3 * qnpmax ) * r1_rday * tgfunc3(:,:,:)
@@ -118,4 +121,5 @@
 
       ! compute the day length depending on latitude and the day
+      ! Astronomical parameterization taken from HAMOCC3
       zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
       zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
@@ -131,5 +135,11 @@
       END DO
 
-         ! Impact of the day duration on phytoplankton growth
+      ! Impact of the day duration and light intermittency on phytoplankton growth
+      ! Intermittency is supposed to have a similar effect on production as 
+      ! day length (Shatwell et al., 2012). The correcting factor is zmxl_fac. 
+      ! zmxl_chl is the fractional day length and is used to compute the mean
+      ! PAR during daytime. The effect of mixing is computed using the 
+      ! absolute light level definition of the euphotic zone
+      ! ------------------------------------------------------------------------- 
       DO jk = 1, jpkm1
          DO jj = 1 ,jpj
@@ -154,18 +164,22 @@
       ! Maximum light intensity
       zdaylen(:,:) = MAX(1., zstrn(:,:)) / 24.
-      WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
-
+
+      ! Computation of the P-I slope for nanos, picos and diatoms
+      ! The formulation proposed by Geider et al. (1997) has been used.
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                  ! Computation of the P-I slope for nanos and diatoms
                   ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
                   zadap       = xadap * ztn / ( 2.+ ztn )
-                  !
+                  ! Nanophytoplankton
                   zpislopeadn(ji,jj,jk) = pislopen * trb(ji,jj,jk,jpnch)    &
                   &                       /( trb(ji,jj,jk,jpphy) * 12. + rtrn)
+
+                  ! Picophytoplankton
                   zpislopeadp(ji,jj,jk) = pislopep * ( 1. + zadap * EXP( -0.25 * epico(ji,jj,jk) ) )   &
                   &                       * trb(ji,jj,jk,jppch) /( trb(ji,jj,jk,jppic) * 12. + rtrn)
+
+                  ! Diatoms
                   zpislopeadd(ji,jj,jk) = pisloped * trb(ji,jj,jk,jpdch)    &
                      &                    /( trb(ji,jj,jk,jpdia) * 12. + rtrn)
@@ -176,4 +190,5 @@
 
                   ! Computation of production function for Carbon
+                  ! Actual light levels are used here 
                   !  ---------------------------------------------
                   zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) / zmxl_chl(ji,jj,jk) )  )
@@ -181,6 +196,9 @@
                   zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) / zmxl_chl(ji,jj,jk))  )
 
-                  ! Computation of production function for Chlorophyll
-                  !  -------------------------------------------------
+                  !  Computation of production function for Chlorophyll
+                  !  Mean light level in the mixed layer (when appropriate)
+                  !  is used here (acclimation is in general slower than 
+                  !  the characteristic time scales of vertical mixing)
+                  !  ------------------------------------------------------
                   zpislopen = zpislopen * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
                   zpisloped = zpisloped * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
@@ -199,9 +217,13 @@
 
                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)
+                  ! A parameterization derived from Flynn (2003) is used for the control
+                  ! when Si is not limiting which is similar to the parameterisation
+                  ! proposed by Gurney and Davidson (1999).
+                  ! -----------------------------------------------------------------------
                   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) )
@@ -225,5 +247,7 @@
       END DO
 
-      !  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
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 13232)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/trcini_pisces.F90	(revision 13233)
@@ -74,9 +74,9 @@
       USE p4zpoc          !  Remineralization of organic particles
       USE p4zligand       !  Remineralization of organic ligands
-      USE p5zlim          !  Co-limitations of differents nutrients
-      USE p5zprod         !  Growth rate of the 2 phyto groups
-      USE p5zmicro        !  Sources and sinks of microzooplankton
-      USE p5zmeso         !  Sources and sinks of mesozooplankton
-      USE p5zmort         !  Mortality terms for phytoplankton
+      USE p5zlim          !  Co-limitations of differents nutrients (QUOTA)
+      USE p5zprod         !  Growth rate of the 3 phyto groups (QUOTA)
+      USE p5zmicro        !  Sources and sinks of microzooplankton (QUOTA)
+      USE p5zmeso         !  Sources and sinks of mesozooplankton (QUOTA)
+      USE p5zmort         !  Mortality terms for phytoplankton (QUOTA)
       !
       REAL(wp), SAVE ::   sco2   =  2.312e-3_wp
@@ -253,4 +253,7 @@
 
 
+      ! Initialization of the different PISCES modules
+      ! Mainly corresponds to the namelist use
+      ! ----------------------------------------------
       CALL p4z_sink_init         !  vertical flux of particulate organic matter
       CALL p4z_opt_init          !  Optic: PAR in the water column
@@ -271,9 +274,9 @@
          & CALL p4z_ligand_init  !  remineralisation of organic ligands
 
-      IF( ln_p4z ) THEN
+      IF( ln_p4z ) THEN ! PISCES-std
          CALL p4z_mort_init      !  phytoplankton mortality 
          CALL p4z_micro_init     !  microzooplankton
          CALL p4z_meso_init      !  mesozooplankton
-      ELSE
+      ELSE ! PISCES-QUOTA
          CALL p5z_mort_init      !  phytoplankton mortality 
          CALL p5z_micro_init     !  microzooplankton
