Index: NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/namelist_pisces_ref
===================================================================
--- NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/namelist_pisces_ref	(revision 12682)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/namelist_pisces_ref	(revision 12759)
@@ -69,11 +69,11 @@
    concdfer   =  3.E-9    ! Iron half saturation for diatoms
    concbfe    =  1.E-11   ! Iron half-saturation for DOC remin.
-   concbnh4   =  2.E-8    ! NH4 half saturation for DOC remin.
+   concbnh4   =  6.E-8    ! NH4 half saturation for DOC remin.
    concbno3   =  2.E-7    ! Nitrate half saturation for DOC remin.
    xsizedia   =  1.E-6    ! Minimum size criteria for diatoms
    xsizephy   =  1.E-6    ! Minimum size criteria for phyto
    xsizern    =  3.0      ! Size ratio for nanophytoplankton
-   xsizerd    =  3.0      ! Size ratio for diatoms
-   xksi1      =  4.E-6    ! half saturation constant for Si uptake
+   xsizerd    =  4.0      ! Size ratio for diatoms
+   xksi1      =  8.E-6    ! half saturation constant for Si uptake
    xksi2      =  20E-6    ! half saturation constant for Si/C
    xkdoc      =  417.E-6  ! half-saturation constant of DOC remineralization
@@ -160,5 +160,5 @@
    fecnm      =  40E-6    ! Maximum Fe/C in nanophytoplankton
    fecdm      =  40E-6    ! Maximum Fe/C in diatoms
-   grosip     =  0.13     ! mean Si/C ratio
+   grosip     =  0.11     ! mean Si/C ratio
 /
 !-----------------------------------------------------------------------
@@ -203,7 +203,7 @@
 !-----------------------------------------------------------------------
    part2      =  0.75     ! part of calcite not dissolved in mesozoo guts
-   grazrat2   =  0.75     ! maximal mesozoo grazing rate
-   resrat2    =  0.005    ! exsudation rate of mesozooplankton
-   mzrat2     =  0.03     ! mesozooplankton mortality rate
+   grazrat2   =  0.55     ! maximal mesozoo grazing rate
+   resrat2    =  0.005    ! Linear mortality rate of mesozooplankton
+   mzrat2     =  0.025    ! mesozooplankton mortality rate
    xpref2d    =  1.       ! mesozoo preference for diatoms
    xpref2n    =  0.3      ! mesozoo preference for nanophyto.
@@ -230,5 +230,5 @@
    grazrat2    = 0.75    ! maximal mesozoo grazing rate
    bmetexc2    = .true.  ! Metabolic use of excess carbon 
-   resrat2     = 0.005   ! exsudation rate of mesozooplankton
+   resrat2     = 0.005   ! Linear mortality rate of mesozooplankton
    mzrat2      = 0.03    ! mesozooplankton mortality rate
    xpref2d     = 1.      ! meso preference for diatoms
@@ -259,10 +259,10 @@
 !-----------------------------------------------------------------------
    part       =  0.5      ! part of calcite not dissolved in microzoo guts
-   grazrat    =  3.0      ! maximal zoo grazing rate
-   resrat     =  0.03     ! exsudation rate of zooplankton
-   mzrat      =  0.004    ! zooplankton mortality rate
+   grazrat    =  2.0      ! maximal zoo grazing rate
+   resrat     =  0.03     ! Linear mortality rate of zooplankton
+   mzrat      =  0.005    ! zooplankton mortality rate
    xprefc     =  0.1      ! Microzoo preference for POM
    xprefn     =  1.       ! Microzoo preference for Nanophyto
-   xprefd     =  0.6      ! Microzoo preference for Diatoms
+   xprefd     =  0.8      ! Microzoo preference for Diatoms
    xthreshdia =  1.E-8    ! Diatoms feeding threshold for microzooplankton 
    xthreshphy =  1.E-8    ! Nanophyto feeding threshold for microzooplankton 
@@ -270,6 +270,6 @@
    xthresh    =  3.E-7    ! Food threshold for feeding
    xkgraz     =  20.E-6   ! half sturation constant for grazing
-   epsher     =  0.3      ! Efficiency of microzoo growth
-   epshermin  =  0.3      ! Minimum efficiency of microzoo growth
+   epsher     =  0.35     ! Efficiency of microzoo growth
+   epshermin  =  0.35     ! Minimum efficiency of microzoo growth
    sigma1     =  0.6      ! Fraction of microzoo excretion as DOM
    unass      =  0.3      ! non assimilated fraction of phyto by zoo
@@ -281,5 +281,5 @@
    grazrat    =  3.0      ! maximal zoo grazing rate
    bmetexc    =  .true.   ! Metabolic use of excess carbon
-   resrat     =  0.03     ! exsudation rate of zooplankton
+   resrat     =  0.03     ! Linear mortality rate of zooplankton
    mzrat      =  0.004    ! zooplankton mortality rate
    xprefc     =  0.1      ! Microzoo preference for POM
@@ -332,5 +332,5 @@
 &nampispoc     !   parameters for organic particles
 !-----------------------------------------------------------------------
-   xremip    =  0.035     ! remineralisation rate of POC
+   xremip    =  0.03      ! remineralisation rate of POC
    jcpoc     =  15        ! Number of lability classes
    rshape    =  1.0       ! Shape of the gamma function
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 12682)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 12759)
@@ -107,7 +107,7 @@
 
       ! Update of the size of the different phytoplankton groups
+      sized(:,:,:) = sizeda(:,:,:)
+      sizen(:,:,:) = sizena(:,:,:)
       IF (ln_p5z) THEN
-         sized(:,:,:) = sizeda(:,:,:)
-         sizen(:,:,:) = sizena(:,:,:)
          sizep(:,:,:) = sizepa(:,:,:)
       ENDIF
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 12682)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlim.F90	(revision 12759)
@@ -61,4 +61,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concdfe    !: Limitation of diatoms uptake of Fe
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concnfe    !: Limitation of Nano uptake of Fe
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanofer   !: Limitation of Fe uptake by nanophyto
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatfer   !: Limitation of Fe uptake by diatoms
 
    ! Coefficient for iron limitation following Flynn and Hipkin (1999)
@@ -86,12 +88,14 @@
       !
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zno3, zferlim
-      REAL(wp) ::   zconcd, zconcd2, zconcn, zconcn2
+      REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zno3, zferlim, zcoef
       REAL(wp) ::   z1_trbdia, z1_trbphy, ztem1, ztem2, zetot1, zetot2
       REAL(wp) ::   zdenom, zratio, zironmin
       REAL(wp) ::   zconc1d, zconc1dnh4, zconc0n, zconc0nnh4   
+      REAL(wp) ::   fananof, fadiatf, znutlim, zfalim
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_lim')
+      !
+      sizena(:,:,:) = 0.0  ;  sizeda(:,:,:) = 0.0
       !
       DO jk = 1, jpkm1
@@ -110,18 +114,27 @@
                ! that increasing biomass is made of generally bigger cells
                !------------------------------------------------------------
-               zconcd   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
-               zconcd2  = trb(ji,jj,jk,jpdia) - zconcd
-               zconcn   = MAX( 0.e0 , trb(ji,jj,jk,jpphy) - xsizephy )
-               zconcn2  = trb(ji,jj,jk,jpphy) - zconcn
                z1_trbphy   = 1. / ( trb(ji,jj,jk,jpphy) + rtrn )
                z1_trbdia   = 1. / ( trb(ji,jj,jk,jpdia) + rtrn )
 
-               concdfe(ji,jj,jk) = MAX( concdfer, ( zconcd2 * concdfer + concdfer * xsizerd * zconcd ) * z1_trbdia )
-               zconc1d           = MAX( concdno3, ( zconcd2 * concdno3 + concdno3 * xsizerd * zconcd ) * z1_trbdia )
-               zconc1dnh4        = MAX( concdnh4, ( zconcd2 * concdnh4 + concdnh4 * xsizerd * zconcd ) * z1_trbdia )
-
-               concnfe(ji,jj,jk) = MAX( concnfer, ( zconcn2 * concnfer + concnfer * xsizern * zconcn ) * z1_trbphy )
-               zconc0n           = MAX( concnno3, ( zconcn2 * concnno3 + concnno3 * xsizern * zconcn ) * z1_trbphy )
-               zconc0nnh4        = MAX( concnnh4, ( zconcn2 * concnnh4 + concnnh4 * xsizern * zconcn ) * z1_trbphy )
+               concnfe(ji,jj,jk) = concnfer * sizen(ji,jj,jk)**0.81
+               zconc0n           = concnno3 * sizen(ji,jj,jk)**0.81
+               zconc0nnh4        = concnnh4 * sizen(ji,jj,jk)**0.81
+
+               concdfe(ji,jj,jk) = concdfer * sized(ji,jj,jk)**0.81 
+               zconc1d           = concdno3 * sized(ji,jj,jk)**0.81 
+               zconc1dnh4        = concdnh4 * sized(ji,jj,jk)**0.81  
+
+               ! Computation of the optimal allocation parameters
+               ! Based on the different papers by Pahlow et al., and 
+               ! Smith et al.
+               ! ---------------------------------------------------
+
+               ! Nanophytoplankton
+               znutlim = biron(ji,jj,jk) / concnfe(ji,jj,jk)
+               fananof = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+
+               ! Diatoms
+               znutlim = biron(ji,jj,jk) / concdfe(ji,jj,jk)
+               fadiatf = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
 
                ! Michaelis-Menten Limitation term by nutrients of
@@ -144,4 +157,9 @@
                ! Michaelis-Menten Limitation term by nutrients: Nanophyto
                ! --------------------------------------------------------
+               ! Limitation of Fe uptake
+               zfalim = (1.-fananof) / fananof
+               xnanofer(ji,jj,jk) = (1. - fananof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * concnfe(ji,jj,jk) )
+
+               ! Limitation of nanophytoplankton growth
                zdenom = 1. /  ( zconc0n * zconc0nnh4 + zconc0nnh4 * trb(ji,jj,jk,jpno3) + zconc0n * trb(ji,jj,jk,jpnh4) )
                xnanono3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * zconc0nnh4 * zdenom
@@ -163,4 +181,9 @@
                !   Michaelis-Menten Limitation term by nutrients : Diatoms
                !   -------------------------------------------------------
+               ! Limitation of Fe uptake
+               zfalim = (1.-fadiatf) / fadiatf
+               xdiatfer(ji,jj,jk) = (1. - fadiatf) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * concdfe(ji,jj,jk) )
+
+               ! Limitation of diatoms growth
                zdenom   = 1. / ( zconc1d * zconc1dnh4 + zconc1dnh4 * trb(ji,jj,jk,jpno3) + zconc1d * trb(ji,jj,jk,jpnh4) )
                xdiatno3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * zconc1dnh4 * zdenom
@@ -184,4 +207,20 @@
          END DO
       END DO
+
+      ! Size estimation of phytoplankton based on total biomass
+      ! Assumes that larger biomass implies addition of larger cells
+      ! ------------------------------------------------------------
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zcoef = trb(ji,jj,jk,jpphy) - MIN(xsizephy, trb(ji,jj,jk,jpphy) )
+               sizena(ji,jj,jk) = 1. + ( xsizern -1.0 ) * zcoef / ( xsizephy + zcoef )
+               zcoef = trb(ji,jj,jk,jpdia) - MIN(xsizedia, trb(ji,jj,jk,jpdia) )
+               sizeda(ji,jj,jk) = 1. + ( xsizerd - 1.0 ) * zcoef / ( xsizedia + zcoef )
+
+            END DO
+         END DO
+      END DO
+
 
       ! Compute the fraction of nanophytoplankton that is made of calcifiers
@@ -235,4 +274,6 @@
         IF( iom_use( "LNFe"    ) )   CALL iom_put( "LNFe"   , xlimnfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
         IF( iom_use( "LDFe"    ) )   CALL iom_put( "LDFe"   , xlimdfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "SIZEN"   ) )   CALL iom_put( "SIZEN"  , sizen(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "SIZED"   ) )   CALL iom_put( "SIZED"  , sized(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
       ENDIF
       !
@@ -315,4 +356,5 @@
          &      xnanonh4(jpi,jpj,jpk), xdiatnh4(jpi,jpj,jpk),       &
          &      xnanopo4(jpi,jpj,jpk), xdiatpo4(jpi,jpj,jpk),       &
+         &      xnanofer(jpi,jpj,jpk), xdiatfer(jpi,jpj,jpk),       &
          &      xlimphy (jpi,jpj,jpk), xlimdia (jpi,jpj,jpk),       &
          &      xlimnfe (jpi,jpj,jpk), xlimdfe (jpi,jpj,jpk),       &
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 12682)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 12759)
@@ -79,7 +79,8 @@
       REAL(wp) :: zmortzgoc, zfrac, zfracfe, zratio, zratio2, zfracal, zgrazcal
       REAL(wp) :: zepsherf, zepshert, zepsherq, zepsherv, zgrarsig, zgraztotc, zgraztotn, zgraztotf
-      REAL(wp) :: zmigreltime, zprcaca, zmortz, zgrasrat, zgrasratn
+      REAL(wp) :: zmigreltime, zprcaca, zmortz, zgrasratf, zgrasratn
       REAL(wp) :: zrespz, ztortz, zgrazd, zgrazz, zgrazpof, zgrazn, zgrazpoc, zgraznf, zgrazf
       REAL(wp) :: zgrazfffp, zgrazfffg, zgrazffep, zgrazffeg, zrum, zcodel, zargu, zval
+      REAL(wp) :: zsigma, zdiffdn, ztmp1, ztmp2, ztmp3, ztmp4, ztmptot
       CHARACTER (len=25) :: charout
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing, zfezoo2
@@ -153,12 +154,38 @@
                zgraze2   = grazrat2 * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpmes) * (1. - nitrfac(ji,jj,jk)) 
 
-               ! The grazing pressure on each prey is computed assuming passive switching. This 
-               ! is equivalent to assuming that mesozooplankton have an opportunistic feeding
-               ! behaviour.
-               ! -----------------------------------------------------------------------------
-               zgrazd    = zgraze2  * xpref2d  * zcompadi  * zdenom2 
-               zgrazz    = zgraze2  * xpref2z  * zcompaz   * zdenom2 
-               zgrazn    = zgraze2  * xpref2n  * zcompaph  * zdenom2 
-               zgrazpoc  = zgraze2  * xpref2c  * zcompapoc * zdenom2 
+               ! An active switching parameterization is used here.
+               ! We don't use the KTW parameterization proposed by 
+               ! Vallina et al. because it tends to produce too steady biomass
+               ! composition and the variance of Chl is too low as it grazes
+               ! too strongly on winning organisms. We use a generalized
+               ! switching parameterization proposed by Morozov and 
+               ! Petrovskii (2013)
+               ! ------------------------------------------------------------  
+               ! The width of the selection window is increased when preys
+               ! have low abundance, .i.e. zooplankton become less specific 
+               ! to avoid starvation.
+               ! ----------------------------------------------------------
+               zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
+               zsigma = 0.5 + 1.0 * zsigma
+               ! Nanophytoplankton and diatoms are the only preys considered
+               ! to be close enough to have potential interference
+               ! -----------------------------------------------------------
+               zdiffdn = exp( -ABS(log(1.5 * sizen(ji,jj,jk) / (5.0 * sized(ji,jj,jk) + rtrn )) )**2 / zsigma**2 )
+               ztmp1 = xpref2n * zcompaph * ( zcompaph + zdiffdn * zcompadi ) / ( 1.0 + zdiffdn )
+               ztmp2 = xpref2c * zcompapoc**2
+               ztmp3 = xpref2d * zcompadi * ( zdiffdn * zcompadi + zcompaph ) / ( 1.0 + zdiffdn )
+               ztmp4 = xpref2z * zcompaz**2
+               ztmptot = ztmp1 + ztmp2 + ztmp3 + ztmp4 + rtrn
+               ztmp1 = ztmp1 / ztmptot
+               ztmp2 = ztmp2 / ztmptot
+               ztmp3 = ztmp3 / ztmptot
+               ztmp4 = ztmp4 / ztmptot
+
+               !   Mesozooplankton regular grazing on the different preys
+               !   ------------------------------------------------------
+               zgrazd    = zgraze2  * ztmp3 * zdenom
+               zgrazn    = zgraze2  * ztmp1 * zdenom
+               zgrazpoc  = zgraze2  * ztmp2 * zdenom
+               zgrazz    = zgraze2  * ztmp4 * zdenom
 
                zgraznf   = zgrazn   * trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) + rtrn)
@@ -215,7 +242,7 @@
                ! Fulton, 2012)
                ! -----------------------------------------------------------------------------------
-               zgrasrat  =  ( zgraztotf + rtrn )/ ( zgraztotc + rtrn )
+               zgrasratf =  ( zgraztotf + rtrn )/ ( zgraztotc + rtrn )
                zgrasratn =  ( zgraztotn + rtrn )/ ( zgraztotc + rtrn )
-               zepshert  = MIN( 1., zgrasratn, zgrasrat / ferat3)
+               zepshert  = MIN( 1., zgrasratn, zgrasratf / ferat3)
                zbeta     = MAX(0., (epsher2 - epsher2min) )
                zepsherf  = epsher2min + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta ) 
@@ -263,5 +290,5 @@
                zgrarem(ji,jj,jk) = zgraztotc * ( 1. - zepsherv - unass2 ) &
                &         + ( 1. - epsher2 - unass2 ) / ( 1. - epsher2 ) * ztortz
-               zgraref(ji,jj,jk) = zgraztotc * MAX( 0. , ( 1. - unass2 ) * zgrasrat - ferat3 * zepsherv )    &
+               zgraref(ji,jj,jk) = zgraztotc * MAX( 0. , ( 1. - unass2 ) * zgrasratf - ferat3 * zepsherv )    &
                &         + ferat3 * ( ( 1. - epsher2 - unass2 ) /( 1. - epsher2 ) * ztortz )
                zgrapoc(ji,jj,jk) = zgraztotc * unass2 + zmortzgoc
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 12682)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 12759)
@@ -68,11 +68,10 @@
       INTEGER  :: ji, jj, jk
       REAL(wp) :: zcompadi, zcompaz , zcompaph, zcompapoc
-      REAL(wp) :: zgraze  , zdenom, zdenom2
-      REAL(wp) :: zfact   , zfood, zfoodlim, zbeta
+      REAL(wp) :: zgraze  , zdenom, zdenom2, zfact, zfood, zfoodlim, zbeta
       REAL(wp) :: zepsherf, zepshert, zepsherq, zepsherv, zgrarsig, zgraztotc, zgraztotn, zgraztotf
       REAL(wp) :: zgrarem, zgrafer, zgrapoc, zprcaca, zmortz
-      REAL(wp) :: zrespz, ztortz, zgrasrat, zgrasratn
-      REAL(wp) :: zgrazp, zgrazm, zgrazsd
-      REAL(wp) :: zgrazmf, zgrazsf, zgrazpf
+      REAL(wp) :: zrespz, ztortz, zgrasratf, zgrasratn
+      REAL(wp) :: zgraznc, zgrazpoc, zgrazdc, zgrazpof, zgrazdf, zgraznf
+      REAL(wp) :: zsigma, zdiffdn, ztmp1, ztmp2, ztmp3, ztmptot, zproport
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing, zfezoo
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d, zzligprod
@@ -92,4 +91,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 ))
 
                !  linear mortality of mesozooplankton
@@ -113,5 +116,5 @@
                !   big for microzooplankton.
                !   --------------------------------------------------------
-               zcompadi  = MIN( MAX( ( trb(ji,jj,jk,jpdia) - xthreshdia ), 0.e0 ), xsizedia )
+               zcompadi  = zproport * MAX( ( trb(ji,jj,jk,jpdia) - xthreshdia ), 0.e0 )
                zcompaph  = MAX( ( trb(ji,jj,jk,jpphy) - xthreshphy ), 0.e0 )
                zcompapoc = MAX( ( trb(ji,jj,jk,jppoc) - xthreshpoc ), 0.e0 )
@@ -130,19 +133,38 @@
                zgraze    = grazrat * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpzoo) * (1. - nitrfac(ji,jj,jk))
 
-               ! The grazing pressure on each prey is computed assuming passive switching. This 
-               ! is equivalent to assuming that microzooplankton have an opportunistic feeding
-               ! behaviour.
-               ! -----------------------------------------------------------------------------
-               zgrazp    = zgraze  * xprefn * zcompaph  * zdenom2 
-               zgrazm    = zgraze  * xprefc * zcompapoc * zdenom2 
-               zgrazsd   = zgraze  * xprefd * zcompadi  * zdenom2 
-
-               zgrazpf   = zgrazp  * trb(ji,jj,jk,jpnfe) / (trb(ji,jj,jk,jpphy) + rtrn)
-               zgrazmf   = zgrazm  * trb(ji,jj,jk,jpsfe) / (trb(ji,jj,jk,jppoc) + rtrn)
-               zgrazsf   = zgrazsd * trb(ji,jj,jk,jpdfe) / (trb(ji,jj,jk,jpdia) + rtrn)
-               !
-               zgraztotc = zgrazp  + zgrazm  + zgrazsd 
-               zgraztotf = zgrazpf + zgrazsf + zgrazmf 
-               zgraztotn = zgrazp * quotan(ji,jj,jk) + zgrazm + zgrazsd * quotad(ji,jj,jk)
+               ! An active switching parameterization is used here.
+               ! We don't use the KTW parameterization proposed by 
+               ! Vallina et al. because it tends to produce too steady biomass
+               ! composition and the variance of Chl is too low as it grazes
+               ! too strongly on winning organisms. We use a generalized
+               ! switching parameterization proposed by Morozov and 
+               ! Petrovskii (2013)
+               ! ------------------------------------------------------------  
+               ! The width of the selection window is increased when preys
+               ! have low abundance, .i.e. zooplankton become less specific 
+               ! to avoid starvation.
+               ! ----------------------------------------------------------
+               zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
+               zsigma = 0.5 + 1.0*zsigma
+               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 )
+               ztmp2 = xprefd * zcompadi * ( zdiffdn * zcompaph + zcompadi ) / ( 1.0 + zdiffdn )
+               ztmp3 = xprefc * zcompapoc**2
+               ztmptot = ztmp1 + ztmp2 + ztmp3 + rtrn
+               ztmp1 = ztmp1 / ztmptot
+               ztmp2 = ztmp2 / ztmptot
+               ztmp3 = ztmp3 / ztmptot
+
+               zgraznc   = zgraze   * ztmp1 * zdenom
+               zgrazdc   = zgraze   * ztmp2 * zdenom
+               zgrazpoc  = zgraze   * ztmp3 * zdenom
+
+               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)
+               !
+               zgraztotc = zgraznc + zgrazpoc + zgrazdc
+               zgraztotf = zgraznf + zgrazdf  + zgrazpof 
+               zgraztotn = zgraznc * quotan(ji,jj,jk) + zgrazpoc + zgrazdc * quotad(ji,jj,jk)
 
                ! Grazing by microzooplankton
@@ -157,12 +179,12 @@
                ! Fulton, 2012)
                ! -----------------------------------------------------------------------------
-               zgrasrat  = ( zgraztotf + rtrn ) / ( zgraztotc + rtrn )
+               zgrasratf = ( zgraztotf + rtrn ) / ( zgraztotc + rtrn )
                zgrasratn = ( zgraztotn + rtrn ) / ( zgraztotc + rtrn )
-               zepshert  =  MIN( 1., zgrasratn, zgrasrat / ferat3)
+               zepshert  =  MIN( 1., zgrasratn, zgrasratf / ferat3)
                zbeta     = MAX(0., (epsher - epshermin) )
                zepsherf  = epshermin + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
                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 ) 
+               zgrafer   = zgraztotc * MAX( 0. , ( 1. - unass ) * zgrasratf - ferat3 * zepsherv ) 
                zgrarem   = zgraztotc * ( 1. - zepsherv - unass )
                zgrapoc   = zgraztotc * unass
@@ -190,20 +212,20 @@
                zmortz = ztortz + zrespz
                tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) - zmortz + zepsherv * zgraztotc 
-               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgrazp
-               tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazsd
-               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazp  * trb(ji,jj,jk,jpnch)/(trb(ji,jj,jk,jpphy)+rtrn)
-               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazsd * trb(ji,jj,jk,jpdch)/(trb(ji,jj,jk,jpdia)+rtrn)
-               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) - zgrazsd * trb(ji,jj,jk,jpdsi)/(trb(ji,jj,jk,jpdia)+rtrn)
-               tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) + zgrazsd * trb(ji,jj,jk,jpdsi)/(trb(ji,jj,jk,jpdia)+rtrn)
-               tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zgrazpf
-               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazsf
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortz - zgrazm
+               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgraznc
+               tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazdc
+               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgraznc * trb(ji,jj,jk,jpnch)/(trb(ji,jj,jk,jpphy)+rtrn)
+               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazdc * trb(ji,jj,jk,jpdch)/(trb(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) - zgrazdc * trb(ji,jj,jk,jpdsi)/(trb(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) + zgrazdc * trb(ji,jj,jk,jpdsi)/(trb(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zgraznf
+               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazdf
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortz - zgrazpoc
                prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zmortz
-               conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazm
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz - zgrazmf
+               conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazpoc
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz - zgrazpof
                !
                ! Calcite remineralization due to zooplankton activity
                ! part of the ingested calcite is dissolving in the acidic gut
-               zprcaca = xfracal(ji,jj,jk) * zgrazp
+               zprcaca = xfracal(ji,jj,jk) * zgraznc
                prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
                !
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 12682)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 12759)
@@ -68,12 +68,11 @@
       REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap, zlim, zsilfac2, zsiborn
       REAL(wp) ::   zprod, zproreg, zproreg2, zprochln, zprochld
-      REAL(wp) ::   zmaxday, zdocprod, zpislopen, zpisloped
-      REAL(wp) ::   zmxltst, zmxlday
+      REAL(wp) ::   zdocprod, zpislopen, zpisloped, zfact
+      REAL(wp) ::   zratiosi, zmaxsi, zlimfac, zsizetmp
       REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup, chlcnm_n, chlcdm_n
-      REAL(wp) ::   zfact
       CHARACTER (len=25) :: charout
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zw2d
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zstrn, zmixnano, zmixdiat
+      REAL(wp), DIMENSION(jpi,jpj    ) :: zstrn
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprmaxn,zprmaxd
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadd, zysopt  
@@ -129,5 +128,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
@@ -153,5 +152,5 @@
                   zconctemp   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
                   zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp
-                  !
+
                   ! The initial slope of the PI curve can be increased for nano
                   ! to account for photadaptation, for instance in the DCM
@@ -197,18 +196,18 @@
                 zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) )   &
                 &      * zprmaxn(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn )
-                quotan(ji,jj,jk) = MIN( 1., 0.2 + 0.8 * zval )
+                quotan(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
                 zval = MIN( xdiatpo4(ji,jj,jk), ( xdiatnh4(ji,jj,jk) + xdiatno3(ji,jj,jk) ) )   &
                 &      * zprmaxd(ji,jj,jk) / ( zprdia(ji,jj,jk) + rtrn )
-                quotad(ji,jj,jk) = MIN( 1., 0.2 + 0.8 * zval )
-            END DO
-         END DO
-      END DO
-
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+                quotad(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
+            END DO
+         END DO
+      END DO
+
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                    ! Si/C of diatoms
                    ! ------------------------
@@ -217,14 +216,20 @@
                    ! to mimic the very high ratios observed in the Southern Ocean (zsilfac2)
                    ! -----------------------------------------------------------------------
-                  zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
-                  zsilim = MIN( zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) )
-                  zsilfac = 4.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 )
-                  ELSE
-                    zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 )
-                  ENDIF
-                  zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2
+                 zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
+                 zsilim = xlimdia(ji,jj,jk) * zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn )
+                 zsiborn = trb(ji,jj,1,jpsil) * trb(ji,jj,1,jpsil) * trb(ji,jj,1,jpsil)
+                 IF (gphit(ji,jj) < -30 ) THEN
+                   zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
+                 ELSE
+                   zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 )
+                 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.7 * zmaxsi
+                 ENDIF
               ENDIF
             END DO
@@ -255,6 +260,14 @@
                   !  New production (uptake of NO3)
                   zpronewn(ji,jj,jk)  = zprorcan(ji,jj,jk)* xnanono3(ji,jj,jk) / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn )
-                  !
-                  zratio = trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * fecnm + rtrn )
+
+                  ! Size computation
+                  ! Size is made a function of the limitation of of phytoplankton growth
+                  ! Strongly limited cells are supposed to be smaller. sizena is the 
+                  ! size at time step t+1 and is thus updated at the end of the 
+                  ! current time step
+                  ! --------------------------------------------------------------------
+                  zlimfac = xlimphy(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 ) )
 
                   !  Iron uptake rates of nanophytoplankton. Upregulation  
@@ -262,9 +275,10 @@
                   !  formulation used in quota formulations. Uptake is downregulated
                   !  when the quota is close to the maximum quota 
-                  zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
+                  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) ) ) 
                   zprofen(ji,jj,jk) = fecnm * zprmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
-                  &             * ( 4. - 4.5 * xlimnfe(ji,jj,jk) / ( xlimnfe(ji,jj,jk) + 0.5 ) )    &
-                  &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concnfe(ji,jj,jk) )  &
-                  &             * zmax * trb(ji,jj,jk,jpphy) * rfact2
+                  &          * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn + xnanono3(ji,jj,jk)  &
+                  &          + 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)
                   zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trb(ji,jj,jk,jpdia) * rfact2
@@ -272,4 +286,14 @@
                   ! New production (uptake of NO3)
                   zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk) / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn )
+
+                  ! Size computation
+                  ! Size is made a function of the limitation of of phytoplankton growth
+                  ! Strongly limited cells are supposed to be smaller. sizeda is
+                  ! size at time step t+1 and is thus updated at the end of the 
+                  ! current time step. 
+                  ! --------------------------------------------------------------------
+                  zlimfac = zprchld(ji,jj,jk) * xlimdia(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 ) )
 
                   !  Iron uptake rates of nanophytoplankton. Upregulation  
@@ -277,10 +301,10 @@
                   !  formulation used in quota formulations. Uptake is downregulated
                   !  when the quota is close to the maximum quota 
-                  zratio = trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * fecdm + rtrn )
-                  zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
-                  zprofed(ji,jj,jk) = fecdm * zprmaxd(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
-                  &             * ( 4. - 4.5 * xlimdfe(ji,jj,jk) / ( xlimdfe(ji,jj,jk) + 0.5 ) )    &
-                  &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concdfe(ji,jj,jk) )  &
-                  &             * zmax * trb(ji,jj,jk,jpdia) * rfact2
+                  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) ) ) 
+                  zprofed(ji,jj,jk) = fecdm * zprmaxd(ji,jj,jk) * (1.0 - fr_i(ji,jj) )  &
+                  &          * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn + xdiatno3(ji,jj,jk)  &
+                  &          + xdiatnh4(ji,jj,jk) ) * (1. - xdiatfer(ji,jj,jk) ) )   &
+                  &          * xdiatfer(ji,jj,jk) * zmax * trb(ji,jj,jk,jpdia) * rfact2
                ENDIF
             END DO
@@ -336,5 +360,6 @@
                  tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorcad(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) + zdocprod
                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) &
@@ -343,5 +368,6 @@
                  zfeup = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(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)
                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk) ) &
@@ -401,5 +427,5 @@
           ENDIF
           IF( iom_use( "PBSi" ) )  THEN
-              zw3d(:,:,:) = zprorcad(:,:,:) * zfact * tmask(:,:,:) * zysopt(:,:,:) ! biogenic silica production
+              zw3d(:,:,:) = zprmaxd(:,:,:) * 1.E3 * tmask(:,:,:) * zysopt(:,:,:) * trb(:,:,:,jpdia) ! biogenic silica production
               CALL iom_put( "PBSi"  , zw3d )
           ENDIF
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 12682)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 12759)
@@ -320,16 +320,24 @@
          ENDIF
          !
+         ! PISCES size proxy
+         IF( iom_varid( numrtr, 'sized', ldstop = .FALSE. ) > 0 ) THEN
+            CALL iom_get( numrtr, jpdom_autoglo, 'sized' , sized(:,:,:)  )
+         ELSE
+            sized(:,:,:) = 1.
+         ENDIF
+         IF( iom_varid( numrtr, 'sizen', ldstop = .FALSE. ) > 0 ) THEN
+            CALL iom_get( numrtr, jpdom_autoglo, 'sizen' , sizen(:,:,:)  )
+         ELSE
+            sizen(:,:,:) = 1.
+         ENDIF
+
          ! PISCES-QUOTA specific part
          IF( ln_p5z ) THEN
             ! Read the size of the different phytoplankton groups
             ! If not in the restart file, they are set to 1
-            IF( iom_varid( numrtr, 'sized', ldstop = .FALSE. ) > 0 ) THEN
+            IF( iom_varid( numrtr, 'sizep', ldstop = .FALSE. ) > 0 ) THEN
                CALL iom_get( numrtr, jpdom_autoglo, 'sizep' , sizep(:,:,:)  )
-               CALL iom_get( numrtr, jpdom_autoglo, 'sizen' , sizen(:,:,:)  )
-               CALL iom_get( numrtr, jpdom_autoglo, 'sized' , sized(:,:,:)  )
             ELSE
                sizep(:,:,:) = 1.
-               sizen(:,:,:) = 1.
-               sized(:,:,:) = 1.
             ENDIF
         ENDIF
@@ -346,8 +354,8 @@
          CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) ) ! Si max concentration
          CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum ) ! Cumulative CO2 flux
+         CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sizen(:,:,:) )  ! Size of nanophytoplankton
+         CALL iom_rstput( kt, nitrst, numrtw, 'sized', sized(:,:,:) )  ! Size of diatoms
          IF( ln_p5z ) THEN
             CALL iom_rstput( kt, nitrst, numrtw, 'sizep', sizep(:,:,:) )  ! Size of picophytoplankton
-            CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sizen(:,:,:) )  ! Size of nanophytoplankton
-            CALL iom_rstput( kt, nitrst, numrtw, 'sized', sized(:,:,:) )  ! Size of diatoms
          ENDIF
       ENDIF
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 12682)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 12759)
@@ -79,7 +79,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicodop   !: Limitation of DOP uptake by picophyto
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatdop   !: Limitation of DOP uptake by diatoms
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanofer   !: Limitation of Fe uptake by nanophyto
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicofer   !: Limitation of Fe uptake by picophyto
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatfer   !: Limitation of Fe uptake by diatoms
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimpic    !: Limitation of picophyto PP by nutrients
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimpics   !: Limitation of picophyto PP by nutrients
@@ -440,5 +438,4 @@
                ! ------------------------------------------------------------
                zcoef = trb(ji,jj,jk,jpdia) - MIN(xsizedia, trb(ji,jj,jk,jpdia) )
-               sized(ji,jj,jk) = 1. + ( xsizerd - 1.0 ) * zcoef / ( xsizedia + zcoef )
                sizeda(ji,jj,jk) = 1. + ( xsizerd - 1.0 ) * zcoef / ( xsizedia + zcoef )
 
@@ -639,5 +636,4 @@
          &      xpicopo4(jpi,jpj,jpk), xpicodop(jpi,jpj,jpk),       &
          &      xnanodop(jpi,jpj,jpk), xdiatdop(jpi,jpj,jpk),       &
-         &      xnanofer(jpi,jpj,jpk), xdiatfer(jpi,jpj,jpk),       &
          &      xpicofer(jpi,jpj,jpk), xlimpfe (jpi,jpj,jpk),       &
          &      fvnuptk (jpi,jpj,jpk), fvduptk (jpi,jpj,jpk),       &
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 12682)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 12759)
@@ -194,8 +194,8 @@
                ! -----------------------------------------------------------
                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 )
+               ztmp1 = xpref2n * zcompaph * ( zcompaph + zdiffdn * zcompadi ) / (1.0 + zdiffdn)
                ztmp2 = xpref2m * zcompames**2
                ztmp3 = xpref2c * zcompapoc**2
-               ztmp4 = xpref2d * zcompadi * ( zdiffdn * zcompadi + zcompaph )
+               ztmp4 = xpref2d * zcompadi * ( zdiffdn * zcompadi + zcompaph ) / (1.0 + zdiffdn)
                ztmp5 = xpref2z * zcompaz**2
                ztmptot = ztmp1 + ztmp2 + ztmp3 + ztmp4 + ztmp5 + rtrn
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 12682)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/sms_pisces.F90	(revision 12759)
@@ -132,5 +132,5 @@
       !!----------------------------------------------------------------------
       USE lib_mpp , ONLY: ctl_stop
-      INTEGER ::   ierr(10)        ! Local variables
+      INTEGER ::   ierr(11)        ! Local variables
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -165,13 +165,16 @@
          !* Temperature dependency of SMS terms
          ALLOCATE( tgfunc (jpi,jpj,jpk) , tgfunc2(jpi,jpj,jpk),     &
-            &      tgfunc3(jpi,jpj,jpk) , STAT=ierr(6) )
+            &      tgfunc3(jpi,jpj,jpk) ,                         STAT=ierr(6) )
          !
          !* Sinking speed
-         ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4 (jpi,jpj,jpk),     &
-            &                             STAT=ierr(7) )   
+         ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4 (jpi,jpj,jpk),   STAT=ierr(7) )
+
+         !*  Size of phytoplankton cells
+         ALLOCATE( sizen (jpi,jpj,jpk), sized (jpi,jpj,jpk),        &
+           &       sizena(jpi,jpj,jpk), sizeda(jpi,jpj,jpk),      STAT=ierr(8) )
          ! 
          !* Prognostic ligand
          IF( ln_ligand ) THEN
-           ALLOCATE( plig(jpi,jpj,jpk)  ,                         STAT=ierr(8) )
+           ALLOCATE( plig(jpi,jpj,jpk)  ,                         STAT=ierr(9) )
          ENDIF
       ENDIF
@@ -179,10 +182,8 @@
       IF( ln_p5z ) THEN
          ! PISCES-QUOTA specific part      
-         ALLOCATE( epico(jpi,jpj,jpk)   , epicom(jpi,jpj,jpk) ,   STAT=ierr(9) ) 
+         ALLOCATE( epico(jpi,jpj,jpk)   , epicom(jpi,jpj,jpk) ,   STAT=ierr(10) ) 
 
          !*  Size of phytoplankton cells
-         ALLOCATE( sizen(jpi,jpj,jpk), sizep(jpi,jpj,jpk),         &
-           &       sized(jpi,jpj,jpk), sizena(jpi,jpj,jpk),        &
-           &       sizepa(jpi,jpj,jpk), sizeda(jpi,jpj,jpk),       STAT=ierr(10) )
+         ALLOCATE( sizep(jpi,jpj,jpk), sizepa(jpi,jpj,jpk),       STAT=ierr(11) )
       ENDIF
       !
