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 12348)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 12349)
@@ -105,4 +105,10 @@
       IF( ln_ligand )  &
       & CALL p4z_ligand( kt, knt )
+      ! Update of the size of the organisms
+      IF (ln_p5z) THEN
+         sized(:,:,:) = sizeda(:,:,:)
+         sizen(:,:,:) = sizena(:,:,:)
+         sizep(:,:,:) = sizepa(:,:,:)
+      ENDIF
       !                                                             !
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
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 12348)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zint.F90	(revision 12349)
@@ -45,4 +45,5 @@
       tgfunc (:,:,:) = EXP( 0.063913 * tsn(:,:,:,jp_tem) )
       tgfunc2(:,:,:) = EXP( 0.07608  * tsn(:,:,:,jp_tem) )
+      tgfunc3(:,:,:) = EXP( 0.0825   * tsn(:,:,:,jp_tem) )
 
       ! Computation of the silicon dependant half saturation  constant for silica uptake
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 12348)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 12349)
@@ -66,5 +66,5 @@
       REAL(wp) :: zfact   , zfood, zfoodlim, zproport, zbeta
       REAL(wp) :: zmortzgoc, zfrac, zfracfe, zratio, zratio2, zfracal, zgrazcal
-      REAL(wp) :: zepsherf, zepshert, zepsherv, zgrarsig, zgraztotc, zgraztotn, zgraztotf
+      REAL(wp) :: zepsherf, zepshert, zepsherq, zepsherv, zgrarsig, zgraztotc, zgraztotn, zgraztotf
       REAL(wp) :: zgrarem2, zgrafer2, zgrapoc2, zprcaca, zmortz, zgrasrat, zgrasratn
       REAL(wp) :: zrespz, ztortz, zgrazd, zgrazz, zgrazpof
@@ -171,6 +171,5 @@
                zbeta     = MAX(0., (epsher2 - epsher2min) )
                zepsherf  = epsher2min + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta ) 
-               zepsherv  = zepsherf * zepshert 
-
+               zepsherq  = 0.5 + (1.0 - 0.5) * zepshert * ( 1.0 + 1.0 ) / ( zepshert + 1.0 )
                zgrarem2  = zgraztotc * ( 1. - zepsherv - unass2 ) &
                &         + ( 1. - epsher2 - unass2 ) / ( 1. - epsher2 ) * ztortz
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 12348)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 12349)
@@ -64,5 +64,5 @@
       REAL(wp) :: zgraze  , zdenom, zdenom2
       REAL(wp) :: zfact   , zfood, zfoodlim, zbeta
-      REAL(wp) :: zepsherf, zepshert, zepsherv, zgrarsig, zgraztotc, zgraztotn, zgraztotf
+      REAL(wp) :: zepsherf, zepshert, zepsherq, zepsherv, zgrarsig, zgraztotc, zgraztotn, zgraztotf
       REAL(wp) :: zgrarem, zgrafer, zgrapoc, zprcaca, zmortz
       REAL(wp) :: zrespz, ztortz, zgrasrat, zgrasratn
@@ -131,6 +131,6 @@
                zbeta     = MAX(0., (epsher - epshermin) )
                zepsherf  = epshermin + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
-               zepsherv  = zepsherf * zepshert 
-
+               zepsherq  = 0.5 + (1.0 - 0.5) * zepshert * ( 1.0 + 1.0 ) / ( zepshert + 1.0 )
+               zepsherv  = zepsherf * zepshert * zepsherq
                zgrafer   = zgraztotc * MAX( 0. , ( 1. - unass ) * zgrasrat - ferat3 * zepsherv ) 
                zgrarem   = zgraztotc * ( 1. - zepsherv - unass )
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 12348)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 12349)
@@ -26,34 +26,34 @@
 
    !! * Shared module variables
-   REAL(wp), PUBLIC ::  concpno3    !:  NO3, PO4 half saturation   
-   REAL(wp), PUBLIC ::  concpnh4    !:  NH4 half saturation for phyto  
-   REAL(wp), PUBLIC ::  concnpo4    !:  NH4 half saturation for diatoms
-   REAL(wp), PUBLIC ::  concppo4    !:  NH4 half saturation for diatoms
-   REAL(wp), PUBLIC ::  concdpo4    !:  NH4 half saturation for diatoms
-   REAL(wp), PUBLIC ::  concpfer    !:  Iron half saturation for nanophyto 
+   REAL(wp), PUBLIC ::  concpno3    !:  NO3 half saturation for picophyto  
+   REAL(wp), PUBLIC ::  concpnh4    !:  NH4 half saturation for picophyto
+   REAL(wp), PUBLIC ::  concnpo4    !:  PO4 half saturation for nanophyto
+   REAL(wp), PUBLIC ::  concppo4    !:  PO4 half saturation for picophyto
+   REAL(wp), PUBLIC ::  concdpo4    !:  PO4 half saturation for diatoms
+   REAL(wp), PUBLIC ::  concpfer    !:  Iron half saturation for picophyto
    REAL(wp), PUBLIC ::  concbpo4    !:  PO4 half saturation for bacteria
-   REAL(wp), PUBLIC ::  xsizepic    !:  Minimum size criteria for diatoms
-   REAL(wp), PUBLIC ::  xsizerp     !:  Size ratio for nanophytoplankton
+   REAL(wp), PUBLIC ::  xsizepic    !:  Minimum size criteria for picophyto
+   REAL(wp), PUBLIC ::  xsizerp     !:  Size ratio for picophytoplankton
    REAL(wp), PUBLIC ::  qfnopt      !:  optimal Fe quota for nanophyto
-   REAL(wp), PUBLIC ::  qfpopt      !:  optimal Fe quota for nanophyto
+   REAL(wp), PUBLIC ::  qfpopt      !:  optimal Fe quota for picophyto
    REAL(wp), PUBLIC ::  qfdopt      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qnnmin      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qnnmax      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qpnmin      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qpnmax      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qnpmin      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qnpmax      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qppmin      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qppmax      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qndmin      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qndmax      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qpdmin      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qpdmax      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qfnmax      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qfpmax      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  qfdmax      !:  optimal Fe quota for diatoms
-   REAL(wp), PUBLIC ::  zpsinh4
-   REAL(wp), PUBLIC ::  zpsino3
-   REAL(wp), PUBLIC ::  zpsiuptk
+   REAL(wp), PUBLIC ::  qnnmin      !:  minimum N  quota for nanophyto
+   REAL(wp), PUBLIC ::  qnnmax      !:  maximum N quota for nanophyto
+   REAL(wp), PUBLIC ::  qpnmin      !:  minimum P quota for nanophyto
+   REAL(wp), PUBLIC ::  qpnmax      !:  maximum P quota for nanophyto
+   REAL(wp), PUBLIC ::  qnpmin      !:  minimum N quota for nanophyto
+   REAL(wp), PUBLIC ::  qnpmax      !:  maximum N quota for nanophyto
+   REAL(wp), PUBLIC ::  qppmin      !:  minimum P quota for nanophyto
+   REAL(wp), PUBLIC ::  qppmax      !:  maximum P quota for nanophyto
+   REAL(wp), PUBLIC ::  qndmin      !:  minimum N quota for diatoms
+   REAL(wp), PUBLIC ::  qndmax      !:  maximum N quota for diatoms
+   REAL(wp), PUBLIC ::  qpdmin      !:  minimum P quota for diatoms
+   REAL(wp), PUBLIC ::  qpdmax      !:  maximum P quota for diatoms
+   REAL(wp), PUBLIC ::  qfnmax      !:  maximum Fe quota for nanophyto
+   REAL(wp), PUBLIC ::  qfpmax      !:  maximum Fe quota for picophyto
+   REAL(wp), PUBLIC ::  qfdmax      !:  maximum Fe quota for diatoms
+   REAL(wp), PUBLIC ::  zpsinh4     !:  respiration cost of NH4 assimilation
+   REAL(wp), PUBLIC ::  zpsino3     !:  respiration cost of NO3 assimilation
+   REAL(wp), PUBLIC ::  zpsiuptk    !:  Mean respiration cost
 
    !!*  Allometric variations of the quotas
@@ -71,5 +71,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpdmax    !: ???
 
-   !!* Phytoplankton limitation terms
+   !!* Phytoplankton nutrient limitation terms
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicono3   !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpiconh4   !: ???
@@ -82,4 +82,7 @@
    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
@@ -104,7 +107,14 @@
       !!
       !! ** Purpose :   Compute the co-limitations by the various nutrients
-      !!              for the various phytoplankton species
+      !!                for the various phytoplankton species. Quota based
+      !!                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. 
       !!
-      !! ** Method  : - ???
+      !! ** Method  : Quota based approach. The quota model is derived from 
+      !!              theoretical models by Pahlow and Oschlies (2009) and 
+      !!              Flynn (2001). Various adaptations from several publications
+      !!              by these authors have been also adopted.
       !!---------------------------------------------------------------------
       !
@@ -114,5 +124,5 @@
       REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zno3, zferlim
       REAL(wp) ::   z1_trndia, z1_trnpic, z1_trnphy, ztem1, ztem2, zetot1
-      REAL(wp) ::   zratio, zration, zratiof, znutlim, zfalim
+      REAL(wp) ::   zratio, zration, zratiof, znutlim, zfalim, zzpsiuptk
       REAL(wp) ::   zconc1d, zconc1dnh4, zconc0n, zconc0nnh4, zconc0npo4, zconc0dpo4
       REAL(wp) ::   zconc0p, zconc0pnh4, zconc0ppo4, zconcpfe, zconcnfe, zconcdfe
@@ -122,4 +132,6 @@
       REAL(wp) ::   zfvn, zfvp, zfvf, zsizen, zsizep, zsized, znanochl, zpicochl, zdiatchl
       REAL(wp) ::   zqfemn, zqfemp, zqfemd, zbactno3, zbactnh4
+      REAL(wp) ::   zlim1f, zsizetmp
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: xlimnpn, xlimnpp, xlimnpd
       !!---------------------------------------------------------------------
       !
@@ -127,4 +139,5 @@
       !
       zratchl = 6.0
+      sizena(:,:,:) = 0.0  ;  sizepa(:,:,:) = 0.0  ;  sizeda(:,:,:) = 0.0
       !
       DO jk = 1, jpkm1
@@ -157,5 +170,5 @@
                zconc0dpo4        = concdpo4 * zsized
 
-               zsizep            = 1.
+               zsizep            = sizep(ji,jj,jk)**0.81
                zconcpfe          = concpfer * zsizep
                zconc0p           = concpno3 * zsizep
@@ -163,5 +176,5 @@
                zconc0ppo4        = concppo4 * zsizep
 
-               zsizen            = 1.
+               zsizen            = sizen(ji,jj,jk)**0.81
                zconcnfe          = concnfer * zsizen
                zconc0n           = concnno3 * zsizen
@@ -172,10 +185,10 @@
                ! From Talmy et al. (2014) and Maranon et al. (2013)
                ! -------------------------------------------------------
-               xqnnmin(ji,jj,jk) = qnnmin
+               xqnnmin(ji,jj,jk) = qnnmin * sizen(ji,jj,jk)**(-0.3)
                xqnnmax(ji,jj,jk) = qnnmax
-               xqndmin(ji,jj,jk) = qndmin * sized(ji,jj,jk)**(-0.27) 
+               xqndmin(ji,jj,jk) = qndmin * sized(ji,jj,jk)**(-0.3)
                xqndmax(ji,jj,jk) = qndmax
-               xqnpmin(ji,jj,jk) = qnpmin
-               xqnpmax(ji,jj,jk) = qnpmax
+               xqnpmin(ji,jj,jk) = qnpmin * sizep(ji,jj,jk)**(-0.48)
+               xqnpmax(ji,jj,jk) = qnpmax * sizep(ji,jj,jk)**(-0.21)
 
                ! Computation of the optimal allocation parameters
@@ -236,14 +249,21 @@
                !
                zration = trb(ji,jj,jk,jpnph) * z1_trnphy
-               zration = MIN(xqnnmax(ji,jj,jk), MAX( 2. * xqnnmin(ji,jj,jk), zration ))
-               fvnuptk(ji,jj,jk) = 1. / zpsiuptk * rno3 * 2. * xqnnmin(ji,jj,jk) / (zration + rtrn)  &
+               zration = MIN(xqnnmax(ji,jj,jk), MAX( xqnnmin(ji,jj,jk), zration ))
+               zzpsiuptk = xqnnmin(ji,jj,jk) * rno3 / zpsiuptk**2
+               fvnuptk(ji,jj,jk) = 1. / zzpsiuptk * xqnnmin(ji,jj,jk) / (zration + rtrn)  &
                &                   * MAX(0., (1. - zratchl * znanochl / 12. ) )
                !
-               zlim1    = max(0., (zration - 2. * xqnnmin(ji,jj,jk) )  &
-               &          / (xqnnmax(ji,jj,jk) - 2. * xqnnmin(ji,jj,jk) ) ) * xqnnmax(ji,jj,jk)  &
+               zlim1  = max(0., (zration - xqnnmin(ji,jj,jk) )  &
+               &          / (xqnnmax(ji,jj,jk) - xqnnmin(ji,jj,jk) ) ) * xqnnmax(ji,jj,jk)  &
                &          / (zration + rtrn)
-               zlim3    = MAX( 0.,( zratiof - zqfemn ) / qfnopt ) 
-               xlimnfe(ji,jj,jk) = MIN( 1., zlim3 )
-               xlimphy(ji,jj,jk) = MIN( 1., zlim1, zlim3 )
+               !  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 )
+               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
@@ -267,14 +287,21 @@
                !
                zration   = trb(ji,jj,jk,jpnpi) * z1_trnpic
-               zration = MIN(xqnpmax(ji,jj,jk), MAX( 2. * xqnpmin(ji,jj,jk), zration ))
-               fvpuptk(ji,jj,jk) = 1. / zpsiuptk * rno3 * 2. * xqnpmin(ji,jj,jk) / (zration + rtrn)  &
+               zration = MIN(xqnpmax(ji,jj,jk), MAX( xqnpmin(ji,jj,jk), zration ))
+               zzpsiuptk = xqnpmin(ji,jj,jk) * rno3 / zpsiuptk**2
+               fvpuptk(ji,jj,jk) = 1. / zzpsiuptk * xqnpmin(ji,jj,jk) / (zration + rtrn)  &
                &                   * MAX(0., (1. - zratchl * zpicochl / 12. ) ) 
                !
-               zlim1    = max(0., (zration - 2. * xqnpmin(ji,jj,jk) )  &
-               &          / (xqnpmax(ji,jj,jk) - 2. * xqnpmin(ji,jj,jk) ) ) * xqnpmax(ji,jj,jk)  &
+               zlim1    = max(0., (zration - xqnpmin(ji,jj,jk) )  &
+               &          / (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
+               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 )
-               xlimpfe(ji,jj,jk) = MIN( 1., zlim3 )
-               xlimpic(ji,jj,jk) = MIN( 1., zlim1, zlim3 )
+               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
@@ -298,16 +325,24 @@
                !
                zration   = trb(ji,jj,jk,jpndi) * z1_trndia
-               zration = MIN(xqndmax(ji,jj,jk), MAX( 2. * xqndmin(ji,jj,jk), zration ))
-               fvduptk(ji,jj,jk) = 1. / zpsiuptk * rno3 * 2. * xqndmin(ji,jj,jk) / (zration + rtrn)   &
+               zration   = MIN(xqndmax(ji,jj,jk), MAX( xqndmin(ji,jj,jk), zration ))
+               zzpsiuptk = xqndmin(ji,jj,jk) * rno3 / zpsiuptk**2
+               fvduptk(ji,jj,jk) = 1. / zzpsiuptk * xqndmin(ji,jj,jk) / (zration + rtrn)   &
                &                   * MAX(0., (1. - zratchl * zdiatchl / 12. ) ) 
                !
-               zlim1    = max(0., (zration - 2. * xqndmin(ji,jj,jk) )    &
-               &          / (xqndmax(ji,jj,jk) - 2. * xqndmin(ji,jj,jk) ) )   &
+               zlim1    = max(0., (zration - xqndmin(ji,jj,jk) )    &
+               &          / (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
+               zlim1f   = max(0., (1.077 - xqndmin(ji,jj,jk) )    &
+               &          / (xqndmax(ji,jj,jk) - xqndmin(ji,jj,jk) ) )   &
+               &          * xqndmax(ji,jj,jk)
                zlim3    = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi(ji,jj) )
                zlim4    = MAX( 0., ( zratiof - zqfemd ) / qfdopt )
                xlimdfe(ji,jj,jk) = MIN( 1., zlim4 )
                xlimdia(ji,jj,jk) = MIN( 1., zlim1, zlim3, zlim4 )
+               xlimdias(ji,jj,jk) = MIN (1.0, zlim1 / (zlim1f + rtrn ), zlim3, zlim4 )
                xlimsi(ji,jj,jk)  = MIN( zlim1, zlim4 )
+               xlimnpd(ji,jj,jk) = MIN( 1., zlim1 )
             END DO
          END DO
@@ -324,45 +359,51 @@
                ! Size estimation of nanophytoplankton
                ! ------------------------------------
-               zfvn = 2. * fvnuptk(ji,jj,jk)
-               sizen(ji,jj,jk) = MAX(1., MIN(xsizern, 1.0 / ( MAX(rtrn, zfvn) ) ) )
-
+               zcoef = trb(ji,jj,jk,jpphy) - MIN(xsizephy, trb(ji,jj,jk,jpphy) )
+               sizena(ji,jj,jk) = 1. + ( xsizern -1.0 ) * zcoef / ( xsizephy + zcoef )
                ! N/P ratio of nanophytoplankton
                ! ------------------------------
-               zfuptk = 0.23 * zfvn
-               zrpho = 2.24 * trb(ji,jj,jk,jpnch) / ( trb(ji,jj,jk,jpnph) * rno3 * 15. + rtrn )
-               zrass = 1. - 0.2 - zrpho - zfuptk
-               xqpnmax(ji,jj,jk) = ( zfuptk + zrpho ) * 0.0128 * 16. + zrass * 1./ 7.2 * 16.
-               xqpnmax(ji,jj,jk) = xqpnmax(ji,jj,jk) * trb(ji,jj,jk,jpnph) / ( trb(ji,jj,jk,jpphy) + rtrn ) + 0.13
-               xqpnmin(ji,jj,jk) = 0.13 + 0.23 * 0.0128 * 16.
+               zfuptk = 0.2 + 0.12 / ( 3.0 * sizen(ji,jj,jk) + rtrn )
+               zrpho  = 1.54 * trb(ji,jj,jk,jpnch) / ( trb(ji,jj,jk,jpnph) * rno3 * 14. + rtrn )
+               zrass = MAX(0.62/4., ( 1. - zrpho - zfuptk ) * xlimnpn(ji,jj,jk) )
+               xqpnmin(ji,jj,jk) = ( 0.0 + 0.0078 + 0.62/4. * 0.0783 * xqnnmin(ji,jj,jk) ) * 16.
+               xqpnmax(ji,jj,jk) = ( zrpho * 0.0128 + zrass * 0.0783 ) * 16.
+               xqpnmax(ji,jj,jk) = xqpnmax(ji,jj,jk) * trb(ji,jj,jk,jpnph) / ( trb(ji,jj,jk,jpphy) + rtrn )  &
+               &      + (0.033 + 0.0078 ) * 16.
+               xqpnmax(ji,jj,jk) = MIN( qpnmax, xqpnmax(ji,jj,jk) )
+
 
                ! Size estimation of picophytoplankton
                ! ------------------------------------
-               zfvn = 2. * fvpuptk(ji,jj,jk)
-               sizep(ji,jj,jk) = MAX(1., MIN(xsizerp, 1.0 / ( MAX(rtrn, zfvn) ) ) )
+               zcoef = trb(ji,jj,jk,jppic) - MIN(xsizepic, trb(ji,jj,jk,jppic) )
+               sizepa(ji,jj,jk) = 1. + ( xsizerp -1.0 ) * zcoef / ( xsizepic + zcoef )
 
                ! N/P ratio of picophytoplankton
                ! ------------------------------
-               zfuptk = 0.35 * zfvn
-               zrpho = 2.24 * trb(ji,jj,jk,jppch) / ( trb(ji,jj,jk,jpnpi) * rno3 * 15. + rtrn )
-               zrass = 1. - 0.4 - zrpho - zfuptk
-               xqppmax(ji,jj,jk) =  (zrpho + zfuptk) * 0.0128 * 16. + zrass * 1./ 9. * 16.
-               xqppmax(ji,jj,jk) = xqppmax(ji,jj,jk) * trb(ji,jj,jk,jpnpi) / ( trb(ji,jj,jk,jppic) + rtrn ) + 0.13
-               xqppmin(ji,jj,jk) = 0.13
+               zfuptk = 0.2 + 0.12 / ( 0.5 * sizep(ji,jj,jk) + rtrn )
+               zrpho = 1.54 * trb(ji,jj,jk,jppch) / ( trb(ji,jj,jk,jpnpi) * rno3 * 14. + rtrn )
+               zrass = MAX(0.4/4., ( 1. - zrpho - zfuptk ) * xlimnpp(ji,jj,jk) )
+               xqppmin(ji,jj,jk) = ( (0.0 + 0.0078 ) + 0.4/4. * 0.0517 * xqnpmin(ji,jj,jk) ) * 16.
+               xqppmax(ji,jj,jk) = ( zrpho * 0.0128 + zrass * 0.0517 ) * 16.
+               xqppmax(ji,jj,jk) = xqppmax(ji,jj,jk) * trb(ji,jj,jk,jpnpi) / ( trb(ji,jj,jk,jppic) + rtrn ) &
+               &      +  (0.033 + 0.0078 ) * 16
+               xqppmax(ji,jj,jk) = MIN( qppmax, xqppmax(ji,jj,jk) )
 
                ! Size estimation of diatoms
                ! --------------------------
-               zfvn = 2. * fvduptk(ji,jj,jk)
-               sized(ji,jj,jk) = MAX(1., MIN(xsizerd, 1.0 / ( MAX(rtrn, zfvn) ) ) )
                zcoef = trb(ji,jj,jk,jpdia) - MIN(xsizedia, trb(ji,jj,jk,jpdia) )
-               sized(ji,jj,jk) = 1. + xsizerd * zcoef *1E6 / ( 1. + zcoef * 1E6 )
+               sized(ji,jj,jk) = 1. + ( xsizerd - 1.0 ) * zcoef / ( xsizedia + zcoef )
+               sizeda(ji,jj,jk) = 1. + ( xsizerd - 1.0 ) * zcoef / ( xsizedia + zcoef )
 
                ! N/P ratio of diatoms
                ! --------------------
-               zfuptk = 0.2 * zfvn
-               zrpho = 2.24 * trb(ji,jj,jk,jpdch) / ( trb(ji,jj,jk,jpndi) * rno3 * 15. + rtrn )
-               zrass = 1. - 0.2 - zrpho - zfuptk
-               xqpdmax(ji,jj,jk) = ( zfuptk + zrpho ) * 0.0128 * 16. + zrass * 1./ 7.2 * 16.
-               xqpdmax(ji,jj,jk) = xqpdmax(ji,jj,jk) * trb(ji,jj,jk,jpndi) / ( trb(ji,jj,jk,jpdia) + rtrn ) + 0.13
-               xqpdmin(ji,jj,jk) = 0.13 + 0.2 * 0.0128 * 16.
+               zfuptk = 0.2 + 0.12 / ( 5.0 * sized(ji,jj,jk) + rtrn )
+               zrpho = 1.54 * trb(ji,jj,jk,jpdch) / ( trb(ji,jj,jk,jpndi) * rno3 * 14. + rtrn )
+               zrass = MAX(0.66/4., ( 1. - zrpho - zfuptk ) * xlimnpd(ji,jj,jk) )
+
+               xqpdmin(ji,jj,jk) = ( ( 0.0 + 0.0078 ) + 0.66/4. * 0.0783 *  xqndmin(ji,jj,jk) ) * 16.
+               xqpdmax(ji,jj,jk) = ( zrpho * 0.0128 + zrass * 0.0783 ) * 16.
+               xqpdmax(ji,jj,jk) = xqpdmax(ji,jj,jk) * trb(ji,jj,jk,jpndi) / ( trb(ji,jj,jk,jpdia) + rtrn ) &
+               &      + ( 0.0078 + 0.033 ) * 16.
+               xqpdmax(ji,jj,jk) = MIN(qpdmax, xqpdmax(ji,jj,jk) )
 
             END DO
@@ -384,6 +425,5 @@
                zetot1 = MAX( 0., etot(ji,jj,jk) - 1.) / ( 4. + etot(ji,jj,jk) ) * 20. / ( 20. + etot(ji,jj,jk) ) 
 
-!               xfracal(ji,jj,jk) = caco3r * MIN( zlim1, zlim2, zlim3 )                  &
-               xfracal(ji,jj,jk) = caco3r                 &
+               xfracal(ji,jj,jk) = caco3r * MIN( zlim1, zlim2, zlim3 )    &
                &                   * ztem1 / ( 1. + ztem1 ) * MAX( 1., trb(ji,jj,jk,jpphy)*1E6 )   &
                   &                * ( 1. + EXP(-ztem2 * ztem2 / 25. ) )         &
@@ -522,7 +562,7 @@
       ENDIF
       !
-      zpsino3 = 2.3 * rno3
-      zpsinh4 = 1.8 * rno3
-      zpsiuptk = 2.3 * rno3
+      zpsino3  = 2.3 * rno3
+      zpsinh4  = 1.8 * rno3
+      zpsiuptk = 1.0 / 6.625
       !
       nitrfac (:,:,:) = 0._wp
@@ -547,4 +587,6 @@
          &      xpicofer(jpi,jpj,jpk), xlimpfe (jpi,jpj,jpk),       &
          &      fvnuptk (jpi,jpj,jpk), fvduptk (jpi,jpj,jpk),       &
+         &      xlimphys(jpi,jpj,jpk), xlimdias(jpi,jpj,jpk),       &
+         &      xlimpics(jpi,jpj,jpk),                              &
          &      fvpuptk (jpi,jpj,jpk), xlimpic (jpi,jpj,jpk),    STAT=ierr(1) )
          !
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 12348)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 12349)
@@ -72,5 +72,5 @@
       REAL(wp) :: zgraze2, zdenom, zfact, zfood, zfoodlim, zproport
       REAL(wp) :: zmortzgoc, zfracc, zfracn, zfracp, zfracfe, zratio, zratio2
-      REAL(wp) :: zepsherf, zepshert, zepsherv, zrespirc, zrespirn, zrespirp, zbasresb, zbasresi
+      REAL(wp) :: zepsherf, zepshert, zepsherq, zepsherv, zrespirc, zrespirn, zrespirp, zbasresb, zbasresi
       REAL(wp) :: zgraztotc, zgraztotn, zgraztotp, zgraztotf, zbasresn, zbasresp, zbasresf
       REAL(wp) :: zgradoc, zgradon, zgradop, zgratmp, zgradoct, zgradont, zgrareft, zgradopt
@@ -85,5 +85,5 @@
       REAL(wp) :: zgrazffnp, zgrazffng, zgrazffpp, zgrazffpg
       CHARACTER (len=25) :: charout
-      REAL(wp) :: zrfact2, zmetexcess
+      REAL(wp) :: zrfact2, zmetexcess, zsigma, zdiffdn
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing, zfezoo2
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d, zz2ligprod
@@ -140,15 +140,18 @@
                !   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
+               !   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. Thus, instead of a square
-               !   a 1.5 power value is used which decreases the pressure on the
-               !   most abundant species
+               !   too strongly on winning organisms. We use a generalized
+               !   switching parameterization proposed by Morozov and 
+               !   Petrovskii (2013)
                !   ------------------------------------------------------------  
-               ztmp1 = xpref2n * zcompaph**1.5
-               ztmp2 = xpref2m * zcompames**1.5
-               ztmp3 = xpref2c * zcompapoc**1.5
-               ztmp4 = xpref2d * zcompadi**1.5
-               ztmp5 = xpref2z * zcompaz**1.5
+               zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
+               zsigma = 0.5 + 1.0 * zsigma
+               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 )
+               ztmp2 = xpref2m * zcompames**2
+               ztmp3 = xpref2c * zcompapoc**2
+               ztmp4 = xpref2d * zcompadi * ( zdiffdn * zcompadi + zcompaph )
+               ztmp5 = xpref2z * zcompaz**2
                ztmptot = ztmp1 + ztmp2 + ztmp3 + ztmp4 + ztmp5 + rtrn
                ztmp1 = ztmp1 / ztmptot
@@ -238,6 +241,6 @@
                zbeta     = MAX(0., (epsher2 - epsher2min) )
                zepsherf  = epsher2min + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
-               zepsherv  = zepsherf * zepshert
-
+               zepsherq  = 0.5 + (1.0 - 0.5) * zepshert * ( 1.0 + 1.0 ) / ( zepshert + 1.0 )
+               zepsherv  = zepsherf * zepshert * zepsherq
                !   Respiration of mesozooplankton
                !   Excess carbon in the food is used preferentially
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 12348)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmicro.F90	(revision 12349)
@@ -75,5 +75,5 @@
       REAL(wp) :: zcompapi, zgraze  , zdenom, zfact, zfood, zfoodlim
       REAL(wp) :: ztmp1, ztmp2, ztmp3, ztmp4, ztmp5, ztmptot
-      REAL(wp) :: zepsherf, zepshert, zepsherv, zrespirc, zrespirn, zrespirp, zbasresb, zbasresi
+      REAL(wp) :: zepsherf, zepshert, zepsherq, zepsherv, zrespirc, zrespirn, zrespirp, zbasresb, zbasresi
       REAL(wp) :: zgraztotc, zgraztotn, zgraztotp, zgraztotf, zbasresn, zbasresp, zbasresf
       REAL(wp) :: zgradoc, zgradon, zgradop, zgraref, zgradoct, zgradont, zgradopt, zgrareft
@@ -84,4 +84,5 @@
       REAL(wp) :: zgrazdc, zgrazdn, zgrazdp, zgrazdf, zgraznf, zgrazz
       REAL(wp) :: zgrazpc, zgrazpn, zgrazpp, zgrazpf, zbeta, zrfact2, zmetexcess
+      REAL(wp) :: zsigma, zdiffdn, zdiffpn, zdiffdp, zproport, zproport2
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing, zfezoo
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d, zzligprod
@@ -104,5 +105,8 @@
                zcompaz = MAX( ( trb(ji,jj,jk,jpzoo) - 1.e-9 ), 0.e0 )
                zfact   = xstep * tgfunc2(ji,jj,jk) * zcompaz
-
+               ! Proportion of nano and diatoms that are within the size range
+               ! accessible to microzooplankton. 
+               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
                !   -----------------------------------------------------
@@ -118,6 +122,6 @@
                !   A threshold can be specified in the namelist
                !   --------------------------------------------
-               zcompadi  = MIN( MAX( ( trb(ji,jj,jk,jpdia) - xthreshdia ), 0.e0 ), xsizedia )
-               zcompaph  = MAX( ( trb(ji,jj,jk,jpphy) - xthreshphy ), 0.e0 )
+               zcompadi  = zproport * MAX( ( trb(ji,jj,jk,jpdia) - xthreshdia ), 0.e0 )
+               zcompaph  = zproport2 * MAX( ( trb(ji,jj,jk,jpphy) - xthreshphy ), 0.e0 )
                zcompaz   = MAX( ( trb(ji,jj,jk,jpzoo) - xthreshzoo ), 0.e0 )
                zcompapi  = MAX( ( trb(ji,jj,jk,jppic) - xthreshpic ), 0.e0 )
@@ -140,9 +144,14 @@
                !   most abundant species
                !   ------------------------------------------------------------  
-               ztmp1 = xprefn * zcompaph**1.5
-               ztmp2 = xprefp * zcompapi**1.5
-               ztmp3 = xprefc * zcompapoc**1.5
-               ztmp4 = xprefd * zcompadi**1.5
-               ztmp5 = xprefz * zcompaz**1.5
+               zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
+               zsigma = 0.5 + 1.0*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)
+               zdiffdp = exp( -ABS(log(0.5 * sizep(ji,jj,jk) / (5.0 * sized(ji,jj,jk) + rtrn )) )**2 / zsigma**2)
+               ztmp1 = xprefn * zcompaph * ( zcompaph + zdiffdn * zcompadi + zdiffpn * zcompapi ) / ( 1.0 + zdiffdn + zdiffpn )
+               ztmp2 = xprefp * zcompapi * ( zcompapi + zdiffpn * zcompaph + zdiffdp * zcompadi ) / ( 1.0 + zdiffpn + zdiffdp )
+               ztmp3 = xprefc * zcompapoc**2
+               ztmp4 = xprefd * zcompadi * ( zdiffdp * zcompapi + zdiffdn * zcompaph + zcompadi ) / ( 1.0 + zdiffdn + zdiffdp )
+               ztmp5 = xprefz * zcompaz**2
                ztmptot = ztmp1 + ztmp2 + ztmp3 + ztmp4 + ztmp5 + rtrn
                ztmp1 = ztmp1 / ztmptot
@@ -192,5 +201,6 @@
                zbeta     = MAX( 0., (epsher - epshermin) )
                zepsherf  = epshermin + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
-               zepsherv  = zepsherf * zepshert
+               zepsherq  = 0.5 + (1.0 - 0.5) * zepshert * ( 1.0 + 1.0 ) / ( zepshert + 1.0 )
+               zepsherv  = zepsherf * zepshert * zepsherq
 
                !   Respiration of microzooplankton
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 12348)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zprod.F90	(revision 12349)
@@ -71,15 +71,15 @@
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zsilfac, znanotot, zpicotot, zdiattot, zconctemp, zconctemp2
-      REAL(wp) ::   zration, zratiop, zratiof, zmax, zmax2, zsilim, ztn, zadap
+      REAL(wp) ::   zration, zratiop, zratiof, zmax, zsilim, ztn, zadap
       REAL(wp) ::   zpronmax, zpropmax, zprofmax, zrat
       REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zprontot, zproptot, zprodtot
       REAL(wp) ::   zprnutmax, zdocprod, zprochln, zprochld, zprochlp
-      REAL(wp) ::   zpislopen, zpislopep, zpisloped, thetannm_n, thetandm_n, thetanpm_n
+      REAL(wp) ::   zpislopen, zpislopep, zpisloped
       REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
-      REAL(wp) ::   zfact, zrfact2
+      REAL(wp) ::   zfact, zrfact2, zmaxsi, zratiosi, zsizetmp, zlimfac
       CHARACTER (len=25) :: charout
       REAL(wp), DIMENSION(jpi,jpj    ) :: zmixnano, zmixpico, zmixdiat, zstrn
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadp, zpislopeadd
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprnut, zprmaxp, zprmaxn, zprmaxd
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprnut, zprnutp, zprmaxp, zprmaxn, zprmaxd
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprbio, zprpic, zprdia, zysopt
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprchln, zprchlp, zprchld
@@ -111,8 +111,9 @@
 
       ! Computation of the optimal production
-      zprnut (:,:,:) = 0.65_wp * r1_rday * tgfunc(:,:,:)
-      zprmaxn(:,:,:) = ( 0.65_wp * (1. + zpsino3 * qnpmax ) ) * r1_rday * tgfunc(:,:,:)
-      zprmaxp(:,:,:) = 0.5 / 0.65 * zprmaxn(:,:,:) 
-      zprmaxd(:,:,:) = zprmaxn(:,:,:) 
+      zprnut (:,:,:) = 0.6_wp * (1.0 + zpsino3 * qnnmax ) * r1_rday * tgfunc(:,:,:)
+      zprnutp(:,:,:) =  0.6_wp * (1. + zpsino3 * qnpmax ) * r1_rday * tgfunc3(:,:,:)
+      zprmaxn(:,:,:) = ( 0.6_wp * (1. + zpsino3 * qnnmax ) ) * r1_rday * tgfunc(:,:,:)
+      zprmaxd(:,:,:) = ( 0.6_wp * (1. + zpsino3 * qndmax ) ) * r1_rday * tgfunc(:,:,:)
+      zprmaxp(:,:,:) = ( 0.4_wp * (1. + zpsino3 * qnpmax ) ) * r1_rday * tgfunc3(:,:,:)
 
       ! compute the day length depending on latitude and the day
@@ -140,5 +141,5 @@
                   ENDIF
                   zmxl_chl(ji,jj,jk) = zval / 24.
-                  zmxl_fac(ji,jj,jk) = 1.5 * zval / ( 12. + zval )
+                  zmxl_fac(ji,jj,jk) = 1.0 - exp( -0.26 * zval )
                ENDIF
             END DO
@@ -176,7 +177,7 @@
                   ! Computation of production function for Carbon
                   !  ---------------------------------------------
-                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
-                  zprpic(ji,jj,jk) = zprpic(ji,jj,jk) * ( 1.- EXP( -zpislopep * epico(ji,jj,jk) )  )
-                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
+                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) / zmxl_chl(ji,jj,jk) )  )
+                  zprpic(ji,jj,jk) = zprpic(ji,jj,jk) * ( 1.- EXP( -zpislopep * epico(ji,jj,jk) / zmxl_chl(ji,jj,jk))  )
+                  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
@@ -197,5 +198,5 @@
             DO ji = 1, jpi
 
-                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                   !    Si/C of diatoms
                   !    ------------------------
@@ -205,13 +206,19 @@
                   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) )
-                  zsilfac = 3.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim - 0.5 ) )  ) + 1.e0
                   zsiborn = trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil)
                   IF (gphit(ji,jj) < -30 ) THEN
-                    zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
+                    zsilfac2 = 1. + 1. * zsiborn / ( zsiborn + xksi2**3 )
                   ELSE
-                    zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 )
+                    zsilfac2 = 1.
                   ENDIF
-                  zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2
-              ENDIF
+                  zratiosi = 1.0 - trb(ji,jj,jk,jpdsi) / ( trb(ji,jj,jk,jpdia) + rtrn ) / ( zsilfac2 * grosip * 3.0 + rtrn )
+                  zratiosi = MAX(0., MIN(1.0, zratiosi) )
+                  zmaxsi  = (1.0 + 0.1**4) * zratiosi**4 / ( zratiosi**4 + 0.1**4 )
+                  IF ( xlimsi(ji,jj,jk) /= xlimdia(ji,jj,jk) ) THEN
+                     zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zmaxsi
+                  ELSE
+                     zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zsilim**0.75 * zmaxsi
+                  ENDIF
+               ENDIF
             END DO
          END DO
@@ -226,4 +233,5 @@
                zprdia(ji,jj,jk)  = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
                zprnut(ji,jj,jk)  = zprnut(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
+               zprnutp(ji,jj,jk) = zprnutp(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
             END DO
          END DO
@@ -237,5 +245,9 @@
                   !  production terms for nanophyto.
                   zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trb(ji,jj,jk,jpphy) * rfact2
-                  !
+                  ! size computation
+                  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
                   zration = trb(ji,jj,jk,jpnph) / ( trb(ji,jj,jk,jpphy) + rtrn )
                   zratiop = trb(ji,jj,jk,jppph) / ( trb(ji,jj,jk,jpphy) + rtrn )
@@ -243,22 +255,22 @@
                   zprnutmax = zprnut(ji,jj,jk) * fvnuptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jpphy) * rfact2
                   ! Uptake of nitrogen
-                  zrat = MIN( 1., zration / (xqnnmax(ji,jj,jk) + rtrn) ) 
-                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
+                  zrat = 1.0 - MIN( 1., zration / (xqnnmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
                   zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqpnmin(ji,jj,jk) )   &
                   &          / ( xqpnmax(ji,jj,jk) - xqpnmin(ji,jj,jk) + rtrn ), xlimnfe(ji,jj,jk) ) )
-                  zpronewn(ji,jj,jk) = zpronmax * zdaylen(ji,jj) * xnanono3(ji,jj,jk)
+                  zpronewn(ji,jj,jk) = zpronmax * xnanono3(ji,jj,jk)
                   zproregn(ji,jj,jk) = zpronmax * xnanonh4(ji,jj,jk)
                   ! Uptake of phosphorus
-                  zrat = MIN( 1., zratiop / (xqpnmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
-                  zpropmax = zprnutmax * zmax * xlimnfe(ji,jj,jk)
+                  zrat = 1.0 - MIN( 1., zratiop / (xqpnmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
+                  zpropmax = zprnutmax * zmax * xlimnfe(ji,jj,jk) * 16. / 10.
                   zpropo4n(ji,jj,jk) = zpropmax * xnanopo4(ji,jj,jk)
                   zprodopn(ji,jj,jk) = zpropmax * xnanodop(ji,jj,jk)
                   ! Uptake of iron
-                  zrat = MIN( 1., zratiof / qfnmax )
-                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
-                  zprofmax = zprnutmax * qfnmax * zmax
-                  zprofen(ji,jj,jk) = zprofmax * xnanofer(ji,jj,jk) * ( 3. - 2.4 * xlimnfe(ji,jj,jk)    &
-                  &          / ( xlimnfe(ji,jj,jk) + 0.2 ) ) * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn  &
+                  zrat = 1.0 - MIN( 1., zratiof / qfnmax )
+                  zmax = MAX(0., MIN(1., zrat**2/ (0.05**2 + zrat**2) ) )
+                  zprofmax = zprnutmax * qfnmax * zmax 
+                  zprofen(ji,jj,jk) = zprofmax * xnanofer(ji,jj,jk)    &
+                  &          * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn  &
                   &          + xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) ) * (1. - xnanofer(ji,jj,jk) ) )
                ENDIF
@@ -274,28 +286,32 @@
                   !  production terms for picophyto.
                   zprorcap(ji,jj,jk) = zprpic(ji,jj,jk)  * xlimpic(ji,jj,jk) * trb(ji,jj,jk,jppic) * rfact2
-                  !
+                  ! Size computation
+                  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
                   zration = trb(ji,jj,jk,jpnpi) / ( trb(ji,jj,jk,jppic) + rtrn )
                   zratiop = trb(ji,jj,jk,jpppi) / ( trb(ji,jj,jk,jppic) + rtrn )
                   zratiof = trb(ji,jj,jk,jppfe) / ( trb(ji,jj,jk,jppic) + rtrn )
-                  zprnutmax = zprnut(ji,jj,jk) * fvpuptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jppic) * rfact2
+                  zprnutmax = zprnutp(ji,jj,jk) * fvpuptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jppic) * rfact2
                   ! Uptake of nitrogen
-                  zrat = MIN( 1., zration / (xqnpmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
+                  zrat = 1.0 - MIN( 1., zration / (xqnpmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zrat**2/ (0.05**2 + zrat**2) ) )
                   zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqppmin(ji,jj,jk) )   &
                   &          / ( xqppmax(ji,jj,jk) - xqppmin(ji,jj,jk) + rtrn ), xlimpfe(ji,jj,jk) ) )
-                  zpronewp(ji,jj,jk) = zpronmax * zdaylen(ji,jj) * xpicono3(ji,jj,jk) 
+                  zpronewp(ji,jj,jk) = zpronmax * xpicono3(ji,jj,jk) 
                   zproregp(ji,jj,jk) = zpronmax * xpiconh4(ji,jj,jk)
                   ! Uptake of phosphorus
-                  zrat = MIN( 1., zratiop / (xqppmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
-                  zpropmax = zprnutmax * zmax * xlimpfe(ji,jj,jk)
+                  zrat = 1.0 - MIN( 1., zratiop / (xqppmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
+                  zpropmax = zprnutmax * zmax * xlimpfe(ji,jj,jk) * 16./10.
                   zpropo4p(ji,jj,jk) = zpropmax * xpicopo4(ji,jj,jk)
                   zprodopp(ji,jj,jk) = zpropmax * xpicodop(ji,jj,jk)
                   ! Uptake of iron
-                  zrat = MIN( 1., zratiof / qfpmax )
-                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
+                  zrat = 1.0 - MIN( 1., zratiof / qfpmax )
+                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
                   zprofmax = zprnutmax * qfpmax * zmax
-                  zprofep(ji,jj,jk) = zprofmax * xpicofer(ji,jj,jk) * ( 3. - 2.4 * xlimpfe(ji,jj,jk)   &
-                  &          / ( xlimpfe(ji,jj,jk) + 0.2 ) ) * (1. + 0.8 * xpicono3(ji,jj,jk) / ( rtrn   &
+                  zprofep(ji,jj,jk) = zprofmax * xpicofer(ji,jj,jk)  &
+                  &          * (1. + 0.8 * xpicono3(ji,jj,jk) / ( rtrn   &
                   &          + xpicono3(ji,jj,jk) + xpiconh4(ji,jj,jk) ) * (1. - xpicofer(ji,jj,jk) ) )
                ENDIF
@@ -311,7 +327,9 @@
                   !  production terms for diatomees
                   zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trb(ji,jj,jk,jpdia) * rfact2
-                  ! Computation of the respiration term according to pahlow 
-                  ! & oschlies (2013)
-                  !
+                  ! size computation
+                  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
                   zration = trb(ji,jj,jk,jpndi) / ( trb(ji,jj,jk,jpdia) + rtrn )
                   zratiop = trb(ji,jj,jk,jppdi) / ( trb(ji,jj,jk,jpdia) + rtrn )
@@ -319,22 +337,22 @@
                   zprnutmax = zprnut(ji,jj,jk) * fvduptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jpdia) * rfact2
                   ! Uptake of nitrogen
-                  zrat = MIN( 1., zration / (xqndmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05)) 
+                  zrat = 1.0 - MIN( 1., zration / (xqndmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
                   zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqpdmin(ji,jj,jk) )   &
                   &          / ( xqpdmax(ji,jj,jk) - xqpdmin(ji,jj,jk) + rtrn ), xlimdfe(ji,jj,jk) ) )
-                  zpronewd(ji,jj,jk) = zpronmax * zdaylen(ji,jj) * xdiatno3(ji,jj,jk)
+                  zpronewd(ji,jj,jk) = zpronmax * xdiatno3(ji,jj,jk)
                   zproregd(ji,jj,jk) = zpronmax * xdiatnh4(ji,jj,jk)
                   ! Uptake of phosphorus
-                  zrat = MIN( 1., zratiop / (xqpdmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05)) 
-                  zpropmax = zprnutmax * zmax * xlimdfe(ji,jj,jk)
+                  zrat = 1.0 - MIN( 1., zratiop / (xqpdmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zrat**2/ (0.05**2 + zrat**2) ) )
+                  zpropmax = zprnutmax * zmax * xlimdfe(ji,jj,jk) * 16./10.
                   zpropo4d(ji,jj,jk) = zpropmax * xdiatpo4(ji,jj,jk)
                   zprodopd(ji,jj,jk) = zpropmax * xdiatdop(ji,jj,jk)
                   ! Uptake of iron
-                  zrat = MIN( 1., zratiof / qfdmax )
-                  zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
+                  zrat = 1.0 - MIN( 1., zratiof / qfdmax )
+                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
                   zprofmax = zprnutmax * qfdmax * zmax
-                  zprofed(ji,jj,jk) = zprofmax * xdiatfer(ji,jj,jk) * ( 3. - 2.4 * xlimdfe(ji,jj,jk)     &
-                  &          / ( xlimdfe(ji,jj,jk) + 0.2 ) ) * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn   &
+                  zprofed(ji,jj,jk) = zprofmax * xdiatfer(ji,jj,jk)    &
+                  &          * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn   &
                   &          + xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) ) * (1. - xdiatfer(ji,jj,jk) ) )
                ENDIF
@@ -350,21 +368,15 @@
                   znanotot = enanom(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
                   zprod = rday * (zpronewn(ji,jj,jk) + zproregn(ji,jj,jk)) * zprchln(ji,jj,jk) * xlimphy(ji,jj,jk)
-                  thetannm_n   = MIN ( thetannm, ( thetannm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
-                  &               * (1. - 1.14 / 43.4 * 20.))
-                  zprochln = thetannm_n * zprod / ( zpislopeadn(ji,jj,jk) * znanotot + rtrn )
+                  zprochln = thetannm * zprod / ( zpislopeadn(ji,jj,jk) * znanotot + rtrn )
                   zprochln = MAX(zprochln, chlcmin * 12. * zprorcan (ji,jj,jk) )
                      !  production terms for picophyto. ( chlorophyll )
                   zpicotot = epicom(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
                   zprod = rday * (zpronewp(ji,jj,jk) + zproregp(ji,jj,jk)) * zprchlp(ji,jj,jk) * xlimpic(ji,jj,jk)
-                  thetanpm_n   = MIN ( thetanpm, ( thetanpm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
-                  &               * (1. - 1.14 / 43.4 * 20.))
-                  zprochlp = thetanpm_n * zprod / ( zpislopeadp(ji,jj,jk) * zpicotot + rtrn )
+                  zprochlp = thetanpm * zprod / ( zpislopeadp(ji,jj,jk) * zpicotot + rtrn )
                   zprochlp = MAX(zprochlp, chlcmin * 12. * zprorcap(ji,jj,jk) )
                   !  production terms for diatomees ( 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)
-                  thetandm_n   = MIN ( thetandm, ( thetandm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
-                  &               * (1. - 1.14 / 43.4 * 20.))
-                  zprochld = thetandm_n * zprod / ( zpislopeadd(ji,jj,jk) * zdiattot + rtrn )
+                  zprochld = thetandm * zprod / ( zpislopeadd(ji,jj,jk) * zdiattot + rtrn )
                   zprochld = MAX(zprochld, chlcmin * 12. * zprorcad(ji,jj,jk) )
                   !   Update the arrays TRA which contain the Chla sources and sinks
@@ -416,5 +428,5 @@
               &                     + zprodopd(ji,jj,jk) * texcretd
               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcretd
-              tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcretd
+              tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk) * rfact2 * trb(ji,jj,jk,jpdia)
               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)  &
               &                     + excretp * zprorcap(ji,jj,jk)
@@ -430,5 +442,5 @@
               zfeup = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) + texcretp * zprofep(ji,jj,jk)
               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zfeup
-              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - texcretd * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
+              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk) * rfact2 * trb(ji,jj,jk,jpdia)
               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk) - zprorcap(ji,jj,jk)  &
               &                     + zpsino3 * zpronewn(ji,jj,jk) + zpsinh4 * zproregn(ji,jj,jk)   &
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 12348)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/sms_pisces.F90	(revision 12349)
@@ -103,8 +103,10 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   consgoc    !: Calcite production
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   blim       !: bacterial production factor
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizen      !: size of diatoms 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizep      !: size of diatoms 
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizen      !: size of nanophyto
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizep      !: size of picophyto
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sized      !: size of diatoms 
-
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizena     !: size of nanophytoplankton, after
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizepa     !: size of picophyto, after
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizeda     !: size of diatomss, after
 
    !!* Variable for chemistry of the CO2 cycle
@@ -120,4 +122,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tgfunc    !: Temp. dependancy of various biological rates
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tgfunc2   !: Temp. dependancy of mesozooplankton rates
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tgfunc3   !: Temp. dependancy of picophytoplankton rates
 
    !!----------------------------------------------------------------------
@@ -165,5 +168,6 @@
          !
          !* Temperature dependancy of SMS terms
-         ALLOCATE( tgfunc(jpi,jpj,jpk)  , tgfunc2(jpi,jpj,jpk),   STAT=ierr(6) )
+         ALLOCATE( tgfunc (jpi,jpj,jpk) , tgfunc2(jpi,jpj,jpk),     &
+            &      tgfunc3(jpi,jpj,jpk) , STAT=ierr(6) )
          !
          !* Sinkong speed
@@ -182,5 +186,6 @@
          !*  Size of phytoplankton cells
          ALLOCATE( sizen(jpi,jpj,jpk), sizep(jpi,jpj,jpk),         &
-           &       sized(jpi,jpj,jpk),                            STAT=ierr(10) )
+           &       sized(jpi,jpj,jpk), sizena(jpi,jpj,jpk),        &
+           &       sizepa(jpi,jpj,jpk), sizeda(jpi,jpj,jpk),       STAT=ierr(10) )
       ENDIF
       !
