Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/field_def_nemo-pisces.xml
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/field_def_nemo-pisces.xml	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/field_def_nemo-pisces.xml	(revision 14276)
@@ -112,4 +112,5 @@
        <field id="DOM_e3t"     long_name="DOM * e3t"                                unit="mmol-N/m2"  > DOM * e3t </field >
 
+
      </field_group>
 
@@ -136,4 +137,5 @@
        <field id="SedpH"           long_name="PH"                                       unit="1"          />
        <field id="SedCO3por"       long_name="Bicarbonates"                             unit="mol/m3" />
+       <field id="Sedligand"       long_name="Ligands Concentration"                    unit="mol/m3" />
      </field_group>
 
@@ -150,5 +152,4 @@
        <field id="FlxSO4"      long_name="SO4 sediment flux"                       unit="mol/cm2/s"     />
        <field id="FlxFe2"      long_name="Fe2+ sediment flux"                      unit="mol/cm2/s"     />
-       <field id="Flxtot"      long_name="Sediment net burial rate"                unit="cm/s"     />
        <field id="dzdep"       long_name="Sedimentation rate"                      unit="cm/s"     />
        <field id="sflxclay"    long_name="Clay sedimentation rate"                 unit="g/m2/s"     />
@@ -157,4 +158,12 @@
        <field id="sflxpoc"     long_name="POC Sedimentation rate"                  unit="mol/m2/s"     />
        <field id="sflxfeo"     long_name="Fe(OH)3 Sedimentation rate"              unit="mol/m2/s"     />
+       <field id="FlxBSi"      long_name="BSi burial flux"                         unit="g/cm2/s"     />
+       <field id="FlxClay"     long_name="Clay burial flux"                        unit="g/cm2/s"     />
+       <field id="FlxPOC"      long_name="POC burial flux"                         unit="g/cm2/s"     />
+       <field id="FlxCaCO3"    long_name="Calcite burial flux"                     unit="g/cm2/s"     />
+       <field id="FlxPOS"      long_name="POS burial flux"                         unit="g/cm2/s"     />
+       <field id="FlxPOR"      long_name="POR burial flux"                         unit="g/cm2/s"     />
+       <field id="FlxFeO"      long_name="FeO burial flux"                         unit="g/cm2/s"     />
+       <field id="FlxFeS"      long_name="FeS burial flux"                         unit="g/cm2/s"     />
 </field_group>
 
@@ -206,4 +215,7 @@
        <field id="SIZEP"       long_name="Mean relative size of picophyto."        unit="-"          grid_ref="grid_T_3D" />
        <field id="SIZED"       long_name="Mean relative size of diatoms"           unit="-"          grid_ref="grid_T_3D" />
+       <field id="RASSD"       long_name="Mean relative size of diatoms"           unit="-"          grid_ref="grid_T_3D" />
+       <field id="RASSN"       long_name="Mean relative size of diatoms"           unit="-"          grid_ref="grid_T_3D" />
+       <field id="RASSP"       long_name="Mean relative size of diatoms"           unit="-"          grid_ref="grid_T_3D" />
        <field id="Fe3"         long_name="Iron III concentration"                  unit="nmol/m3"    grid_ref="grid_T_3D" />
        <field id="FeL1"        long_name="Complexed Iron concentration with L1"    unit="nmol/m3"    grid_ref="grid_T_3D" />
@@ -240,4 +252,5 @@
        <field id="BACT"        long_name="Bacterial Biomass"                       unit="mmol/m3"  grid_ref="grid_T_3D" />
        <field id="FEBACT"      long_name="Bacterial uptake of Fe"                  unit="molFe/m3/s"  grid_ref="grid_T_3D" />
+       <field id="FEPREC"      long_name="precipitation of Fe"                     unit="molFe/m3/s"  grid_ref="grid_T_3D" />
        <field id="LPRODR"      long_name="OM remineralisation ligand production rate" unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
        <field id="LPRODP"      long_name="phytoplankton ligand production rate"    unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/namelist_pisces_ref	(revision 14276)
@@ -14,6 +14,6 @@
 !-----------------------------------------------------------------------
   ln_p2z    = .false.        !  LOBSTER model used
-  ln_p4z    = .true.         !  PISCES model used
-  ln_p5z    = .false.        !  PISCES QUOTA model used
+  ln_p4z    = .false.        !  PISCES model used
+  ln_p5z    = .true.         !  PISCES QUOTA model used
   ln_ligand = .false.        !  Enable  organic ligands
   ln_sediment = .false.      !  Enable sediment module
@@ -44,18 +44,19 @@
 &nampisbio     !   biological parameters
 !-----------------------------------------------------------------------
-   nrdttrc    =  1        ! time step frequency for biology
-   wsbio      =  2.       ! POC sinking speed
-   xkmort     =  2.E-7    ! half saturation constant for mortality
-   ferat3     =  10.E-6   ! Fe/C in zooplankton 
-   wsbio2     =  50.      ! Big particles sinking speed
-   wsbio2max  =  50.      ! Big particles maximum sinking speed
-   wsbio2scale =  5000.    ! Big particles length scale of sinking
+   nrdttrc     =  1       ! time step frequency for biology
+   wsbio       =  2.      ! POC sinking speed
+   xkmort      =  2.E-7   ! half saturation constant for mortality
+   feratz      =  10.E-6  ! Fe/C in zooplankton 
+   feratm      =  15.E-6  ! Fe/C in mesozooplankton
+   wsbio2      =  50.     ! Big particles sinking speed
+   wsbio2max   =  50.     ! Big particles maximum sinking speed
+   wsbio2scale =  5000.   ! Big particles length scale of sinking
 !                         !  ln_ligand enabled
-   ldocp      =  1.E-4    ! Phyto ligand production per unit doc 
-   ldocz      =  1.E-4    ! Zoo ligand production per unit doc 
-   lthet      =  1.0      ! Proportional loss of ligands due to Fe uptake 
+   ldocp       =  1.E-4   ! Phyto ligand production per unit doc 
+   ldocz       =  1.E-4   ! Zoo ligand production per unit doc 
+   lthet       =  1.0     ! Proportional loss of ligands due to Fe uptake 
 !                         !  ln_p5z enabled
-   no3rat3    =  0.151    ! N/C ratio in zooplankton
-   po4rat3    =  0.00943  ! P/C ratio in zooplankton
+   no3rat3     =  0.151   ! N/C ratio in zooplankton
+   po4rat3     =  0.00944 ! P/C ratio in zooplankton
 /
 !-----------------------------------------------------------------------
@@ -64,11 +65,11 @@
    concnno3   =  1.e-6    ! Nitrate half saturation of nanophytoplankton
    concdno3   =  3.E-6    ! Nitrate half saturation for diatoms
-   concnnh4   =  1.E-7    ! NH4 half saturation for phyto
-   concdnh4   =  3.E-7    ! NH4 half saturation for diatoms
-   concnfer   =  1.E-9    ! Iron half saturation for phyto
-   concdfer   =  3.E-9    ! Iron half saturation for diatoms
-   concbfe    =  1.E-11   ! Iron 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.
+   concnnh4   =  1.E-6    ! NH4 half saturation for phyto
+   concdnh4   =  3.E-6    ! NH4 half saturation for diatoms
+   concnfer   =  3.E-9    ! Iron half saturation for phyto
+   concdfer   =  9.E-9    ! Iron half saturation for diatoms
+   concbfe    =  3.E-11   ! Iron half-saturation for DOC remin.
+   concbnh4   =  3.E-7    ! NH4 half saturation for DOC remin.
+   concbno3   =  3.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
@@ -78,7 +79,7 @@
    xksi2      =  20E-6    ! half saturation constant for Si/C
    xkdoc      =  417.E-6  ! half-saturation constant of DOC remineralization
-   qnfelim    =  7.E-6    ! Optimal quota of phyto
-   qdfelim    =  7.E-6    ! Optimal quota of diatoms
-   caco3r     =  0.3      ! mean rain ratio
+   qnfelim    =  10.E-6   ! Optimal quota of phyto
+   qdfelim    =  10.E-6   ! Optimal quota of diatoms
+   caco3r     =  0.28     ! mean rain ratio
    oxymin     =  1.E-6    ! Half-saturation constant for anoxia
 /
@@ -86,20 +87,20 @@
 &namp5zlim     !   parameters for nutrient limitations PISCES QUOTA    - ln_p5z
 !-----------------------------------------------------------------------
-   concnno3   =  3e-6     ! Nitrate half saturation of nanophytoplankton
-   concpno3   =  8e-7     ! Nitrate half saturation of picophytoplankton
-   concdno3   =  5E-6     ! Phosphate half saturation for diatoms
-   concnnh4   =  1E-6     ! NH4 half saturation for phyto
-   concpnh4   =  2.3E-7   ! NH4 half saturation for picophytoplankton
-   concdnh4   =  1.7E-6   ! NH4 half saturation for diatoms
-   concnpo4   =  3E-6     ! PO4 half saturation for phyto
-   concppo4   =  8E-7     ! PO4 half saturation for picophytoplankton
-   concdpo4   =  5E-6     ! PO4 half saturation for diatoms
-   concnfer   =  3E-9     ! Iron half saturation for phyto
-   concpfer   =  8E-10    ! Iron half saturation for picophytoplankton
-   concdfer   =  5E-9     ! Iron half saturation for diatoms
-   concbfe    =  1.E-11   ! Half-saturation for Fe limitation of Bacteria
-   concbnh4   =  1.E-7    ! NH4 half saturation for phyto
-   concbno3   =  5.E-7    ! Phosphate half saturation for diatoms
-   concbpo4   =  1E-7     ! Phosphate half saturation for bacteria
+   concnno3   =  2e-6     ! Nitrate half saturation of nanophytoplankton
+   concpno3   =  7e-7     ! Nitrate half saturation of picophytoplankton
+   concdno3   =  3E-6     ! Phosphate half saturation for diatoms
+   concnnh4   =  2E-6     ! NH4 half saturation for phyto
+   concpnh4   =  7E-7     ! NH4 half saturation for picophytoplankton
+   concdnh4   =  3E-6     ! NH4 half saturation for diatoms
+   concnpo4   =  2E-6     ! PO4 half saturation for phyto
+   concppo4   =  7E-7     ! PO4 half saturation for picophytoplankton
+   concdpo4   =  3E-6     ! PO4 half saturation for diatoms
+   concnfer   =  6E-9     ! Iron half saturation for phyto
+   concpfer   =  2E-9     ! Iron half saturation for picophytoplankton
+   concdfer   =  9E-9     ! Iron half saturation for diatoms
+   concbfe    =  3E-11    ! Half-saturation for Fe limitation of Bacteria
+   concbnh4   =  4.E-7    ! NH4 half saturation for phyto
+   concbno3   =  4.E-7    ! Phosphate half saturation for diatoms
+   concbpo4   =  4.E-7    ! Phosphate half saturation for bacteria
    xsizedia   =  1.E-6    ! Minimum size criteria for diatoms
    xsizephy   =  1.E-6    ! Minimum size criteria for phyto
@@ -117,22 +118,22 @@
 &namp5zquota   !   parameters for nutrient limitations PISCES quota    - ln_p5z
 !-----------------------------------------------------------------------
-   qfnopt     =  7.E-6    ! Optimal Fe quota of nanophyto
-   qfpopt     =  7.E-6    ! Optimal Fe quota of picophyto
-   qfdopt     =  7.E-6    ! Optimal quota of diatoms
+   qfnopt     =  10.E-6   ! Optimal Fe quota of nanophyto
+   qfpopt     =  10.E-6   ! Optimal Fe quota of picophyto
+   qfdopt     =  10.E-6   ! Optimal quota of diatoms
    qnnmin     =  0.61     ! Minimal N quota for nano
    qnnmax     =  1.25     ! Maximal N quota for nano
    qpnmin     =  0.24     ! Minimal P quota for nano
    qpnmax     =  1.35     ! Maximal P quota for nano
-   qnpmin     =  1.04     ! Minimal N quota for pico
-   qnpmax     =  1.525    ! Maximal N quota for pico
+   qnpmin     =  1.02     ! Minimal N quota for pico
+   qnpmax     =  1.39     ! Maximal N quota for pico
    qppmin     =  0.19     ! Minimal P quota for pico
    qppmax     =  1.15     ! Maximal P quota for pico
-   qndmin     =  0.53     ! Minimal N quota for diatoms
+   qndmin     =  0.51     ! Minimal N quota for diatoms
    qndmax     =  1.25     ! Maximal N quota for diatoms
    qpdmin     =  0.24     ! Minimal P quota for diatoms
    qpdmax     =  1.525    ! Maximal P quota for diatoms
-   qfnmax     =  80E-6    ! Maximal Fe quota for nano
-   qfpmax     =  80E-6    ! Maximal Fe quota for pico
-   qfdmax     =  80E-6    ! Maximal Fe quota for diatoms
+   qfnmax     =  60E-6    ! Maximal Fe quota for nano
+   qfpmax     =  60E-6    ! Maximal Fe quota for pico
+   qfdmax     =  60E-6    ! Maximal Fe quota for diatoms
 /
 !-----------------------------------------------------------------------
@@ -149,6 +150,6 @@
 &namp4zprod    !   parameters for phytoplankton growth for PISCES std  - ln_p4z
 !-----------------------------------------------------------------------
-   pislopen   =  2.       ! P-I slope
-   pisloped   =  2.       ! P-I slope  for diatoms
+   pislopen   =  2.5      ! P-I slope
+   pisloped   =  2.5      ! P-I slope  for diatoms
    xadap      =  0.       ! Adaptation factor to low light
    excretn    =  0.05     ! excretion ratio of phytoplankton
@@ -157,15 +158,15 @@
    chlcnm     =  0.033    ! Maximum Chl/C in nanophytoplankton
    chlcdm     =  0.05     ! Maximum Chl/C in diatoms
-   chlcmin    =  0.004    ! Minimum Chl/c in phytoplankton
-   fecnm      =  40E-6    ! Maximum Fe/C in nanophytoplankton
-   fecdm      =  40E-6    ! Maximum Fe/C in diatoms
-   grosip     =  0.11     ! mean Si/C ratio
+   chlcmin    =  0.003    ! Minimum Chl/c in phytoplankton
+   fecnm      =  60E-6    ! Maximum Fe/C in nanophytoplankton
+   fecdm      =  60E-6    ! Maximum Fe/C in diatoms
+   grosip     =  0.13     ! mean Si/C ratio
 /
 !-----------------------------------------------------------------------
 &namp5zprod    !   parameters for phytoplankton growth for PISCES quota- ln_p5z
 !-----------------------------------------------------------------------
-   pislopen   =  4.       ! P-I slope of nanophytoplankton
-   pislopep   =  4.       ! P-I slope for picophytoplankton
-   pisloped   =  4.       ! P-I slope  for diatoms
+   pislopen   =  5        ! P-I slope of nanophytoplankton
+   pislopep   =  5        ! P-I slope for picophytoplankton
+   pisloped   =  5        ! P-I slope  for diatoms
    excretn    =  0.05     ! excretion ratio of phytoplankton
    excretp    =  0.05     ! excretion ratio of picophytoplankton
@@ -176,6 +177,6 @@
    thetanpm   =  0.3      ! Maximum Chl/N in picophytoplankton
    thetandm   =  0.4      ! Maximum Chl/N in diatoms
-   chlcmin    =  0.004    ! Minimum Chl/c in phytoplankton
-   grosip     =  0.11     ! mean Si/C ratio
+   chlcmin    =  0.003    ! Minimum Chl/c in phytoplankton
+   grosip     =  0.12     ! mean Si/C ratio
 /
 !-----------------------------------------------------------------------
@@ -183,6 +184,5 @@
 !-----------------------------------------------------------------------
    wchln      =  0.01     ! quadratic mortality of phytoplankton
-   wchld      =  0.01     ! maximum quadratic mortality of diatoms
-   wchldm     =  0.03     ! maximum quadratic mortality of diatoms
+   wchld      =  0.03     ! maximum quadratic mortality of diatoms
    mpratn     =  0.01     ! phytoplankton mortality rate
    mpratd     =  0.01     ! Diatoms mortality rate
@@ -193,6 +193,5 @@
    wchln      =  0.01     ! quadratic mortality of nanophytoplankton
    wchlp      =  0.01     ! quadratic mortality of picophytoplankton
-   wchld      =  0.01     ! maximum quadratic mortality of diatoms
-   wchldm     =  0.03     ! maximum quadratic mortality of diatoms
+   wchld      =  0.03     ! maximum quadratic mortality of diatoms
    mpratn     =  0.01     ! nanophytoplankton mortality rate
    mpratp     =  0.01     ! picophytoplankton mortality rate
@@ -205,58 +204,64 @@
    grazrat2    =  0.5      ! maximal mesozoo grazing rate
    resrat2     =  0.005    ! exsudation rate of mesozooplankton
-   mzrat2      =  0.02     ! mesozooplankton mortality rate
+   mzrat2      =  0.01     ! mesozooplankton mortality rate
    xpref2d     =  1.       ! mesozoo preference for diatoms
    xpref2n     =  0.3      ! mesozoo preference for nanophyto.
    xpref2z     =  1.       ! mesozoo preference for microzoo.
    xpref2c     =  0.3      ! mesozoo preference for poc
-   xthresh2zoo = 1E-8      ! zoo feeding threshold for mesozooplankton 
-   xthresh2dia = 1E-8      ! diatoms feeding threshold for mesozooplankton 
-   xthresh2phy = 1E-8      ! nanophyto feeding threshold for mesozooplankton 
-   xthresh2poc = 1E-8      ! poc feeding threshold for mesozooplankton 
+   xthresh2zoo =  1E-8     ! zoo feeding threshold for mesozooplankton 
+   xthresh2dia =  1E-8     ! diatoms feeding threshold for mesozooplankton 
+   xthresh2phy =  1E-8     ! nanophyto feeding threshold for mesozooplankton 
+   xthresh2poc =  1E-8     ! poc feeding threshold for mesozooplankton 
    xthresh2    =  3E-7     ! Food threshold for grazing
    xkgraz2     =  20.E-6   ! half saturation constant for meso grazing
-   epsher2     =  0.35     ! Efficicency of Mesozoo growth
-   epsher2min  =  0.35     ! Minimum efficiency of mesozoo growth
+   epsher2     =  0.4      ! Efficicency of Mesozoo growth
+   epsher2min  =  0.4      ! Minimum efficiency of mesozoo growth
    sigma2      =  0.6      ! Fraction of mesozoo excretion as DOM
    unass2      =  0.3      ! non assimilated fraction of P by mesozoo
-   grazflux    =  3.e3     ! flux-feeding rate
-   ln_dvm_meso = .false.   ! Activates DVM for mesozooplankton
-   xfracmig    = 0.3       ! Fraction of mesozooplankton performing DVM
+   grazflux    =  2.e3     ! flux-feeding rate
+   xsigma2     =  0.5      ! Predation window size
+   xsigma2del  =  1.0      ! Predation window size scaling
+   ln_dvm_meso =  .false.  ! Activates DVM for mesozooplankton
+   xfracmig    =  0.3      ! Fraction of mesozooplankton performing DVM
 /
 !-----------------------------------------------------------------------
 &namp5zmes     !   parameters for mesozooplankton
 !-----------------------------------------------------------------------
-   part2      =  0.75     ! part of calcite not dissolved in mesozoo guts
-   grazrat2   =  0.75     ! maximal mesozoo grazing rate
-   bmetexc2   =  .true.   ! Metabolic use of excess carbon 
-   resrat2    =  0.005    ! exsudation rate of mesozooplankton
-   mzrat2     =  0.03     ! mesozooplankton mortality rate
-   xpref2d    =  1.       ! meso preference for diatoms
-   xpref2n    =  0.3      ! meso preference for nano
-   xpref2z    =  1.       ! meso preference for zoo
-   xpref2m    =  0.1      ! meso preference for zoo
-   xpref2c    =  0.3      ! meso preference for poc
-   xthresh2zoo = 1E-8     ! zoo feeding threshold for mesozooplankton
-   xthresh2dia = 1E-8     ! diatoms feeding threshold for mesozooplankton
-   xthresh2phy = 1E-8     ! nanophyto feeding threshold for mesozooplankton
-   xthresh2mes = 1E-8     ! meso feeding threshold for mesozooplankton
-   xthresh2poc = 1E-8     ! poc feeding threshold for mesozooplankton
-   xthresh2    =  3E-7    ! Food threshold for grazing
-   xkgraz2     =  20.E-6  ! half sturation constant for meso grazing
-   epsher2     =  0.5     ! Efficicency of Mesozoo growth
-   epsher2min  =  0.5     ! Minimum efficiency of mesozoo growth
-   ssigma2     =  0.5     ! Fraction excreted as semi-labile DOM
-   srespir2    =  0.2     ! Active respiration
-   unass2c     =  0.3     ! non assimilated fraction of C by mesozoo
-   unass2n     =  0.3     ! non assimilated fraction of N by mesozoo
-   unass2p     =  0.3     ! non assimilated fraction of P by mesozoo
-   grazflux   =  3.e3     ! flux-feeding rate
+   part2       =  0.75     ! part of calcite not dissolved in mesozoo guts
+   grazrat2    =  0.5      ! maximal mesozoo grazing rate
+   bmetexc2    =  .true.   ! Metabolic use of excess carbon 
+   resrat2     =  0.005    ! exsudation rate of mesozooplankton
+   mzrat2      =  0.01     ! mesozooplankton mortality rate
+   xpref2d     =  1.       ! meso preference for diatoms
+   xpref2n     =  0.3      ! meso preference for nano
+   xpref2z     =  1.       ! meso preference for zoo
+   xpref2m     =  0.       ! meso preference for zoo
+   xpref2c     =  0.3      ! meso preference for poc
+   xthresh2zoo =  1E-8     ! zoo feeding threshold for mesozooplankton
+   xthresh2dia =  1E-8     ! diatoms feeding threshold for mesozooplankton
+   xthresh2phy =  1E-8     ! nanophyto feeding threshold for mesozooplankton
+   xthresh2mes =  1E-8     ! meso feeding threshold for mesozooplankton
+   xthresh2poc =  1E-8     ! poc feeding threshold for mesozooplankton
+   xthresh2    =  3E-7     ! Food threshold for grazing
+   xkgraz2     =  20.E-6   ! half sturation constant for meso grazing
+   epsher2     =  0.5      ! Efficicency of Mesozoo growth
+   epsher2min  =  0.5      ! Minimum efficiency of mesozoo growth
+   ssigma2     =  0.5      ! Fraction excreted as semi-labile DOM
+   srespir2    =  0.2      ! Active respiration
+   unass2c     =  0.3      ! non assimilated fraction of C by mesozoo
+   unass2n     =  0.3      ! non assimilated fraction of N by mesozoo
+   unass2p     =  0.3      ! non assimilated fraction of P by mesozoo
+   xsigma2     =  0.5      ! Predation window size
+   xsigma2del  =  1.0      ! Predation window size scaling
+   grazflux    =  2.e3     ! flux-feeding rate
+   ln_dvm_meso =  .false.  ! Activates DVM for mesozooplankton
+   xfracmig    =  0.25     ! Fraction of mesozooplankton performing DVM
 /
 !-----------------------------------------------------------------------
 &namp4zzoo     !   parameters for microzooplankton for PISCES std      - ln_p4z
 !-----------------------------------------------------------------------
-   part       =  0.5      ! part of calcite not dissolved in microzoo guts
+   part       =  0.75     ! part of calcite not dissolved in microzoo guts
    grazrat    =  2.0      ! maximal zoo grazing rate
-   resrat     =  0.03     ! Linear mortality rate of zooplankton
+   resrat     =  0.02     ! Linear mortality rate of zooplankton
    mzrat      =  0.005    ! zooplankton mortality rate
    xprefc     =  0.1      ! Microzoo preference for POM
@@ -268,22 +273,24 @@
    xthresh    =  3.E-7    ! Food threshold for feeding
    xkgraz     =  20.E-6   ! half sturation constant for grazing
-   epsher     =  0.35     ! Efficiency of microzoo growth
-   epshermin  =  0.35     ! Minimum efficiency of microzoo growth
+   epsher     =  0.4     ! Efficiency of microzoo growth
+   epshermin  =  0.4     ! Minimum efficiency of microzoo growth
    sigma1     =  0.6      ! Fraction of microzoo excretion as DOM
    unass      =  0.3      ! non assimilated fraction of phyto by zoo
+   xsigma     =  0.5      ! Predation window size
+   xsigmadel  =  1.0      ! Predation window size scaling
 /
 !-----------------------------------------------------------------------
 &namp5zzoo     !   parameters for microzooplankton
 !-----------------------------------------------------------------------
-   part       =  0.5      ! part of calcite not dissolved in microzoo gutsa
-   grazrat    =  3.0      ! maximal zoo grazing rate
+   part       =  0.75     ! part of calcite not dissolved in microzoo gutsa
+   grazrat    =  2.0      ! maximal zoo grazing rate
    bmetexc    =  .true.   ! Metabolic use of excess carbon
-   resrat     =  0.03     ! exsudation rate of zooplankton
-   mzrat      =  0.004    ! zooplankton mortality rate
+   resrat     =  0.02     ! exsudation rate of zooplankton
+   mzrat      =  0.005    ! zooplankton mortality rate
    xprefc     =  0.1      ! Microzoo preference for POM
    xprefn     =  1.0      ! Microzoo preference for Nanophyto
    xprefp     =  1.0      ! Microzoo preference for picophyto
    xprefd     =  0.8      ! Microzoo preference for Diatoms
-   xprefz     =  0.1      ! Microzoo preference for microzooplankton
+   xprefz     =  0.       ! Microzoo preference for microzooplankton
    xthreshdia =  1.E-8    ! Diatoms feeding threshold for microzooplankton
    xthreshphy =  1.E-8    ! Nanophyto feeding threshold for microzooplankton
@@ -300,4 +307,6 @@
    unassn     =  0.3      ! non assimilated fraction of N by zoo
    unassp     =  0.3      ! non assimilated fraction of P by zoo
+   xsigma     =  0.5      ! Predation window size
+   xsigmadel  =  1.0      ! Predation window size scaling
 /
 !-----------------------------------------------------------------------
@@ -305,30 +314,28 @@
 !-----------------------------------------------------------------------
    ln_ligvar =  .false.   ! variable ligand concentration
-   xlam1     =  0.005      ! scavenging rate of Iron by biogenic particles
+   xlam1     =  0.05      ! scavenging rate of Iron by biogenic particles
    xlamdust  =  150.0     ! Scavenging rate of Iron by dust
-   ligand    =  7E-10     ! Ligands concentration 
+   ligand    =  1E-9     ! Ligands concentration 
    kfep      =  0.01      ! Nanoparticle formation rate constant
+   scaveff   =  1.0       ! Fraction of scavenged Fe that goes to POFe
 /
 !-----------------------------------------------------------------------  
 &nampisrem     !   parameters for remineralization
 !-----------------------------------------------------------------------
-   xremik    =  0.3       ! remineralization rate of DOC
    nitrif    =  0.05      ! NH4 nitrification rate
    xsirem    =  0.003     ! remineralization rate of Si
    xsiremlab =  0.03      ! fast remineralization rate of Si
    xsilab    =  0.5       ! Fraction of labile biogenic silica
-   feratb    =  10.E-6    ! Fe/C quota in bacteria
+   feratb    =  30.E-6    ! Fe/C quota in bacteria
    xkferb    =  3E-10     ! Half-saturation constant for bacteria Fe/C
 !                         ! ln_p5z
-   xremikc   =  0.25      ! remineralization rate of DOC
-   xremikn   =  0.3       ! remineralization rate of DON
-   xremikp   =  0.4       ! remineralization rate of DOP
-!   feratb    =  20E-6     ! Bacterial Fe/C ratio
-!   xkferb    =  3E-10     ! Half-saturation constant for bact. Fe/C
+   xremikc   =  0.4       ! remineralization rate of DOC
+   xremikn   =  0.4       ! remineralization rate of DON
+   xremikp   =  0.5       ! remineralization rate of DOP
 /
 !-----------------------------------------------------------------------
 &nampispoc     !   parameters for organic particles
 !-----------------------------------------------------------------------
-   xremip    =  0.03      ! remineralisation rate of POC
+   xremip    =  0.035     ! remineralisation rate of POC
    jcpoc     =  15        ! Number of lability classes
    rshape    =  1.0       ! Shape of the gamma function
@@ -341,6 +348,6 @@
 &nampiscal     !   parameters for Calcite chemistry
 !-----------------------------------------------------------------------
-   kdca       =  6.       ! calcite dissolution rate constant (1/time)
-   nca        =  1.       ! order of dissolution reaction (dimensionless)
+   kdca       =  3.       ! calcite dissolution rate constant (1/time)
+   nca        =  2.       ! order of dissolution reaction (dimensionless)
 /
 !-----------------------------------------------------------------------
@@ -377,5 +384,5 @@
    icefeinput  =  15.e-9   ! Iron concentration in sea ice
    nitrfix     =  1.e-7    ! Nitrogen fixation rate
-   diazolight  =  50.      ! Diazotrophs sensitivity to light (W/m2)
+   diazolight  =  30.      ! Diazotrophs sensitivity to light (W/m2)
    concfediaz  =  1.e-10   ! Diazotrophs half-saturation Cste for Iron
    hratio      =  1.e+7    ! Fe to 3He ratio assumed for vent iron supply 
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/namelist_sediment_ref
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/namelist_sediment_ref	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/namelist_sediment_ref	(revision 14276)
@@ -18,12 +18,12 @@
 &nam_geom     !  Characteristics of the sediment columns
 !-----------------------------------------------------------------------
-  jpksed    = 11       ! Number of vertical layers
-  sedzmin = 0.3        ! Minimum vertical spacing
-  sedhmax = 10.0       ! Maximum depth of the sediment
-  sedkth  = 7.0        ! Default parameters
+  jpksed  = 18         ! Number of vertical layers
+  sedzmin = 0.1        ! Minimum vertical spacing
+  sedhmax = 15.0       ! Maximum depth of the sediment
+  sedkth  = 9.0        ! Default parameters
   sedacr  = 2.0        ! Default parameters
-  porsurf   = 0.9      ! Sediment porosity at the surface
-  porinf    = 0.75     ! Sediment porosity at infinite depth
-  rhox      = 0.33     ! Length scale of porosity variation (cm)
+  porsurf = 0.9        ! Sediment porosity at the surface
+  porinf  = 0.75       ! Sediment porosity at infinite depth
+  rhox    = 0.33       ! Length scale of porosity variation (cm)
 /
 !-----------------------------------------------------------------------
@@ -57,5 +57,6 @@
 !
    seddiag3d(1)   = 'SedpH      ' , 'pH                                     ',  '-      '
-   seddiag3d(2)   = 'SedCO3por  ' , 'Dissolved CO3 cocnentration            ',  'mol/L'
+   seddiag3d(2)   = 'SedCO3por  ' , 'Dissolved CO3 concentration            ',  'mol/L'
+   seddiag3d(3)   = 'Sedligand'   , 'ligand concentration                   ',  'mol/L'
    seddiag2d(1)   = 'FlxSi      ' , 'Silicate flux                          ',  'mol/cm2/s'
    seddiag2d(2)   = 'FlxO2      ' , 'Dissolved Oxygen Flux                  ',  'mol/L'
@@ -67,7 +68,15 @@
    seddiag2d(8)   = 'FlxH2S     ' , 'Total H2S flux                         ',  'mol/cm2/s'
    seddiag2d(9)   = 'FlxSO4     ' , 'Total SO4 flux                         ',  'mol/cm2/s'
-   seddiag2d(10)   = 'FlxFe2     ' , 'Total Fe2+ flux                        ',  'mol/cm2/s'
-   seddiag2d(11)  = 'Flxtot     ' , 'Total net sediment burial rate         ',  'cm/s'
-   seddiag2d(12)  = 'dzdep      ' , 'Sedimentation rate                     ',  'cm/s'
+   seddiag2d(10)  = 'FlxFe2     ' , 'Total Fe2+ flux                        ',  'mol/cm2/s'
+   seddiag2d(11)  = 'FlxBSi     ' , 'Total net BSi burial rate              ',  'g/cm2/s'
+   seddiag2d(12)  = 'FlxClay    ' , 'Total net CLay burial rate             ',  'g/cm2/s'
+   seddiag2d(13)  = 'FlxPOC     ' , 'Total net POC burial rate              ',  'g/cm2/s'
+   seddiag2d(14)  = 'FlxCaCO3   ' , 'Total net CaCO3 burial rate            ',  'g/cm2/s'
+   seddiag2d(15)  = 'FlxPOS     ' , 'Total net POS burial rate              ',  'g/cm2/s'
+   seddiag2d(16)  = 'FlxPOR     ' , 'Total net POR burial rate              ',  'g/cm2/s'
+   seddiag2d(17)  = 'FlxFeO     ' , 'Total net FeO burial rate              ',  'g/cm2/s'
+   seddiag2d(18)  = 'FlxFeS     ' , 'Total net FeS burial rate              ',  'g/cm2/s'
+   seddiag2d(19)  = 'Flxtot     ' , 'Total burial flux                      ',  'mol/cm2/s'
+   seddiag2d(20)  = 'dzdep      ' , 'Sedimentation rate                     ',  'cm/s'
 /
 !-----------------------------------------------------------------------
@@ -76,5 +85,5 @@
    rcopal   =  0.1      ! Reactivity for clay  [l.mol-1.an-1]
    dcoef    =  8.0E-6   ! Diffusion coefficient in sediments (in cm**2/s)
-   rccal    =  0.014    ! Dissolution rate of calcite [in yr-1]
+   rccal    =  1.4      ! Dissolution rate of calcite [in yr-1]
    rcligc   =  0.1      ! Degradation rate of ligands [in yr-1]
    ratligc  =  1.E-4    ! L/C ratio in POC
@@ -92,8 +101,8 @@
    rcorgr   =  1.E-4   ! Reactivity for refractory POC [an-1]
    rcnh4    =  1E7     ! Reactivity for O2/NH4 [l.mol-1.an-1]
-   rch2s    =  1E5     ! Reactivity for O2/H2S [l.mol-1.an-1]
+   rch2s    =  1.6E5   ! Reactivity for O2/H2S [l.mol-1.an-1]
    rcfe2    =  1E8     ! Reactivity for O2/Fe2+ [l.mol-1.an-1]
-   rcfeh2s  =  1E4     ! Reactivity for FEOH/H2S [l.mol-1.an-1]
-   rcfes    =  1E5     ! Reactivity for FE2+/H2S [l.mol-1.an-1]
+   rcfeh2s  =  8E3     ! Reactivity for FEOH/H2S [l.mol-1.an-1]
+   rcfes    =  1E6     ! Reactivity for FE2+/H2S [l.mol-1.an-1]
    rcfeso   =  3E5     ! Reactivity for FES/O2 [l.mol-1.an-1]
    xksedo2  =  1.E-6   ! Half-saturation constant for oxic remin [mol/l]
@@ -105,9 +114,9 @@
 &nam_btb     !   parameters for bioturbation and bioirrigation
 !-----------------------------------------------------------------------
-   dbiot       =  5.        ! Bioturbation coefficient (in cm2/yr)
-   ln_btbz     =  .false.   ! Depth varying bioturbation 
-   dbtbzsc     =  10.0      ! coefficient for btb attenuation. If no variation, depth of the mixed layer
+   dbiot       =  10.       ! Bioturbation coefficient (in cm2/yr)
+   ln_btbz     =  .true.    ! Depth varying bioturbation 
+   dbtbzsc     =  6.0       ! coefficient for btb attenuation. If no variation, depth of the mixed layer
    adsnh4      =  1.3       ! Adsorption coefficient of NH4
-   ln_irrig    =  .true.   ! Bioirrigation in sediment
+   ln_irrig    =  .true.    ! Bioirrigation in sediment
    xirrzsc     =  1.5       ! Vertical scale of the attentuation of bioirrigation [cm]
 / 
@@ -115,5 +124,5 @@
 &nam_rst     !   parameters for restarting the sediment module
 !-----------------------------------------------------------------------
-   ln_rst_sed   =  .false.       !  restart flag
+   ln_rst_sed   =  .false.         !  restart flag
    cn_sedrst_in  = "restart_sed"   !  name of the input restart file name of the sediment module
    cn_sedrst_indir = "."           !  directory from which to read input sediment restarts
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zagg.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zagg.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zagg.F90	(revision 14276)
@@ -75,12 +75,12 @@
                   ! Part I : Coagulation dependent on turbulence
                   ! The stickiness has been assumed to be 0.1
-                  zagg1 = 25.9  * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
-                  zagg2 = 4452. * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
+                  zagg1 = 12.5  * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
+                  zagg2 = 169.7 * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
 
                   ! Part II : Differential settling
                   ! Aggregation of small into large particles
                   ! The stickiness has been assumed to be 0.1
-                  zagg3 =  47.1 * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
-                  zagg4 =  3.3  * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
+                  zagg3 =  8.63  * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
+                  zagg4 =  132.8 * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
 
                   zagg   = zagg1 + zagg2 + zagg3 + zagg4
@@ -92,14 +92,13 @@
                   ! 3rd term is differential settling of DOC-POC
                   ! 1/3 of DOC is supposed to experience aggregation (HMW)
-                  zaggdoc  = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
-                  &            + 2.4 * xstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc)
+                  zaggdoc  = ( ( 12.0 * 0.3 * trb(ji,jj,jk,jpdoc) + 9.05 * trb(ji,jj,jk,jppoc) ) * zfact       &
+                  &            + 2.49 * xstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc)
                   ! transfer of DOC to GOC : 
                   ! 1st term is shear aggregation
-                  ! 2nd term is differential settling 
-                  ! 1/3 of DOC is supposed to experience aggregation (HMW)
-                  zaggdoc2 = ( 3.53E3 * zfact + 0.1 * xstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
+                  ! 1/3 of DOC is supposed to experience aggregation (HMW)
+                  zaggdoc2 = ( 1.94 * zfact + 1.37 * xstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
                   ! tranfer of DOC to POC due to brownian motion
                   ! The temperature dependency has been omitted.
-                  zaggdoc3 =  114. * 0.3 * trb(ji,jj,jk,jpdoc) *xstep * 0.3 * trb(ji,jj,jk,jpdoc)
+                  zaggdoc3 =  ( 127.8 * 0.3 * trb(ji,jj,jk,jpdoc) + 725.7 * trb(ji,jj,jk,jppoc) ) * xstep * 0.3 * trb(ji,jj,jk,jpdoc)
 
                   !  Update the trends
@@ -150,6 +149,6 @@
                   ! 3rd term is differential settling of DOC-POC
                   ! 1/3 of DOC is supposed to experience aggregation (HMW)
-                  zaggtmp = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
-                  &            + 2.4 * xstep * trb(ji,jj,jk,jppoc) )
+                  zaggtmp = ( ( 0.37 * 0.3 * trb(ji,jj,jk,jpdoc) + 20.5 * trb(ji,jj,jk,jppoc) ) * zfact       &
+                  &            + 0.15 * xstep * trb(ji,jj,jk,jppoc) )
                   zaggdoc  = zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
                   zaggdon  = zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
@@ -160,5 +159,5 @@
                   ! 2nd term is differential settling 
                   ! 1/3 of DOC is supposed to experience aggregation (HMW)
-                  zaggtmp = ( 3.53E3 * zfact + 0.1 * xstep ) * trb(ji,jj,jk,jpgoc)
+                  zaggtmp = 655.4 * zfact * trb(ji,jj,jk,jpgoc)
                   zaggdoc2 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
                   zaggdon2 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
@@ -167,5 +166,5 @@
                   ! tranfer of DOC to POC due to brownian motion
                   ! 1/3 of DOC is supposed to experience aggregation (HMW)
-                  zaggtmp = ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) * xstep
+                  zaggtmp = ( 260.2 * 0.3 * trb(ji,jj,jk,jpdoc) +  418.5 * trb(ji,jj,jk,jppoc) ) * xstep
                   zaggdoc3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
                   zaggdon3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 14276)
@@ -110,8 +110,8 @@
 
       ! Update of the size of the different phytoplankton groups
-      sized(:,:,:) = sizeda(:,:,:)
-      sizen(:,:,:) = sizena(:,:,:)
+      sized(:,:,:) = MAX(1.0, sizeda(:,:,:) )
+      sizen(:,:,:) = MAX(1.0, sizena(:,:,:) )
       IF (ln_p5z) THEN
-         sizep(:,:,:) = sizepa(:,:,:)
+         sizep(:,:,:) = MAX(1.0, sizepa(:,:,:) )
       ENDIF
       !                                                             !
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zfechem.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 14276)
@@ -30,4 +30,5 @@
    REAL(wp), PUBLIC ::   ligand       !: ligand concentration in the ocean 
    REAL(wp), PUBLIC ::   kfep         !: rate constant for nanoparticle formation
+   REAL(wp), PUBLIC ::   scaveff      !: Fraction of scavenged iron that is considered as being subject to solubilization
 
    !!----------------------------------------------------------------------
@@ -50,12 +51,12 @@
       !
       INTEGER  ::   ji, jj, jk, jic, jn
-      REAL(wp) ::   zdep, zlam1a, zlam1b, zlamfac
-      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll, fe3sol, zligco
-      REAL(wp) ::   zdenom1, zscave, zaggdfea, zaggdfeb, zcoag, ztrc, zdust
-      REAL(wp) ::   zdenom2, ztfe, zhplus, zxlam, zaggliga, zaggligb
+      REAL(wp) ::   zlam1a, zlam1b
+      REAL(wp) ::   zkeq, zfesatur, zfecoll, fe3sol, zligco
+      REAL(wp) ::   zscave, zaggdfea, zaggdfeb, ztrc, zdust, zklight
+      REAL(wp) ::   ztfe, zhplus, zxlam, zaggliga, zaggligb
       REAL(wp) ::   zrfact2
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zTL1, zFe3, ztotlig, precip, zFeL1
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zcoll3d, zscav3d, zlcoll3d
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zTL1, zFe3, ztotlig, precip, precipno3, zFeL1
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zcoll3d, zscav3d, zlcoll3d, zprecip3d
       !!---------------------------------------------------------------------
       !
@@ -66,12 +67,12 @@
 
       ! Total ligand concentration : Ligands can be chosen to be constant or variable
-      ! Parameterization from Tagliabue and Voelker (2011)
+      ! Parameterization from Pham and Ito (2018)
       ! -------------------------------------------------
       IF( ln_ligvar ) THEN
-         ztotlig(:,:,:) =  0.09 * trb(:,:,:,jpdoc) * 1E6 + ligand * 1E9
+         ztotlig(:,:,:) =  0.09 * 0.667 * trb(:,:,:,jpdoc) * 1E6 + ligand * 1E9 + MAX(0., chemo2(:,:,:) - trb(:,:,:,jpoxy) ) / 400.E-6
          ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
       ELSE
         IF( ln_ligand ) THEN  ;   ztotlig(:,:,:) = trb(:,:,:,jplgw) * 1E9
-        ELSE                  ;   ztotlig(:,:,:) = ligand * 1E9
+        ELSE                  ;   ztotlig(:,:,:) = ligand * 1E9 
         ENDIF
       ENDIF
@@ -79,5 +80,5 @@
       ! ------------------------------------------------------------
       !  from Aumont and Bopp (2006)
-      ! This model is based on one ligand and Fe' 
+      ! This model is based on one ligand, Fe2+ and Fe3+ 
       ! Chemistry is supposed to be fast enough to be at equilibrium
       ! ------------------------------------------------------------
@@ -87,17 +88,18 @@
                zTL1(ji,jj,jk)  = ztotlig(ji,jj,jk)
                zkeq            = fekeq(ji,jj,jk)
+               zklight         = 4.77E-7 * etot(ji,jj,jk) * 0.5 / 10**-6.3
                zfesatur        = zTL1(ji,jj,jk) * 1E-9
-               ztfe            = trb(ji,jj,jk,jpfer) 
+               ztfe            = (1.0 + zklight) * trb(ji,jj,jk,jpfer) 
                ! Fe' is the root of a 2nd order polynom
-               zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe )               &
-                  &              + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2       &
+               zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq + zklight + consfe3(ji,jj,jk)/10**-6.3 - zkeq * trb(ji,jj,jk,jpfer) )               &
+                  &              + SQRT( ( 1. + zfesatur * zkeq + zklight + consfe3(ji,jj,jk)/10**-6.3 - zkeq * trb(ji,jj,jk,jpfer) )**2       &
                   &              + 4. * ztfe * zkeq) ) / ( 2. * zkeq )
-               zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9
-               zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) )
+               zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) - zFe3(ji,jj,jk) )
            END DO
          END DO
       END DO
-         !
-
+      !
+      plig(:,:,:) =  MAX( 0., ( zFeL1(:,:,:) / ( trb(:,:,:,jpfer) + rtrn ) ) )
+      !
       zdust = 0.         ! if no dust available
       DO jk = 1, jpkm1
@@ -113,32 +115,17 @@
                &         + fesol(ji,jj,jk,5) / zhplus )
                !
-               zfeequi = zFe3(ji,jj,jk) * 1E-9
-               zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
+               zfecoll = 0.5 * zFeL1(ji,jj,jk)
                ! precipitation of Fe3+, creation of nanoparticles
-               precip(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * 1E-9 - fe3sol ) ) * kfep * xstep
+               precip(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) - fe3sol ) ) * kfep * xstep * ( 1.0 - nitrfac(ji,jj,jk) ) 
+               ! Precipitation of Fe2+ due to oxidation by NO3 (Croot et al., 2019)
+               ! This occurs in anoxic waters only
+               precipno3(ji,jj,jk) = 2.0 * 130.0 * trb(ji,jj,jk,jpno3) * nitrfac(ji,jj,jk) * xstep * zFe3(ji,jj,jk)
                !
                ztrc   = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 
+               ztrc = MAX( rtrn, ztrc )
                IF( ln_dust )  zdust  = dust(ji,jj) / ( wdust / rday ) * tmask(ji,jj,jk)
-               IF (ln_ligand) THEN
-                  zxlam  = xlam1 * MAX( 1.E-3, EXP(-2 * etot(ji,jj,jk) / 10. ) * (1. - EXP(-2 * trb(ji,jj,jk,jpoxy) / 100.E-6 ) ))
-               ELSE
-                  zxlam  = xlam1 * 1.0
-               ENDIF
-               zlam1b = 3.e-5 + xlamdust * zdust + zxlam * ztrc
-               zscave = zfeequi * zlam1b * xstep
-
-               ! Compute the different ratios for scavenging of iron
-               ! to later allocate scavenged iron to the different organic pools
-               ! ---------------------------------------------------------
-               zdenom1 = zxlam * trb(ji,jj,jk,jppoc) / zlam1b
-               zdenom2 = zxlam * trb(ji,jj,jk,jpgoc) / zlam1b
-
-               ! Increased scavenging for very high iron concentrations found near the coasts 
-               ! due to increased lithogenic particles and let say it is unknown processes (precipitation, ...)
-               !  -----------------------------------------------------------
-               zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. )
-               zlamfac = MIN( 1.  , zlamfac )
-               zdep    = MIN( 1., 1000. / gdept_n(ji,jj,jk) )
-               zcoag   = 1E-4 * ( 1. - zlamfac ) * zdep * xstep * trb(ji,jj,jk,jpfer)
+               zxlam  = MAX( 1.E-3, (1. - EXP(-2 * trb(ji,jj,jk,jpoxy) / 100.E-6 ) ))
+               zlam1b = 3.e-5 + ( xlamdust * zdust + xlam1 * ztrc ) * zxlam
+               zscave = zFe3(ji,jj,jk) * zlam1b * xstep
 
                ! Compute the coagulation of colloidal iron. This parameterization 
@@ -146,17 +133,32 @@
                ! It requires certainly some more work as it is very poorly constrained.
                ! ----------------------------------------------------------------
-               zlam1a   = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
-                   &    + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) )
+               zlam1a   = ( 12.0  * 0.3 * trb(ji,jj,jk,jpdoc) + 9.05  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
+                   &    + ( 2.49  * trb(ji,jj,jk,jppoc) )     &
+                   &    + ( 127.8 * 0.3 * trb(ji,jj,jk,jpdoc) + 725.7 * trb(ji,jj,jk,jppoc) )
                zaggdfea = zlam1a * xstep * zfecoll
                !
-               zlam1b   = 3.53E3 * trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk)
+               zlam1b   = ( 1.94 * xdiss(ji,jj,jk) + 1.37 ) * trb(ji,jj,jk,jpgoc)
                zaggdfeb = zlam1b * xstep * zfecoll
                !
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb &
-               &                     - zcoag - precip(ji,jj,jk)
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea
-               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb
-               zscav3d(ji,jj,jk)   = zscave
+               &                     - precip(ji,jj,jk) - precipno3(ji,jj,jk)
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * scaveff * trb(ji,jj,jk,jppoc) / ztrc
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * scaveff * trb(ji,jj,jk,jpgoc) / ztrc
+
+               ! Precipitated iron is supposed to be permanently lost.
+               ! Scavenged iron is supposed to be released back to seawater
+               ! when POM is solubilized. This is highly uncertain as probably
+               ! a significant part of it may be rescavenged back onto 
+               ! the particles. An efficiency factor is applied that is read
+               ! in the namelist. 
+               ! See for instance Tagliabue et al. (2019).
+               ! Aggregated FeL is considered as biogenic Fe as it 
+               ! probably remains  complexed when the particle is solubilized.
+               ! -------------------------------------------------------------
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zaggdfea
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggdfeb
+               zscav3d(ji,jj,jk)   = zscave 
                zcoll3d(ji,jj,jk)   = zaggdfea + zaggdfeb
+               zprecip3d(ji,jj,jk) = precip(ji,jj,jk) + precipno3(ji,jj,jk)
                !
             END DO
@@ -175,14 +177,15 @@
 
                   ! Coagulation of ligands due to various processes (Brownian, shear, diff. sedimentation
-                  ! Coefficients are taken from the p4zagg
+                  ! Coefficients are taken from p4zagg
                   ! -------------------------------------------------------------------------------------
-                  zlam1a   = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
-                      &    + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) )
+                  zlam1a   = ( 12.0  * 0.3 * trb(ji,jj,jk,jpdoc) + 9.05  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
+                      &    + ( 2.49  * trb(ji,jj,jk,jppoc) )     &
+                      &    + ( 127.8 * 0.3 * trb(ji,jj,jk,jpdoc) + 725.7 * trb(ji,jj,jk,jppoc) )
                   !
-                  zlam1b   = 3.53E3 *   trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk)
-
+                  zlam1b   = ( 1.94 * xdiss(ji,jj,jk) + 1.37 ) * trb(ji,jj,jk,jpgoc)
                   ! 50% of the ligands are supposed to be in the colloidal size fraction
-                  zligco   = 0.5 * trn(ji,jj,jk,jplgw)
-                  zaggliga = zlam1a * xstep * zligco
+                  ! as for FeL
+                  zligco   = 0.5 * trb(ji,jj,jk,jplgw)
+                  zaggliga = zlam1a * xstep * zligco 
                   zaggligb = zlam1b * xstep * zligco
                   tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) - zaggliga - zaggligb
@@ -191,6 +194,4 @@
             END DO
          END DO
-         !
-         plig(:,:,:) =  MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( trb(:,:,:,jpfer) +rtrn ) ) )
          !
       ENDIF
@@ -206,4 +207,5 @@
          IF( iom_use("FESCAV") )  CALL iom_put("FESCAV" , zscav3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
          IF( iom_use("FECOLL") )  CALL iom_put("FECOLL" , zcoll3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
+         IF( iom_use("FEPREC") )  CALL iom_put("FEPREC" , zprecip3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
          IF( iom_use("LGWCOLL"))  CALL iom_put("LGWCOLL", zlcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 )
       ENDIF
@@ -234,5 +236,5 @@
       INTEGER ::   ios   ! Local integer 
       !!
-      NAMELIST/nampisfer/ ln_ligvar, xlam1, xlamdust, ligand, kfep 
+      NAMELIST/nampisfer/ ln_ligvar, xlam1, xlamdust, ligand, kfep, scaveff 
       !!----------------------------------------------------------------------
       !
@@ -258,4 +260,5 @@
          WRITE(numout,*) '      ligand concentration in the ocean         ligand       =', ligand
          WRITE(numout,*) '      rate constant for nanoparticle formation  kfep         =', kfep
+         WRITE(numout,*) '      Scavenged iron that is added to POFe      scaveff      =', scaveff
       ENDIF
       ! 
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zflx.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zflx.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zflx.F90	(revision 14276)
@@ -359,5 +359,5 @@
       IF( ln_presatm ) THEN
          CALL fld_read( kt, 1, sf_patm )               !* input Patm provided at kt + 1/2
-         patm(:,:) = sf_patm(1)%fnow(:,:,1)                        ! atmospheric pressure
+         patm(:,:) = sf_patm(1)%fnow(:,:,1)/101325.0     ! atmospheric pressure
       ENDIF
       !
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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zint.F90	(revision 14276)
@@ -36,5 +36,5 @@
       !
       INTEGER  :: ji, jj                 ! dummy loop indices
-      REAL(wp) :: zvar                   ! local variable
+      REAL(wp) :: zrum, zcodel, zargu, zvar
       !!---------------------------------------------------------------------
       !
@@ -44,9 +44,7 @@
       ! -------------------------------------------
       ! Generic temperature dependence (Eppley, 1972)
-      tgfunc (:,:,:) = EXP( 0.063913 * tsn(:,:,:,jp_tem) )
+      tgfunc (:,:,:) = EXP( 0.0631 * tsn(:,:,:,jp_tem) )
       ! Temperature dependence of mesozooplankton (Buitenhuis et al. (2005))
-      tgfunc2(:,:,:) = EXP( 0.07608  * tsn(:,:,:,jp_tem) )
-      ! Temperature dependence of picophytoplankton (Stawiarsky et al., 2016)
-      tgfunc3(:,:,:) = EXP( 0.0825   * tsn(:,:,:,jp_tem) )
+      tgfunc2(:,:,:) = EXP( 0.0761  * tsn(:,:,:,jp_tem) )
 
       ! Computation of the silicon dependant half saturation  constant for silica uptake
@@ -69,4 +67,19 @@
       ENDIF
       !
+      ! compute the day length depending on latitude and the day
+      ! Astronomical parameterization taken from HAMOCC3
+      zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
+      zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
+
+      ! day length in hours
+      strn(:,:) = 0.
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
+            zargu = MAX( -1., MIN(  1., zargu ) )
+            strn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
+         END DO
+      END DO
+      !
       IF( ln_timing )   CALL timing_stop('p4z_int')
       !
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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlim.F90	(revision 14276)
@@ -63,4 +63,5 @@
    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
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xqfuncfecd, xqfuncfecn
 
    ! Coefficient for iron limitation following Flynn and Hipkin (1999)
@@ -90,27 +91,19 @@
       !
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zno3, zferlim, zcoef
+      REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zcoef
       REAL(wp) ::   z1_trbdia, z1_trbphy, ztem1, ztem2, zetot1, zetot2
-      REAL(wp) ::   zdenom, zratio, zironmin
+      REAL(wp) ::   zdenom, zratio, zironmin, zbactno3, zbactnh4
       REAL(wp) ::   zconc1d, zconc1dnh4, zconc0n, zconc0nnh4   
       REAL(wp) ::   fananof, fadiatf, znutlim, zfalim
+      REAL(wp) ::   znutlimtot, zlimno3, zlimnh4, zbiron
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_lim')
       !
-      sizena(:,:,:) = 0.0  ;  sizeda(:,:,:) = 0.0
+      sizena(:,:,:) = 1.0  ;  sizeda(:,:,:) = 1.0
       !
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
-               
-               ! Tuning of the iron concentration to a minimum level that 
-               ! is set to the detection limit
-               ! --------------------------------------------------------
-               zno3    = trb(ji,jj,jk,jpno3) / 40.e-6
-               zferlim = MAX( 3e-11 * zno3 * zno3, 5e-12 )
-               zferlim = MIN( zferlim, 7e-11 )
-               trb(ji,jj,jk,jpfer) = MAX( trb(ji,jj,jk,jpfer), zferlim )
-
                ! Computation of a variable Ks of diatoms taking into account
                ! that increasing biomass is made of generally bigger cells
@@ -134,9 +127,10 @@
 
                ! Nanophytoplankton
-               znutlim = biron(ji,jj,jk) / concnfe(ji,jj,jk)
+               zbiron = ( 75.0 * ( 1.0 - plig(ji,jj,jk) ) + plig(ji,jj,jk) ) * biron(ji,jj,jk)
+               znutlim = zbiron / 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)
+               znutlim = zbiron / concdfe(ji,jj,jk)
                fadiatf = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
 
@@ -144,9 +138,11 @@
                ! heterotrophic bacteria
                ! -------------------------------------------------
-               zdenom = 1. /  ( concbno3 * concbnh4 + concbnh4 * trb(ji,jj,jk,jpno3) + concbno3 * trb(ji,jj,jk,jpnh4) )
-               xnanono3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * concbnh4 * zdenom
-               xnanonh4(ji,jj,jk) = trb(ji,jj,jk,jpnh4) * concbno3 * zdenom
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( concbno3 + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( concbno3 + trb(ji,jj,jk,jpno3) )
+               znutlimtot = ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) ) / ( concbno3 + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               zbactnh4 = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               zbactno3 = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
-               zlim1    = xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk)
+               zlim1    = zbactno3 + zbactnh4
                zlim2    = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concbnh4 )
                zlim3    = trb(ji,jj,jk,jpfer) / ( concbfe + trb(ji,jj,jk,jpfer) )
@@ -166,10 +162,12 @@
                ! Limitation of Fe uptake (Quota formalism)
                zfalim = (1.-fananof) / fananof
-               xnanofer(ji,jj,jk) = (1. - fananof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * concnfe(ji,jj,jk) )
+               xnanofer(ji,jj,jk) = (1. - fananof) * zbiron / ( zbiron + 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
-               xnanonh4(ji,jj,jk) = trb(ji,jj,jk,jpnh4) * zconc0n    * zdenom
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( zconc0n + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( zconc0n + trb(ji,jj,jk,jpno3) )
+               znutlimtot = ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) ) / ( zconc0n + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               xnanonh4(ji,jj,jk) = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               xnanono3(ji,jj,jk) = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                zlim1    = xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk)
@@ -181,4 +179,5 @@
                ! proposed by Flynn and Hipkin (1999)
                zironmin = xcoef1 * trb(ji,jj,jk,jpnch) * z1_trbphy + xcoef2 * zlim1 + xcoef3 * xnanono3(ji,jj,jk)
+               xqfuncfecn(ji,jj,jk) = zironmin + qnfelim
                zlim3    = MAX( 0.,( zratio - zironmin ) / qnfelim )
                xnanopo4(ji,jj,jk) = zlim2
@@ -190,10 +189,12 @@
                ! Limitation of Fe uptake (Quota formalism)
                zfalim = (1.-fadiatf) / fadiatf
-               xdiatfer(ji,jj,jk) = (1. - fadiatf) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * concdfe(ji,jj,jk) )
+               xdiatfer(ji,jj,jk) = (1. - fadiatf) * zbiron / ( zbiron + 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
-               xdiatnh4(ji,jj,jk) = trb(ji,jj,jk,jpnh4) * zconc1d    * zdenom
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( zconc1d + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( zconc1d + trb(ji,jj,jk,jpno3) )
+               znutlimtot = ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) ) / ( zconc1d + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               xdiatnh4(ji,jj,jk) = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn ) 
+               xdiatno3(ji,jj,jk) = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                zlim1    = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk)
@@ -206,4 +207,5 @@
                ! proposed by Flynn and Hipkin (1999)
                zironmin = xcoef1 * trb(ji,jj,jk,jpdch) * z1_trbdia + xcoef2 * zlim1 + xcoef3 * xdiatno3(ji,jj,jk)
+               xqfuncfecd(ji,jj,jk) = zironmin + qdfelim
                zlim4    = MAX( 0., ( zratio - zironmin ) / qdfelim )
                xdiatpo4(ji,jj,jk) = zlim2
@@ -238,12 +240,11 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zlim1 =  ( trb(ji,jj,jk,jpno3) * concnnh4 + trb(ji,jj,jk,jpnh4) * concnno3 )    &
-                  &   / ( concnno3 * concnnh4 + concnnh4 * trb(ji,jj,jk,jpno3) + concnno3 * trb(ji,jj,jk,jpnh4) ) 
+               zlim1  = xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) 
                zlim2  = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concnnh4 )
-               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  5.E-11   )
-               ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) )
+               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  6.E-11   )
+               ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) + 1.8)
                ztem2  = tsn(ji,jj,jk,jp_tem) - 10.
                zetot1 = MAX( 0., etot_ndcy(ji,jj,jk) - 1.) / ( 4. + etot_ndcy(ji,jj,jk) ) 
-               zetot2 = 30. / ( 30. + etot_ndcy(ji,jj,jk) ) 
+               zetot2 = 30. / ( 30.0 + etot_ndcy(ji,jj,jk) )
 
                xfracal(ji,jj,jk) = caco3r * MIN( zlim1, zlim2, zlim3 )                  &
@@ -375,4 +376,5 @@
          &      xlimbac (jpi,jpj,jpk), xlimbacl(jpi,jpj,jpk),       &
          &      concnfe (jpi,jpj,jpk), concdfe (jpi,jpj,jpk),       &
+         &      xqfuncfecn(jpi,jpj,jpk), xqfuncfecd(jpi,jpj,jpk),   &
          &      xlimsi  (jpi,jpj,jpk), STAT=p4z_lim_alloc )
       !
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlys.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlys.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlys.F90	(revision 14276)
@@ -57,5 +57,5 @@
       !
       INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zdispot, zfact, zcalcon
+      REAL(wp) ::   zdispot, zfact, zcalcon, zkd
       REAL(wp) ::   zomegaca, zexcess, zexcess0
       CHARACTER (len=25) ::   charout
@@ -65,6 +65,6 @@
       IF( ln_timing )  CALL timing_start('p4z_lys')
       !
-      zco3    (:,:,:) = 0.
-      zcaldiss(:,:,:) = 0.
+      zco3    (:,:,:) = 0._wp
+      zcaldiss(:,:,:) = 0._wp     ;     zco3sat(:,:,:) = 0._wp
       zhinit  (:,:,:) = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn )
       !
@@ -106,8 +106,14 @@
                excess(ji,jj,jk) = 1._wp - zomegaca
                zexcess0 = MAX( 0., excess(ji,jj,jk) )
-               zexcess  = zexcess0**nca
+               IF (zomegaca < 0.8) THEN
+                  zexcess  = zexcess0**nca
+                 ! AMOUNT CACO3 THAT RE-ENTERS SOLUTION
+                  zdispot = kdca * zexcess * trb(ji,jj,jk,jpcal)
+               ELSE
+                  zkd = kdca * 0.1**(nca - 0.4)
+                  zexcess  = zexcess0**0.4
+                  zdispot = zkd * zexcess * trb(ji,jj,jk,jpcal)
+               ENDIF
 
-               ! AMOUNT CACO3 THAT RE-ENTERS SOLUTION
-               zdispot = kdca * zexcess * trb(ji,jj,jk,jpcal)
                ! CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
                ! AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 14276)
@@ -89,5 +89,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrarem, zgraref, zgrapoc, zgrapof, zgrabsi
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigrem, zgramigref, zgramigpoc, zgramigpof
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zstrn, zgramigbsi
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigbsi
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zw3d, zz2ligprod
       !!---------------------------------------------------------------------
@@ -236,5 +236,5 @@
                zgraztotn = zgrazdc * quotad(ji,jj,jk) + zgrazz + zgraznc * quotan(ji,jj,jk)   &
                &   + zgrazpoc + zgrazffep + zgrazffeg
-               zgraztotf = zgrazdf + zgraznf + zgrazz * ferat3 + zgrazpof + zgrazfffp + zgrazfffg
+               zgraztotf = zgrazdf + zgraznf + zgrazz * feratz + zgrazpof + zgrazfffp + zgrazfffg
 
                ! Total grazing ( grazing by microzoo is already computed in p4zmicro )
@@ -251,5 +251,5 @@
                zgrasratf =  ( zgraztotf + rtrn )/ ( zgraztotc + rtrn )
                zgrasratn =  ( zgraztotn + rtrn )/ ( zgraztotc + rtrn )
-               zepshert  = MIN( 1., zgrasratn, zgrasratf / ferat3)
+               zepshert  = MIN( 1., zgrasratn, zgrasratf / feratm)
                zbeta     = MAX(0., (epsher2 - epsher2min) )
                ! Food quantity deprivation of GGE
@@ -298,5 +298,5 @@
                zprcaca = part2 * zprcaca
                tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrazcal - zprcaca
-               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * ( zgrazcal + zprcaca )
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * ( zgrazcal - zprcaca )
                tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) - zgrazcal + zprcaca
  
@@ -305,8 +305,8 @@
                zgrarem(ji,jj,jk) = zgraztotc * ( 1. - zepsherv - unass2 ) &
                &         + ( 1. - epsher2 - unass2 ) / ( 1. - epsher2 ) * ztortz
-               zgraref(ji,jj,jk) = zgraztotc * MAX( 0. , ( 1. - unass2 ) * zgrasratf - ferat3 * zepsherv )    &
-               &         + ferat3 * ( ( 1. - epsher2 - unass2 ) /( 1. - epsher2 ) * ztortz )
+               zgraref(ji,jj,jk) = zgraztotc * MAX( 0. , ( 1. - unass2 ) * zgrasratf - feratm * zepsherv )    &
+               &         + feratm * ( ( 1. - epsher2 - unass2 ) /( 1. - epsher2 ) * ztortz )
                zgrapoc(ji,jj,jk) = zgraztotc * unass2 + zmortzgoc
-               zgrapof(ji,jj,jk) = zgraztotf * unass2 + ferat3 * zmortzgoc
+               zgrapof(ji,jj,jk) = zgraztotf * unass2 + feratm * zmortzgoc
             END DO 
          END DO
@@ -320,23 +320,7 @@
          ALLOCATE( zgramigrem(jpi,jpj), zgramigref(jpi,jpj), zgramigpoc(jpi,jpj), zgramigpof(jpi,jpj) )
          ALLOCATE( zgramigbsi(jpi,jpj) )
-         ALLOCATE( zstrn(jpi,jpj) )
          zgramigrem(:,:) = 0.0    ;   zgramigref(:,:) = 0.0
          zgramigpoc(:,:) = 0.0    ;   zgramigpof(:,:) = 0.0
          zgramigbsi(:,:) = 0.0
-
-         ! compute the day length depending on latitude and the day
-         zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
-         zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
-
-         ! day length in hours
-         zstrn(:,:) = 0.
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
-               zargu = MAX( -1., MIN(  1., zargu ) )
-               zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
-               zstrn(ji,jj) = MIN(0.75, MAX( 0.25, zstrn(ji,jj) / 24.) )
-            END DO
-         END DO
 
         ! Compute the amount of materials that will go into vertical migration
@@ -347,5 +331,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zmigreltime = (1. - zstrn(ji,jj))
+                  zmigreltime = (1. - strn(ji,jj))
                   IF ( gdept_n(ji,jj,jk) <= heup(ji,jj) ) THEN
                      zgramigrem(ji,jj) = zgramigrem(ji,jj) + xfracmig * zgrarem(ji,jj,jk) * (1. - zmigreltime )    &
@@ -389,6 +373,4 @@
          ! ------------------------------
          DEALLOCATE( zgramigrem, zgramigref, zgramigpoc, zgramigpof, zgramigbsi )
-         DEALLOCATE( zstrn )
-
       ! End of the ln_dvm_meso part
       ENDIF
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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 14276)
@@ -187,5 +187,5 @@
                zgrasratf = ( zgraztotf + rtrn ) / ( zgraztotc + rtrn )
                zgrasratn = ( zgraztotn + rtrn ) / ( zgraztotc + rtrn )
-               zepshert  =  MIN( 1., zgrasratn, zgrasratf / ferat3)
+               zepshert  =  MIN( 1., zgrasratn, zgrasratf / feratz)
                zbeta     = MAX(0., (epsher - epshermin) )
                ! Food quantity deprivation of the GGE
@@ -196,5 +196,5 @@
                zepsherv  = zepsherf * zepshert * zepsherq
                ! Excretion of Fe
-               zgrafer   = zgraztotc * MAX( 0. , ( 1. - unass ) * zgrasratf - ferat3 * zepsherv ) 
+               zgrafer   = zgraztotc * MAX( 0. , ( 1. - unass ) * zgrasratf - feratz * zepsherv ) 
                ! Excretion of C, N, P
                zgrarem   = zgraztotc * ( 1. - zepsherv - unass )
@@ -236,5 +236,5 @@
                prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zmortz
                conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazpoc
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz - zgrazpof
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + feratz * zmortz - zgrazpof
                !
                ! Calcite remineralization due to zooplankton activity
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmort.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 14276)
@@ -25,5 +25,4 @@
    REAL(wp), PUBLIC ::   wchln    !: Quadratic mortality rate of nanophytoplankton
    REAL(wp), PUBLIC ::   wchld    !: Quadratic mortality rate of diatoms
-   REAL(wp), PUBLIC ::   wchldm   !: Maximum quadratic mortality rate of diatoms
    REAL(wp), PUBLIC ::   mpratn   !: Linear mortality rate of nanophytoplankton
    REAL(wp), PUBLIC ::   mpratd   !: Linear mortality rate of diatoms
@@ -63,7 +62,7 @@
       !!---------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zsizerat, zcompaph
+      REAL(wp) ::   zcompaph
       REAL(wp) ::   zfactfe, zfactch, zprcaca, zfracal
-      REAL(wp) ::   ztortp , zrespp , zmortp 
+      REAL(wp) ::   ztortp , zrespp , zmortp, zlim1, zlim2 
       CHARACTER (len=25) ::   charout
       !!---------------------------------------------------------------------
@@ -76,14 +75,10 @@
             DO ji = 1, jpi
                zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e-8 ), 0.e0 )
-               !     When highly limited by macronutrients, very small cells 
-               !     dominate the community. As a consequence, aggregation
-               !     due to turbulence is negligible. Mortality is also set
-               !     to 0
-               zsizerat = MIN(1., MAX( 0., (quotan(ji,jj,jk) - 0.2) / 0.3) ) * trb(ji,jj,jk,jpphy)
-
                ! Quadratic mortality of nano due to aggregation during
                ! blooms (Doney et al. 1996)
                ! -----------------------------------------------------
-               zrespp = wchln * 1.e6 * xstep * xdiss(ji,jj,jk) * zcompaph * zsizerat 
+               zlim2   = xlimphy(ji,jj,jk) * xlimphy(ji,jj,jk)
+               zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) * trb(ji,jj,jk,jpphy)
+               zrespp  = wchln * 1.e6 * xstep * zlim1 * xdiss(ji,jj,jk) * zcompaph
 
                ! Phytoplankton linear mortality
@@ -91,5 +86,5 @@
                ! extinction of nanophyto in highly limited areas
                ! ----------------------------------------------------
-               ztortp = mpratn * xstep * zcompaph / ( xkmort + trb(ji,jj,jk,jpphy) ) * zsizerat
+               ztortp = mpratn * xstep * zcompaph / ( xkmort + trb(ji,jj,jk,jpphy) ) * trb(ji,jj,jk,jpphy)
 
                zmortp = zrespp + ztortp
@@ -170,5 +165,5 @@
                zlim2   = xlimdia(ji,jj,jk) * xlimdia(ji,jj,jk)
                zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) 
-               zrespp2 = 1.e6 * xstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
+               zrespp2 = 1.e6 * xstep * wchld * zlim1 * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
 
                ! Phytoplankton linear mortality
@@ -228,5 +223,5 @@
       INTEGER ::   ios   ! Local integer
       !
-      NAMELIST/namp4zmort/ wchln, wchld, wchldm, mpratn, mpratd
+      NAMELIST/namp4zmort/ wchln, wchld, mpratn, mpratd
       !!----------------------------------------------------------------------
       !
@@ -249,5 +244,4 @@
          WRITE(numout,*) '      quadratic mortality of phytoplankton        wchln  =', wchln
          WRITE(numout,*) '      maximum quadratic mortality of diatoms      wchld  =', wchld
-         WRITE(numout,*) '      maximum quadratic mortality of diatoms      wchldm =', wchldm
          WRITE(numout,*) '      phytoplankton mortality rate                mpratn =', mpratn
          WRITE(numout,*) '      Diatoms mortality rate                      mpratd =', mpratd
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zopt.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 14276)
@@ -67,5 +67,5 @@
       REAL(wp), DIMENSION(jpi,jpj    ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4
       REAL(wp), DIMENSION(jpi,jpj    ) :: zqsr100, zqsr_corr
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpar, ze0, ze1, ze2, ze3, zchl3d
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze0, ze1, ze2, ze3, zchl3d
       !!---------------------------------------------------------------------
       !
@@ -206,5 +206,5 @@
                  heup(ji,jj) = gdepw_n(ji,jj,jk+1)     ! Euphotic layer depth
               ENDIF
-              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 )  THEN
+              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.10 )  THEN
                  heup_01(ji,jj) = gdepw_n(ji,jj,jk+1)  ! Euphotic layer depth (light level definition)
               ENDIF
@@ -235,6 +235,6 @@
       END DO
       !
-      emoy(:,:,:) = etot(:,:,:)       ! PAR
-      zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle 
+      emoy(:,:,:)  = etot(:,:,:)       ! PAR
+      etotm(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle 
       !
       DO jk = 1, nksrp
@@ -244,5 +244,5 @@
                   z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
                   emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep
-                  zpar (ji,jj,jk) = zetmp2(ji,jj) * z1_dep
+                  etotm(ji,jj,jk) = zetmp2(ji,jj) * z1_dep
                ENDIF
             END DO
@@ -313,5 +313,5 @@
         IF( knt == nrdttrc ) THEN
            IF( iom_use( "Heup"  ) ) CALL iom_put( "Heup" , heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
-           IF( iom_use( "PARDM" ) ) CALL iom_put( "PARDM", zpar(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
+           IF( iom_use( "PARDM" ) ) CALL iom_put( "PARDM", etotm(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
            IF( iom_use( "PAR"   ) ) CALL iom_put( "PAR"  , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
         ENDIF
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 14276)
@@ -67,13 +67,12 @@
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zsilfac, znanotot, zdiattot, zconctemp, zconctemp2
-      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap, zlim, zsilfac2, zsiborn
+      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap, zlim, zsiborn
       REAL(wp) ::   zprod, zproreg, zproreg2, zprochln, zprochld
       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) ::   zratiosi, zmaxsi, zlimfac, zsizetmp, zfecnm, zfecdm
+      REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
       CHARACTER (len=25) :: charout
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zw2d
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zstrn
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprmaxn,zprmaxd
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadd, zysopt  
@@ -94,4 +93,5 @@
       zprbio  (:,:,:) = 0._wp ; zprchld (:,:,:) = 0._wp ; zprchln (:,:,:) = 0._wp 
       zmxl_fac(:,:,:) = 0._wp ; zmxl_chl(:,:,:) = 0._wp 
+      consfe3 (:,:,:) = 0._wp
 
       ! Computation of the maximimum production. Based on a Q10 description
@@ -100,19 +100,4 @@
       zprmaxn(:,:,:) = 0.8_wp * r1_rday * tgfunc(:,:,:)
       zprmaxd(:,:,:) = zprmaxn(:,:,:)
-
-      ! compute the day length depending on latitude and the day
-      ! Astronomical parameterization taken from HAMOCC3
-      zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
-      zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
-
-      ! day length in hours
-      zstrn(:,:) = 0.
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
-            zargu = MAX( -1., MIN(  1., zargu ) )
-            zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
-         END DO
-      END DO
 
       ! Impact of the day duration and light intermittency on phytoplankton growth
@@ -127,5 +112,5 @@
             DO ji = 1, jpi
                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                  zval = MAX( 1., zstrn(ji,jj) )
+                  zval = MAX( 1., strn(ji,jj) )
                   IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
                      zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
@@ -225,7 +210,7 @@
                    ! Si/C of diatoms
                    ! ------------------------
-                   ! Si/C increases with iron stress and silicate availability (zsilfac)
+                   ! Si/C increases with iron stress and silicate availability
                    ! Si/C is arbitrariliy increased for very high Si concentrations
-                   ! to mimic the very high ratios observed in the Southern Ocean (zsilfac2)
+                   ! to mimic the very high ratios observed in the Southern Ocean (zsilfac)
                    ! A parameterization derived from Flynn (2003) is used for the control
                    ! when Si is not limiting which is similar to the parameterisation
@@ -236,15 +221,15 @@
                  zsiborn = trb(ji,jj,1,jpsil)**3
                  IF (gphit(ji,jj) < -30.0 ) THEN
-                   zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
+                   zsilfac = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
                  ELSE
-                   zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 )
+                   zsilfac = 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 = 1.0 - trb(ji,jj,jk,jpdsi) / ( trb(ji,jj,jk,jpdia) + rtrn ) / ( zsilfac * 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
+                    zysopt(ji,jj,jk) = zlim * zsilfac * grosip * 1.0 * zmaxsi
                  ELSE
-                    zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zsilim**0.7 * zmaxsi
+                    zysopt(ji,jj,jk) = zlim * zsilfac * grosip * 1.0 * zsilim**0.7 * zmaxsi
                  ENDIF
               ENDIF
@@ -291,7 +276,8 @@
                   ! do not suggest it for accimated cells. Uptake is
                   ! downregulated when the quota is close to the maximum quota
-                  zratio = 1.0 - MIN(1.0,trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * fecnm + rtrn ) )
+                  zfecnm = xqfuncfecn(ji,jj,jk) + ( fecnm - xqfuncfecn(ji,jj,jk) ) * ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) )
+                  zratio = 1.0 - MIN(1.0,trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * zfecnm + 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) )  &
+                  zprofen(ji,jj,jk) = zfecnm * zprmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
                   &          * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn + xnanono3(ji,jj,jk)  &
                   &          + xnanonh4(ji,jj,jk) ) * (1. - xnanofer(ji,jj,jk) ) )   &
@@ -317,7 +303,8 @@
                   ! do not suggest it for accimated cells. Uptake is
                   ! downregulated when the quota is close to the maximum quota
-                  zratio = 1.0 - MIN(1.0, trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * fecdm + rtrn ) )
+                  zfecdm = xqfuncfecd(ji,jj,jk) + ( fecdm - xqfuncfecd(ji,jj,jk) ) * ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) )
+                  zratio = 1.0 - MIN(1.0, trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * zfecdm + 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) )  &
+                  zprofed(ji,jj,jk) = zfecdm * 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) ) )   &
@@ -339,9 +326,5 @@
                   zprod    = rday * zprorcan(ji,jj,jk) * zprchln(ji,jj,jk) * xlimphy(ji,jj,jk)
                   zprochln = chlcmin * 12. * zprorcan (ji,jj,jk)
-
-                  ! The maximum reachable Chl quota is modulated by temperature
-                  ! following Geider (1987)
-                  chlcnm_n   = MIN ( chlcnm, ( chlcnm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem))) * (1. - 1.14 / 43.4 * 20.))
-                  zprochln = zprochln + (chlcnm_n-chlcmin) * 12. * zprod / &
+                  zprochln = zprochln + (chlcnm - chlcmin) * 12. * zprod / &
                                         & (  zpislopeadn(ji,jj,jk) * znanotot +rtrn)
 
@@ -350,9 +333,5 @@
                   zprod    = rday * zprorcad(ji,jj,jk) * zprchld(ji,jj,jk) * xlimdia(ji,jj,jk)
                   zprochld = chlcmin * 12. * zprorcad(ji,jj,jk)
-
-                  ! The maximum reachable Chl quota is modulated by temperature
-                  ! following Geider (1987)
-                  chlcdm_n   = MIN ( chlcdm, ( chlcdm / (1. - 1.14 / 43.4 * tsn(ji,jj,jk,jp_tem))) * (1. - 1.14 / 43.4 * 20.))
-                  zprochld = zprochld + (chlcdm_n-chlcmin) * 12. * zprod / &
+                  zprochld = zprochld + (chlcdm - chlcmin) * 12. * zprod / &
                                         & ( zpislopeadd(ji,jj,jk) * zdiattot +rtrn )
 
@@ -386,6 +365,8 @@
                  &                   + ( o2ut + o2nit ) * ( zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk) )
                  !
-                 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,jpfer) = tra(ji,jj,jk,jpfer) - ( texcretn * zprofen(ji,jj,jk)    &
+                 &                   + texcretd * zprofed(ji,jj,jk) )
+                 consfe3(ji,jj,jk) = ( texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) )     &
+                 &                    * 75.0 / ( rtrn + ( plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * trb(ji,jj,jk,jpfer) ) / rfact2
                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk)   &
                  &                   * rfact2 * trb(ji,jj,jk,jpdia)
@@ -401,5 +382,5 @@
      ! when ln_ligand is set to .true. in the namelist. Ligand uptake is small 
      ! and based on the FeL model by Morel et al. (2008) and on the study of
-     ! Shaked and Lis (2012)
+     ! Shaked et al. (2020)
      ! -------------------------------------------------------------------------
      IF( ln_ligand ) THEN
@@ -412,8 +393,8 @@
                     zfeup    = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)
                     tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp    &
-                    &       - zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
+                    &       - zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * lthet
                     zpligprod1(ji,jj,jk) = zdocprod * ldocp
                     zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) &
-                    &                      + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
+                    &                      + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * lthet
                  ENDIF
               END DO
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zrem.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zrem.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zrem.F90	(revision 14276)
@@ -35,5 +35,4 @@
    REAL(wp), PUBLIC ::   xremikn    !: remineralisation rate of DON (p5z) 
    REAL(wp), PUBLIC ::   xremikp    !: remineralisation rate of DOP (p5z) 
-   REAL(wp), PUBLIC ::   xremik     !: remineralisation rate of DOC (p4z) 
    REAL(wp), PUBLIC ::   nitrif     !: NH4 nitrification rate 
    REAL(wp), PUBLIC ::   xsirem     !: remineralisation rate of biogenic silica
@@ -71,10 +70,8 @@
       REAL(wp) ::   zremik, zremikc, zremikn, zremikp, zsiremin, zfact 
       REAL(wp) ::   zsatur, zsatur2, znusil, znusil2, zdep, zdepmin, zfactdep
-      REAL(wp) ::   zbactfer, zolimit, zonitr, zrfact2
-      REAL(wp) ::   zammonic, zoxyremc, zoxyremn, zoxyremp
-      REAL(wp) ::   zosil, ztem, zdenitnh4, zolimic, zolimin, zolimip, zdenitrn, zdenitrp
+      REAL(wp) ::   zbactfer, zonitr
+      REAL(wp) ::   zammonic, zoxyremc, zosil, ztem, zdenitnh4, zolimic
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj    ) :: ztempbac
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepbac, zolimi, zdepprod, zfacsi, zfacsib, zdepeff, zfebact
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepbac, zolimi, zfacsi, zfacsib, zdepeff, zfebact
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
       !!---------------------------------------------------------------------
@@ -83,7 +80,5 @@
       !
       ! Initialisation of arrays
-      zdepprod(:,:,:) = 1._wp
       zdepeff (:,:,:) = 0.3_wp
-      ztempbac(:,:)   = 0._wp
       zfacsib(:,:,:)  = xsilab / ( 1.0 - xsilab )
       zfebact(:,:,:)  = 0._wp
@@ -99,12 +94,8 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zdep = MAX( hmld(ji,jj), heup(ji,jj) )
-               IF( gdept_n(ji,jj,jk) < zdep ) THEN
-                  zdepbac(ji,jj,jk) = MIN( 0.7 * ( trb(ji,jj,jk,jpzoo) + 2.* trb(ji,jj,jk,jpmes) ), 4.e-6 )
-                  ztempbac(ji,jj)   = zdepbac(ji,jj,jk)
-               ELSE
+               zdep = MAX( hmld(ji,jj), heup_01(ji,jj) )
+               zdepbac(ji,jj,jk) = 0.6 * ( MAX(0.0, trb(ji,jj,jk,jpzoo) + trb(ji,jj,jk,jpmes) ) * 1.0E6 )**0.6 * 1.E-6
+               IF( gdept_n(ji,jj,jk) >= zdep ) THEN
                   zdepmin = MIN( 1., zdep / gdept_n(ji,jj,jk) )
-                  zdepbac (ji,jj,jk) = zdepmin**0.683 * ztempbac(ji,jj)
-                  zdepprod(ji,jj,jk) = zdepmin**0.273
                   zdepeff (ji,jj,jk) = zdepeff(ji,jj,jk) * zdepmin**0.3
                ENDIF
@@ -113,92 +104,52 @@
       END DO
 
-      IF( ln_p4z ) THEN ! Standard PISCES code
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ! DOC ammonification. Depends on a limitation term of the bacterial activity
-                  ! and on the implicit bacteria concentration
-                  ! --------------------------------------------------------------------------
-                  zremik = xremik * xstep / 1.e-6 * xlimbac(ji,jj,jk) * zdepbac(ji,jj,jk) 
-                  zremik = MAX( zremik, 2.74e-4 * xstep )
-
-                  ! Ammonification in oxic waters with oxygen consumption
-                  ! -----------------------------------------------------
-                  zolimit = zremik * ( 1.- nitrfac(ji,jj,jk) ) * trb(ji,jj,jk,jpdoc) 
-                  zolimi(ji,jj,jk) = MIN( ( trb(ji,jj,jk,jpoxy) - rtrn ) / o2ut, zolimit ) 
-
-                  ! Ammonification in suboxic waters with denitrification
-                  ! -----------------------------------------------------
-                  zammonic = zremik * nitrfac(ji,jj,jk) * trb(ji,jj,jk,jpdoc)
-                  denitr(ji,jj,jk)  = zammonic * ( 1. - nitrfac2(ji,jj,jk) )
-                  denitr(ji,jj,jk)  = MIN( ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenit, denitr(ji,jj,jk) )
-
-                  ! Ammonification in waters depleted in O2 and NO3 based on 
-                  ! other redox processes
-                  ! --------------------------------------------------------
-                  zoxyremc          = zammonic - denitr(ji,jj,jk)
-                  !
-                  zolimi (ji,jj,jk) = MAX( 0.e0, zolimi (ji,jj,jk) )
-                  denitr (ji,jj,jk) = MAX( 0.e0, denitr (ji,jj,jk) )
-                  zoxyremc          = MAX( 0.e0, zoxyremc )
-                  ! Update of the TRA arrays
-                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk) + zoxyremc
-                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk) + zoxyremc
-                  tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - denitr (ji,jj,jk) * rdenit
-                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zolimi (ji,jj,jk) - denitr(ji,jj,jk) - zoxyremc
-                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - zolimi (ji,jj,jk) * o2ut
-                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zolimi (ji,jj,jk) + denitr(ji,jj,jk) + zoxyremc
-                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zolimi(ji,jj,jk) + zoxyremc    &
-                  &                     + ( rdenit + 1.) * denitr(ji,jj,jk) )
-               END DO
-            END DO
-         END DO
-      ELSE ! PISCES-QUOTA
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ! DOC ammonification. Depends on a limitation term of the bacterial activity
-                  ! and on the implicit bacteria concentration
-                  ! ---------------------------------------------------------------
-                  zremik = xstep / 1.e-6 * MAX(0.01, xlimbac(ji,jj,jk)) * zdepbac(ji,jj,jk) 
-                  zremik = MAX( zremik, 2.74e-4 * xstep / xremikc )
-
-                  zremikc = xremikc * zremik
-                  zremikn = xremikn / xremikc
-                  zremikp = xremikp / xremikc
-
-                  ! Ammonification in oxic waters with oxygen consumption
-                  ! -----------------------------------------------------
-                  zolimit = zremikc * ( 1.- nitrfac(ji,jj,jk) ) * trb(ji,jj,jk,jpdoc) 
-                  zolimic = MAX( 0.e0, MIN( ( trb(ji,jj,jk,jpoxy) - rtrn ) / o2ut, zolimit ) ) 
-                  zolimi(ji,jj,jk) = zolimic
-                  zolimin = zremikn * zolimic * trb(ji,jj,jk,jpdon) / ( trb(ji,jj,jk,jpdoc) + rtrn )
-                  zolimip = zremikp * zolimic * trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdoc) + rtrn ) 
-
-                  ! Ammonification in suboxic waters with denitrification
-                  ! -------------------------------------------------------
-                  zammonic = zremikc * nitrfac(ji,jj,jk) * trb(ji,jj,jk,jpdoc)
-                  denitr(ji,jj,jk)  = zammonic * ( 1. - nitrfac2(ji,jj,jk) )
-                  denitr(ji,jj,jk)  = MAX(0., MIN(  ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenit, denitr(ji,jj,jk) ) )
-                  zoxyremc          = MAX(0., zammonic - denitr(ji,jj,jk))
-                  zdenitrn  = zremikn * denitr(ji,jj,jk) * trb(ji,jj,jk,jpdon) / ( trb(ji,jj,jk,jpdoc) + rtrn )
-                  zdenitrp  = zremikp * denitr(ji,jj,jk) * trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdoc) + rtrn )
-                  zoxyremn  = zremikn * zoxyremc * trb(ji,jj,jk,jpdon) / ( trb(ji,jj,jk,jpdoc) + rtrn )
-                  zoxyremp  = zremikp * zoxyremc * trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdoc) + rtrn )
-                  ! Update of the TRA arrays
-                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zolimip + zdenitrp + zoxyremp
-                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zolimin + zdenitrn + zoxyremn
-                  tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - denitr(ji,jj,jk) * rdenit
-                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zolimic - denitr(ji,jj,jk) - zoxyremc
-                  tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) - zolimin - zdenitrn - zoxyremn
-                  tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) - zolimip - zdenitrp - zoxyremp
-                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - zolimic * o2ut
-                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zolimic + denitr(ji,jj,jk) + zoxyremc
-                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zolimin + zoxyremn + ( rdenit + 1.) * zdenitrn )
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! DOC ammonification. Depends on a limitation term of the bacterial activity
+               ! and on the implicit bacteria concentration
+               ! --------------------------------------------------------------------------
+               zremik = xstep / 1.e-6 * xlimbac(ji,jj,jk) * zdepbac(ji,jj,jk) 
+               zremik = MAX( zremik, 2.74e-4 * xstep / xremikc )
+               zremikc = xremikc * zremik
+
+               ! Ammonification in oxic waters with oxygen consumption
+               ! -----------------------------------------------------
+               zolimic = zremikc * ( 1.- nitrfac(ji,jj,jk) ) * trb(ji,jj,jk,jpdoc)
+               zolimic = MAX( 0.e0, MIN( ( trb(ji,jj,jk,jpoxy) - rtrn ) / o2ut, zolimic ) )
+               zolimi(ji,jj,jk) = zolimic
+
+               ! Ammonification in suboxic waters with denitrification
+               ! -----------------------------------------------------
+               zammonic = zremikc * nitrfac(ji,jj,jk) * trb(ji,jj,jk,jpdoc)
+               denitr(ji,jj,jk)  = zammonic * ( 1. - nitrfac2(ji,jj,jk) )
+               denitr(ji,jj,jk)  = MAX(0., MIN(  ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenit, denitr(ji,jj,jk) ) )
+
+               ! Ammonification in waters depleted in O2 and NO3 based on 
+               ! other redox processes
+               ! --------------------------------------------------------
+               zoxyremc          = MAX(0., zammonic - denitr(ji,jj,jk) )
+
+               ! Update of the TRA arrays
+               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - denitr(ji,jj,jk) * rdenit
+               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - ( zolimic + denitr(ji,jj,jk) + zoxyremc )
+               tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - zolimic * o2ut
+               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zolimic + denitr(ji,jj,jk) + zoxyremc
+               IF ( ln_p4z ) THEN ! PISCES-std
+                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zolimic + denitr(ji,jj,jk) + zoxyremc
+                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zolimic + denitr(ji,jj,jk) + zoxyremc
+                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zolimic + zoxyremc + ( rdenit + 1.) * denitr(ji,jj,jk) )
+               ELSE  ! PISCES-QUOTA (p5z)
+                  zremikn = xremikn / xremikc * trb(ji,jj,jk,jpdon) / ( trb(ji,jj,jk,jpdoc) + rtrn )
+                  zremikp = xremikp / xremikc * trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdoc) + rtrn )
+                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zremikp * ( zolimic + denitr(ji,jj,jk) + zoxyremc )
+                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zremikn * ( zolimic + denitr(ji,jj,jk) + zoxyremc )
+                  tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) - zremikn * ( zolimic + denitr(ji,jj,jk) + zoxyremc )
+                  tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) - zremikp * ( zolimic + denitr(ji,jj,jk) + zoxyremc )
+                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zremikn * ( zolimic + zoxyremc + ( rdenit + 1.) * denitr(ji,jj,jk) )
+               ENDIF
+            END DO
+         END DO
+      END DO
 
       DO jk = 1, jpkm1
@@ -211,5 +162,5 @@
                &         / ( 1.+ emoy(ji,jj,jk) ) * ( 1. + fr_i(ji,jj) * emoy(ji,jj,jk) ) 
                zdenitnh4 = nitrif * xstep * trb(ji,jj,jk,jpnh4) * nitrfac(ji,jj,jk)
-               zdenitnh4 = MIN(  ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenita, zdenitnh4 ) 
+               zdenitnh4 = MAX(0., MIN(  ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenita, zdenitnh4 ) ) 
                ! Update of the tracers trends
                ! ----------------------------
@@ -217,5 +168,5 @@
                tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zonitr - rdenita * zdenitnh4
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2nit * zonitr
-               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2 * rno3 * zonitr + rno3 * ( rdenita - 1. ) * zdenitnh4
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * rno3 * zonitr + rno3 * ( rdenita - 1. ) * zdenitnh4
             END DO
          END DO
@@ -236,15 +187,14 @@
                ! studies (especially at Papa) have shown this uptake to be significant
                ! ---------------------------------------------------------------------
-               zbactfer = feratb *  rfact2 * 0.6_wp / rday * tgfunc(ji,jj,jk) * xlimbacl(ji,jj,jk)     &
-                  &              * trb(ji,jj,jk,jpfer) / ( xkferb + trb(ji,jj,jk,jpfer) )    &
-                  &              * zdepeff(ji,jj,jk) * zdepbac(ji,jj,jk)
+               zbactfer = feratb * 0.6_wp * xstep * tgfunc(ji,jj,jk) * xlimbacl(ji,jj,jk) * trb(ji,jj,jk,jpfer)    &
+                  &       / ( xkferb + trb(ji,jj,jk,jpfer) ) * zdepeff(ji,jj,jk) * zdepbac(ji,jj,jk)
 
                ! Only the transfer of iron from its dissolved form to particles
                ! is treated here. The GGE of bacteria supposed to be equal to 
                ! 0.33. This is hard-coded. 
-               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zbactfer*0.33
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zbactfer*0.25
-               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zbactfer*0.08
-               zfebact(ji,jj,jk)   = zbactfer * 0.33
+               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zbactfer*0.18
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zbactfer*0.15
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zbactfer*0.03
+               zfebact(ji,jj,jk)   = zbactfer * 0.18
                blim(ji,jj,jk)      = xlimbacl(ji,jj,jk)  * zdepbac(ji,jj,jk) / 1.e-6
             END DO
@@ -304,5 +254,4 @@
 
       IF( knt == nrdttrc ) THEN
-          zrfact2 = 1.e3 * rfact2r
           ALLOCATE( zw3d(jpi,jpj,jpk) )
           zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
@@ -321,5 +270,5 @@
           ENDIF
           IF( iom_use( "FEBACT" ) )  THEN
-               zw3d(:,:,:) = zfebact(:,:,:) * 1E9 * tmask(:,:,:) * zrfact2   ! Bacterial iron consumption
+               zw3d(:,:,:) = zfebact(:,:,:) * 1E9 * tmask(:,:,:) * zfact   ! Bacterial iron consumption
                CALL iom_put( "FEBACT" , zw3d )
           ENDIF
@@ -345,5 +294,5 @@
       !!
       !!----------------------------------------------------------------------
-      NAMELIST/nampisrem/ xremik, nitrif, xsirem, xsiremlab, xsilab, feratb, xkferb, & 
+      NAMELIST/nampisrem/ nitrif, xsirem, xsiremlab, xsilab, feratb, xkferb, & 
          &                xremikc, xremikn, xremikp
       INTEGER :: ios                 ! Local integer output status for namelist read
@@ -367,5 +316,5 @@
          WRITE(numout,*) '   Namelist parameters for remineralization, nampisrem'
          IF( ln_p4z ) THEN
-            WRITE(numout,*) '      remineralization rate of DOC              xremik    =', xremik
+            WRITE(numout,*) '      remineralization rate of DOC              xremikc   =', xremikc
          ELSE
             WRITE(numout,*) '      remineralization rate of DOC              xremikc   =', xremikc
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsbc.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsbc.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsbc.F90	(revision 14276)
@@ -483,5 +483,5 @@
       !
       sedsilfrac = 0.03     ! percentage of silica loss in the sediments
-      sedcalfrac = 0.6      ! percentage of calcite loss in the sediments
+      sedcalfrac = 0.99     ! percentage of calcite loss in the sediments
       !
    END SUBROUTINE p4z_sbc_init
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 14276)
@@ -66,5 +66,5 @@
       REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep
       REAL(wp) ::  zwstpoc, zwstpon, zwstpop
-      REAL(wp) ::  ztrfer, ztrpo4s, ztrdp, zwdust, zmudia, ztemp
+      REAL(wp) ::  ztrpo4s, ztrdp, zwdust, zmudia, ztemp
       REAL(wp) ::  xdiano3, xdianh4
       !
@@ -73,7 +73,7 @@
       REAL(wp), DIMENSION(jpi,jpj    ) :: zwsbio3, zwsbio4
       REAL(wp), DIMENSION(jpi,jpj    ) :: zsedcal, zsedsi, zsedc
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsoufer, zlight
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zpdep
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zsidep, zironice
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsoufer, zlight, ztrfer
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zpdep, zsidep
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zironice
       !!---------------------------------------------------------------------
       !
@@ -137,5 +137,5 @@
       IF( ln_dust ) THEN
          !                                              
-         ALLOCATE( zsidep(jpi,jpj), zpdep(jpi,jpj,jpk), zirondep(jpi,jpj,jpk) )
+         ALLOCATE( zsidep(jpi,jpj,jpk), zpdep(jpi,jpj,jpk), zirondep(jpi,jpj,jpk) )
 
          ! Iron, P and Si deposition at the surface
@@ -154,6 +154,10 @@
          ! and the solubility are hard coded
          ! ----------------------------------------------------------------
-         zsidep(:,:)   = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 28.1 
-         zpdep (:,:,1) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 31. / po4r 
+         ! Si crustal abundance is about 26.9% by mass and 7.5% is soluble
+         ! (see for instance, Tegen and Kohfeld, 2006).
+         zsidep(:,:,1) = 0.269 * 0.075 * dust(:,:) * rfact2 / e3t_n(:,:,1) / 28.1 
+         ! P Crustal abundance is about 0.1% by mass and about 10% of it soluble
+         ! (Paytan and McLaughlin, 2007).
+         zpdep (:,:,1) = 0.1 * 1.e-3 * dust(:,:) * rfact2 / e3t_n(:,:,1) / 31. / po4r 
 
          ! Iron solubilization of particles in the water column
@@ -168,12 +172,13 @@
          DO jk = 2, jpkm1
             zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -gdept_n(:,:,jk) / (250. * wdust) )
-            zpdep   (:,:,jk) = zirondep(:,:,jk) * 0.38 / po4r
+            zpdep   (:,:,jk) = zirondep(:,:,jk) * 1.e-3 / mfrac * 55.85 / 31.
+            zsidep  (:,:,jk) = zirondep(:,:,jk) * 0.269 / mfrac * 55.85 / 28.1
          END DO
 
          ! Solubilization of particles in the water column (Si, P, Fe)
-         tra(:,:,1,jpsil) = tra(:,:,1,jpsil) + zsidep  (:,:)
          DO jk = 1, jpkm1
             tra(:,:,jk,jppo4) = tra(:,:,jk,jppo4) + zpdep   (:,:,jk)
             tra(:,:,jk,jpfer) = tra(:,:,jk,jpfer) + zirondep(:,:,jk) 
+            tra(:,:,jk,jpsil) = tra(:,:,jk,jpsil) + zsidep  (:,:,jk)
          ENDDO
          ! 
@@ -299,5 +304,5 @@
                  zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
                    &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6
-                 zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2
+                 zbureff(ji,jj) = 0.013 + 0.13 * zflx**2 / ( 7.0 + zflx )**2
               ENDIF
             END DO
@@ -306,8 +311,7 @@
       ENDIF
 
-      ! Fraction of dSi that is remineralized in the sediments. This is 
-      ! set so that the burial in sediments equals the total input of Si
-      ! by rivers and dust (sedsilfrac)
-      ! ----------------------------------------------------------------
+      ! Fraction of dSi that is dissolved in the sediments. This fraction is  
+      ! set to a constant value in p4zsbc
+      ! --------------------------------------------------------------------
       IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac
 
@@ -331,5 +335,7 @@
          ! modeled. The amount of CaCO3 that dissolves in the sediments
          ! is computed using a metamodel constructed from Archer (1996)
-         ! ------------------------------------------------------------
+         ! A minimum set to sedcalfrac is preserved. This value is defined
+         ! in p4zsbc
+         ! ---------------------------------------------------------------
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -341,6 +347,6 @@
                tra(ji,jj,ikt,jpsil) = tra(ji,jj,ikt,jpsil) + zsiloss * zrivsil 
                !
-               zfactcal = MIN( excess(ji,jj,ikt), 0.2 )
-               zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) )
+               zfactcal = MAX(-0.1, MIN( excess(ji,jj,ikt), 0.2 ) )
+               zfactcal = 0.3 + 0.7 * MIN( 1., (0.1 + zfactcal) / ( 0.5 - zfactcal ) )
                zrivalk  = sedcalfrac * zfactcal
                tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0
@@ -394,10 +400,17 @@
                zws4 = zwsbio4(ji,jj) * zdep
                zws3 = zwsbio3(ji,jj) * zdep
+               ! Fraction that is permanently buried in the sediments
                zrivno3 = 1. - zbureff(ji,jj)
                zwstpoc = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3
+               ! Denitrification in the sediments
                zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 )
+               ! Fraction that is not denitrified
                z1pdenit = zwstpoc * zrivno3 - zpdenit
+               ! Oxic remineralization of organic matter in the sediments
                zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) )
+               ! The fraction that cannot be denitrified nor oxidized by O2
+               ! is released back to the water column as DOC
                tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit
+               ! Update of the tracers concentrations
                tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit
                tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit
@@ -408,4 +421,5 @@
                sdenit(ji,jj) = rdenit * zpdenit * e3t_n(ji,jj,ikt)
                zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t_n(ji,jj,ikt)
+               ! PISCES-QUOTA (p5z)
                IF( ln_p5z ) THEN
                   zwstpop              = trb(ji,jj,ikt,jpgop) * zws4 + trb(ji,jj,ikt,jppop) * zws3
@@ -421,8 +435,9 @@
       ! Nitrogen fixation process : light limitation of diazotrophy
       ! Small source of iron from particulate inorganic iron (photochemistry)
+      ! This is a purely adhoc param.
       !----------------------------------------------------------------------
       DO jk = 1, jpkm1
          zlight (:,:,jk) =  ( 1.- EXP( -etot_ndcy(:,:,jk) / diazolight ) ) * ( 1. - fr_i(:,:) ) 
-         zsoufer(:,:,jk) = zlight(:,:,jk) * 2E-11 / ( 2E-11 + biron(:,:,jk) )
+         zsoufer(:,:,jk) = zlight(:,:,jk) * 1E-10 / ( 1E-10 + biron(:,:,jk) )
       ENDDO
 
@@ -438,18 +453,15 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! Potential nitrogen fixation dependant on temperature and iron
+                  ! Potential nitrogen fixation dependant on temperature
                   ztemp = tsn(ji,jj,jk,jp_tem)
-                  zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
-                  xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) )
-                  xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
-                  zlim = ( 1.- xdiano3 - xdianh4 )
-                  ! Nitrogen fixation is almost fully halted when the N 
-                  ! limitation term (xdiano3+xdianh4) is > 0.9
-                  IF( zlim <= 0.1 )   zlim = 0.01
+                  zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) / rno3
+                  ! Nitrogen fixation is inhibited when enough NO3 and/or NH4
+                  zlim = ( 1.- xnanonh4(ji,jj,jk) - xnanono3(ji,jj,jk) )
                   zfact = zlim * rfact2
-                  ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
-                  ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) )
+                  ! Nitrogen fixation limitation by PO4 and Fe
+                  ztrfer(ji,jj,jk) = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
+                  ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( concnnh4 + trb(ji,jj,jk,jppo4) )
                   ztrdp = ztrpo4(ji,jj,jk)
-                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk)
+                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer(ji,jj,jk), ztrdp ) * zlight(ji,jj,jk)
                END DO
             END DO
@@ -460,20 +472,18 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! Potential nitrogen fixation dependant on temperature and iron
+                  ! Potential nitrogen fixation dependant on temperature
                   ztemp = tsn(ji,jj,jk,jp_tem)
                   zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
+                  ! Nitrogen fixation is inhibited when enough NO3 and/or NH4
                   xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) )
                   xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
                   zlim = ( 1.- xdiano3 - xdianh4 )
-
-                  ! Nitrogen fixation is almost fully halted when the N 
-                  ! limitation term (xdiano3+xdianh4) is > 0.9
-                  IF( zlim <= 0.1 )   zlim = 0.01
                   zfact = zlim * rfact2
-                  ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
+                  ! Nitrogen fixation limitation by PO4/DOP and Fe
+                  ztrfer(ji,jj,jk) = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
                   ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) )
                   ztrdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( 1E-6 + trb(ji,jj,jk,jpdop) ) * (1. - ztrpo4(ji,jj,jk))
                   ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk)
-                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk)
+                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer(ji,jj,jk), ztrdp ) * zlight(ji,jj,jk)
                END DO
             END DO
@@ -490,5 +500,5 @@
                   zfact = nitrpot(ji,jj,jk) * nitrfix
                   ! 1/3 of the diazotrophs growth is supposed to be excreted
-                  ! as NH4. 1/3 as DOC and the rest is routed POC and GOC as 
+                  ! as NH4. 1/3 as DOC and the rest is routed to POC/GOC as 
                   ! a result of mortality by predation. Completely adhoc param 
                   tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0
@@ -499,9 +509,11 @@
                   tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0
                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
-                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C
-                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0
-                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
-                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
-                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
+                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C at max
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0
+                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 2.0 / 3.0
+                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 1.0 / 3.0
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.005 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
+                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trb(ji,jj,jk,jppo4) ) &
+                  &                     * 0.001 * trb(ji,jj,jk,jpdoc) * xstep
               END DO
             END DO 
@@ -533,8 +545,8 @@
                   tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0
                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
-                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C
-                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0 
-                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
-                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
+                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C at max
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 
+                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 2.0 / 3.0
+                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 1.0 / 3.0
                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
               END DO
@@ -558,5 +570,6 @@
             IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * zfact ) ! Permanent burial of bSi in sediments
             IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * zfact ) ! Permanent burial of OC in sediments
-            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 ) ! Denitrification in the sediments
+!            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 ) ! Denitrification in the sediments
+            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", zdenit2d(:,:) ) ! Denitrification in the sediments
          ENDIF
       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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 14276)
@@ -64,4 +64,5 @@
       REAL(wp) ::  ztra
       CHARACTER (len=25) :: charout
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jp_pisces) :: ztrbbio
       !!---------------------------------------------------------------------
       !
@@ -102,4 +103,8 @@
          END DO
       ENDIF
+
+      DO jn = jp_pcs0, jp_pcs1              !   Store the tracer concentrations before entering PISCES
+         ztrbbio(:,:,:,jn) = trb(:,:,:,jn)
+      END DO
       !
       IF( ll_sbc ) CALL p4z_sbc( kt )   ! external sources of nutrients 
@@ -147,12 +152,22 @@
             tra(:,:,:,jn) = 0._wp
          END DO
-         ! Euler-forward temporal scheme
-         IF( ln_top_euler ) THEN
-            DO jn = jp_pcs0, jp_pcs1
-               trn(:,:,:,jn) = trb(:,:,:,jn)
-            END DO
-         ENDIF
+         !
       END DO
-
+      !
+#endif
+      !
+      ! If ln_sediment is set to .true. then the sediment module is called
+      IF( ln_sediment ) THEN 
+         !
+         CALL sed_model( kt )     !  Main program of Sediment model
+         !
+      ENDIF
+      !
+      !
+      DO jn = jp_pcs0, jp_pcs1
+         tra(:,:,:,jn) = ( trb(:,:,:,jn) - ztrbbio(:,:,:,jn) ) * rfactr
+         trb(:,:,:,jn) = ztrbbio(:,:,:,jn)
+         ztrbbio(:,:,:,jn) = 0._wp
+      END DO
       !
       IF( l_trdtrc ) THEN
@@ -161,18 +176,4 @@
          END DO
       END IF
-#endif
-      !
-      ! If ln_sediment is set to .true. then the sediment module is called
-      IF( ln_sediment ) THEN 
-         !
-         CALL sed_model( kt )     !  Main program of Sediment model
-         ! Eulor forward temporal scheme
-         IF( ln_top_euler ) THEN
-            DO jn = jp_pcs0, jp_pcs1
-               trn(:,:,:,jn) = trb(:,:,:,jn)
-            END DO
-         ENDIF
-         !
-      ENDIF
       !
       IF( lrst_trc )  CALL p4z_rst( kt, 'WRITE' )  !* Write PISCES informations in restart file 
@@ -199,6 +200,6 @@
       INTEGER :: ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale,    &
-         &                   ldocp, ldocz, lthet, no3rat3, po4rat3
+      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, feratz, feratm, wsbio2, wsbio2max,    &
+         &                wsbio2scale, ldocp, ldocz, lthet, no3rat3, po4rat3
          !
       NAMELIST/nampisdmp/ ln_pisdmp, nn_pisdmp
@@ -226,8 +227,9 @@
          WRITE(numout,*) '      half saturation constant for mortality    xkmort      =', xkmort 
          IF( ln_p5z ) THEN
-            WRITE(numout,*) '      N/C in zooplankton                        no3rat3     =', no3rat3
-            WRITE(numout,*) '      P/C in zooplankton                        po4rat3     =', po4rat3
-         ENDIF
-         WRITE(numout,*) '      Fe/C in zooplankton                       ferat3      =', ferat3
+            WRITE(numout,*) '      N/C in zooplankton                     no3rat3     =', no3rat3
+            WRITE(numout,*) '      P/C in zooplankton                     po4rat3     =', po4rat3
+         ENDIF
+         WRITE(numout,*) '      Fe/C in microzooplankton                  feratz      =', feratz
+         WRITE(numout,*) '      Fe/C in microzooplankton                  feratz      =', feratm
          WRITE(numout,*) '      Big particles sinking speed               wsbio2      =', wsbio2
          WRITE(numout,*) '      Big particles maximum sinking speed       wsbio2max   =', wsbio2max
@@ -313,4 +315,12 @@
          ENDIF
 
+         ! Read the Fe3 consumption term by phytoplankton
+         IF( iom_varid( numrtr, 'Consfe3', ldstop = .FALSE. ) > 0 ) THEN
+            CALL iom_get( numrtr, jpdom_autoglo, 'Consfe3' , consfe3(:,:,:)  )
+         ELSE
+            consfe3(:,:,:) = 0._wp
+         ENDIF
+
+
          ! Read the cumulative total flux. If not in the restart file, it is set to 0          
          IF( iom_varid( numrtr, 'tcflxcum', ldstop = .FALSE. ) > 0 ) THEN  ! cumulative total flux of carbon
@@ -326,4 +336,5 @@
             sized(:,:,:) = 1.
          ENDIF
+         sized(:,:,:) = MAX( 1.0, sized(:,:,:) )
          IF( iom_varid( numrtr, 'sizen', ldstop = .FALSE. ) > 0 ) THEN
             CALL iom_get( numrtr, jpdom_autoglo, 'sizen' , sizen(:,:,:)  )
@@ -331,4 +342,5 @@
             sizen(:,:,:) = 1.
          ENDIF
+         sizen(:,:,:) = MAX( 1.0, sizen(:,:,:) )
 
          ! PISCES-QUOTA specific part
@@ -338,4 +350,5 @@
             IF( iom_varid( numrtr, 'sizep', ldstop = .FALSE. ) > 0 ) THEN
                CALL iom_get( numrtr, jpdom_autoglo, 'sizep' , sizep(:,:,:)  )
+               sizep(:,:,:) = MAX( 1.0, sizep(:,:,:) )
             ELSE
                sizep(:,:,:) = 1.
@@ -353,4 +366,5 @@
          CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )    ! Si 1/2 saturation constant
          CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) ) ! Si max concentration
+         CALL iom_rstput( kt, nitrst, numrtw, 'Consfe3', consfe3(:,:,:) ) ! 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
@@ -377,7 +391,7 @@
       !
       REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
-      REAL(wp) ::  po4mean = 2.165     ! mean value of phosphate
-      REAL(wp) ::  no3mean = 30.90     ! mean value of nitrate
-      REAL(wp) ::  silmean = 91.51     ! mean value of silicate
+      REAL(wp) ::  po4mean = 2.174     ! mean value of phosphate
+      REAL(wp) ::  no3mean = 31.00     ! mean value of nitrate
+      REAL(wp) ::  silmean = 90.33     ! mean value of silicate
       !
       REAL(wp) :: zarea, zalksumn, zpo4sumn, zno3sumn, zsilsumn
@@ -534,5 +548,5 @@
          zwork(:,:,:) =   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe) + trn(:,:,:,jpdfe)   &
             &         +   trn(:,:,:,jpbfe) + trn(:,:,:,jpsfe)                      &
-            &         + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) )  * ferat3    
+            &         + trn(:,:,:,jpzoo) * feratz + trn(:,:,:,jpmes) * feratm
          !
          ferbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )  
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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 14276)
@@ -88,4 +88,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvpuptk    !: Maximum potential uptake rate of picophyto
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvduptk    !: Maximum potential uptake rate of diatoms
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xqfuncfecp !: 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimnpn, xlimnpp, xlimnpd
 
    ! Coefficient for iron limitation following Flynn and Hipkin (1999)
@@ -127,10 +129,10 @@
       REAL(wp) ::   zconc0p, zconc0pnh4, zconc0ppo4, zconcpfe, zconcnfe, zconcdfe
       REAL(wp) ::   fanano, fananop, fananof, fadiat, fadiatp, fadiatf
-      REAL(wp) ::   fapico, fapicop, fapicof
+      REAL(wp) ::   fapico, fapicop, fapicof, zlimpo4, zlimdop
       REAL(wp) ::   zrpho, zrass, zcoef, zfuptk, zratchl
       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
+      REAL(wp) ::   zqfemn, zqfemp, zqfemd, zbactno3, zbactnh4, zbiron
+      REAL(wp) ::   znutlimtot, zlimno3, zlimnh4, zlim1f, zsizetmp
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrassn, zrassp, zrassd
       !!---------------------------------------------------------------------
       !
@@ -143,13 +145,4 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ! 
-               ! Tuning of the iron concentration to a minimum level that
-               ! is set to the detection limit
-               ! --------------------------------------------------------
-               zno3    = trb(ji,jj,jk,jpno3) / 40.e-6
-               zferlim = MAX( 3e-11 * zno3 * zno3, 5e-12 )
-               zferlim = MIN( zferlim, 7e-11 )
-               trb(ji,jj,jk,jpfer) = MAX( trb(ji,jj,jk,jpfer), zferlim )
-
                ! Computation of the Chl/C ratio of each phytoplankton group
                ! ----------------------------------------------------------
@@ -190,10 +183,10 @@
                ! From Talmy et al. (2014) and Maranon et al. (2013)
                ! -------------------------------------------------------
-               xqnnmin(ji,jj,jk) = qnnmin * sizen(ji,jj,jk)**(-0.3)
+               xqnnmin(ji,jj,jk) = qnnmin * sizen(ji,jj,jk)**(-0.36)
                xqnnmax(ji,jj,jk) = qnnmax
-               xqndmin(ji,jj,jk) = qndmin * sized(ji,jj,jk)**(-0.3)
+               xqndmin(ji,jj,jk) = qndmin * sized(ji,jj,jk)**(-0.36)
                xqndmax(ji,jj,jk) = qndmax
-               xqnpmin(ji,jj,jk) = qnpmin * sizep(ji,jj,jk)**(-0.48)
-               xqnpmax(ji,jj,jk) = qnpmax * sizep(ji,jj,jk)**(-0.21)
+               xqnpmin(ji,jj,jk) = qnpmin * sizep(ji,jj,jk)**(-0.36)
+               xqnpmax(ji,jj,jk) = qnpmax
 
                ! Computation of the optimal allocation parameters
@@ -201,5 +194,5 @@
                ! Smith et al.
                ! ---------------------------------------------------
-
+               zbiron = ( 75.0 * ( 1.0 - plig(ji,jj,jk) ) + plig(ji,jj,jk) ) * biron(ji,jj,jk)
                ! Nanophytoplankton
                znutlim = MAX( trb(ji,jj,jk,jpnh4) / zconc0nnh4,    &
@@ -208,5 +201,5 @@
                znutlim = trb(ji,jj,jk,jppo4) / zconc0npo4
                fananop = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
-               znutlim = biron(ji,jj,jk) / zconcnfe
+               znutlim = zbiron / zconcnfe
                fananof = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
 
@@ -217,5 +210,5 @@
                znutlim = trb(ji,jj,jk,jppo4) / zconc0ppo4
                fapicop = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
-               znutlim = biron(ji,jj,jk) / zconcpfe
+               znutlim = zbiron / zconcpfe
                fapicof = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
 
@@ -226,5 +219,5 @@
                znutlim = trb(ji,jj,jk,jppo4) / zconc0dpo4
                fadiatp = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
-               znutlim = biron(ji,jj,jk) / zconcdfe
+               znutlim = zbiron / zconcdfe
                fadiatf = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
 
@@ -232,6 +225,9 @@
                ! heterotrophic bacteria
                ! -------------------------------------------------
-               zbactnh4 = trb(ji,jj,jk,jpnh4) / ( concbnh4 + trb(ji,jj,jk,jpnh4) )
-               zbactno3 = trb(ji,jj,jk,jpno3) / ( concbno3 + trb(ji,jj,jk,jpno3) ) * (1. - zbactnh4)
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( concbno3 + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( concbno3 + trb(ji,jj,jk,jpno3) )
+               znutlimtot = ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) ) / ( concbno3 + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               zbactnh4 = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               zbactno3 = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                zlim1    = zbactno3 + zbactnh4
@@ -250,18 +246,21 @@
                ! Limitation of N based nutrients uptake (NO3 and NH4)
                zfalim = (1.-fanano) / fanano
-               xnanonh4(ji,jj,jk) = (1. - fanano) * trb(ji,jj,jk,jpnh4) / ( zfalim * zconc0nnh4 + trb(ji,jj,jk,jpnh4) )
-               xnanono3(ji,jj,jk) = (1. - fanano) * trb(ji,jj,jk,jpno3) / ( zfalim * zconc0n + trb(ji,jj,jk,jpno3) )  &
-               &                    * (1. - xnanonh4(ji,jj,jk))
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( zconc0n + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( zconc0n + trb(ji,jj,jk,jpno3) )
+               znutlimtot = (1. - fanano) * ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) ) / ( zfalim * zconc0n + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               xnanonh4(ji,jj,jk) = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               xnanono3(ji,jj,jk) = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                ! Limitation of P based nutrients (PO4 and DOP)
                zfalim = (1.-fananop) / fananop
-               xnanopo4(ji,jj,jk) = (1. - fananop) * trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zfalim * zconc0npo4 )
-               xnanodop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + xkdoc )   &
-               &                    * ( 1.0 - xnanopo4(ji,jj,jk) )
-               xnanodop(ji,jj,jk) = 0.
+               zlimpo4 = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zconc0npo4 )
+               zlimdop = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + zconc0npo4 )
+               znutlimtot = (1. - fananop) * ( trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) ) / ( zfalim * zconc0npo4 + trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) )
+               xnanopo4(ji,jj,jk) = znutlimtot * 100.0 * zlimpo4 / ( zlimdop + 100.0 * zlimpo4 + rtrn )
+               xnanodop(ji,jj,jk) = znutlimtot * zlimdop / ( zlimdop + 100.0 * zlimpo4 + rtrn )
                !
                ! Limitation of Fe uptake
                zfalim = (1.-fananof) / fananof
-               xnanofer(ji,jj,jk) = (1. - fananof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcnfe )
+               xnanofer(ji,jj,jk) = (1. - fananof) * zbiron / ( zbiron + zfalim * zconcnfe )
                !
                ! The minimum iron quota depends on the size of PSU, respiration
@@ -270,9 +269,9 @@
                zratiof   = trb(ji,jj,jk,jpnfe) * z1_trnphy
                zqfemn = xcoef1 * znanochl + xcoef2 + xcoef3 * xnanono3(ji,jj,jk)
+               xqfuncfecn(ji,jj,jk) = zqfemn + qfnopt
                !
                zration = trb(ji,jj,jk,jpnph) * z1_trnphy
                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)  &
+               fvnuptk(ji,jj,jk) = 2.5 * zpsiuptk * xqnnmin(ji,jj,jk) / (zration + rtrn)  &
                &                   * MAX(0., (1. - zratchl * znanochl / 12. ) )
                !
@@ -282,5 +281,5 @@
                ! 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) )  &
+               zlim1f = max(0., ( 1.13 - xqnnmin(ji,jj,jk) )  &
                &          / (xqnnmax(ji,jj,jk) - xqnnmin(ji,jj,jk) ) ) * xqnnmax(ji,jj,jk)
                zlim3  = MAX( 0.,( zratiof - zqfemn ) / qfnopt )
@@ -297,28 +296,32 @@
                ! Limitation of N based nutrients uptake (NO3 and NH4) 
                zfalim = (1.-fapico) / fapico 
-               xpiconh4(ji,jj,jk) = (1. - fapico) * trb(ji,jj,jk,jpnh4) / ( zfalim * zconc0pnh4 + trb(ji,jj,jk,jpnh4) )
-               xpicono3(ji,jj,jk) = (1. - fapico) * trb(ji,jj,jk,jpno3) / ( zfalim * zconc0p + trb(ji,jj,jk,jpno3) )  &
-               &                    * (1. - xpiconh4(ji,jj,jk))
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( zconc0p + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( zconc0p + trb(ji,jj,jk,jpno3) )
+               znutlimtot = (1. - fapico) * ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )   &
+               &            / ( zfalim * zconc0p + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               xpiconh4(ji,jj,jk) = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               xpicono3(ji,jj,jk) = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                ! Limitation of P based nutrients uptake (PO4 and DOP)
                zfalim = (1.-fapicop) / fapicop 
-               xpicopo4(ji,jj,jk) = (1. - fapicop) * trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zfalim * zconc0ppo4 )
-               xpicodop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + xkdoc )   &
-               &                    * ( 1.0 - xpicopo4(ji,jj,jk) )
-               xpicodop(ji,jj,jk) = 0.
+               zlimpo4 = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zconc0ppo4 )
+               zlimdop = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + zconc0ppo4 )
+               znutlimtot = (1. - fapicop) * ( trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) ) / ( zfalim * zconc0ppo4 + trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) )
+               xpicopo4(ji,jj,jk) = znutlimtot * 100.0 * zlimpo4 / ( zlimdop + 100.0 * zlimpo4 + rtrn )
+               xpicodop(ji,jj,jk) = znutlimtot * zlimdop / ( zlimdop + 100.0 * zlimpo4 + rtrn )
                !
                zfalim = (1.-fapicof) / fapicof
-               xpicofer(ji,jj,jk) = (1. - fapicof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcpfe )
+               xpicofer(ji,jj,jk) = (1. - fapicof) * zbiron / ( zbiron + zfalim * zconcpfe )
                !
                ! The minimum iron quota depends on the size of PSU, respiration
                ! and the reduction of nitrate following the parameterization 
                ! proposed by Flynn and Hipkin (1999)
-               zratiof   = trb(ji,jj,jk,jppfe) * z1_trnpic
+               zratiof = trb(ji,jj,jk,jppfe) * z1_trnpic
                zqfemp = xcoef1 * zpicochl + xcoef2 + xcoef3 * xpicono3(ji,jj,jk)
+               xqfuncfecp(ji,jj,jk) = zqfemp + qfpopt
                !
                zration   = trb(ji,jj,jk,jpnpi) * z1_trnpic
                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)  &
+               fvpuptk(ji,jj,jk) = 2.5 * zpsiuptk * xqnpmin(ji,jj,jk) / (zration + rtrn)  &
                &                   * MAX(0., (1. - zratchl * zpicochl / 12. ) ) 
                !
@@ -328,5 +331,5 @@
                ! 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) )  &
+               zlim1f   = max(0., (1.29 - xqnpmin(ji,jj,jk) )  &
                &          / (xqnpmax(ji,jj,jk) - xqnpmin(ji,jj,jk) ) ) * xqnpmax(ji,jj,jk)
                zlim3    = MAX( 0.,( zratiof - zqfemp ) / qfpopt )
@@ -344,18 +347,22 @@
                ! Limitation of N based nutrients uptake (NO3 and NH4)
                zfalim = (1.-fadiat) / fadiat 
-               xdiatnh4(ji,jj,jk) = (1. - fadiat) * trb(ji,jj,jk,jpnh4) / ( zfalim * zconc1dnh4 + trb(ji,jj,jk,jpnh4) )
-               xdiatno3(ji,jj,jk) = (1. - fadiat) * trb(ji,jj,jk,jpno3) / ( zfalim * zconc1d + trb(ji,jj,jk,jpno3) )  &
-               &                    * (1. - xdiatnh4(ji,jj,jk))
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( zconc1d + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( zconc1d + trb(ji,jj,jk,jpno3) )
+               znutlimtot = (1.0 - fadiat) * ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )    &
+               &            / ( zfalim * zconc1d + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               xdiatnh4(ji,jj,jk) = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               xdiatno3(ji,jj,jk) = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                ! Limitation of P based nutrients uptake (PO4 and DOP)
                zfalim = (1.-fadiatp) / fadiatp
-               xdiatpo4(ji,jj,jk) = (1. - fadiatp) * trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zfalim * zconc0dpo4 )
-               xdiatdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + xkdoc )  &
-               &                    * ( 1.0 - xdiatpo4(ji,jj,jk) )
-               xdiatdop(ji,jj,jk) = 0.
+               zlimpo4 = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zconc0dpo4 )
+               zlimdop = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + zconc0dpo4 )
+               znutlimtot = (1. - fadiatp) * ( trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) ) / ( zfalim * zconc0dpo4 + trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) )
+               xdiatpo4(ji,jj,jk) = znutlimtot * 100.0 * zlimpo4 / ( zlimdop + 100.0 * zlimpo4 + rtrn )
+               xdiatdop(ji,jj,jk) = znutlimtot * zlimdop / ( zlimdop + 100.0 * zlimpo4 + rtrn )
                !
                ! Limitation of Fe uptake
                zfalim = (1.-fadiatf) / fadiatf
-               xdiatfer(ji,jj,jk) = (1. - fadiatf) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcdfe )
+               xdiatfer(ji,jj,jk) = (1. - fadiatf) * zbiron / ( zbiron + zfalim * zconcdfe )
                !
                ! The minimum iron quota depends on the size of PSU, respiration
@@ -364,9 +371,9 @@
                zratiof   = trb(ji,jj,jk,jpdfe) * z1_trndia
                zqfemd = xcoef1 * zdiatchl + xcoef2 + xcoef3 * xdiatno3(ji,jj,jk)
+               xqfuncfecd(ji,jj,jk) = zqfemd + qfdopt
                !
                zration   = trb(ji,jj,jk,jpndi) * z1_trndia
                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)   &
+               fvduptk(ji,jj,jk) = 2.5 * zpsiuptk * xqndmin(ji,jj,jk) / (zration + rtrn)   &
                &                   * MAX(0., (1. - zratchl * zdiatchl / 12. ) ) 
                !
@@ -376,5 +383,5 @@
                ! 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) )    &
+65             zlim1f   = max(0., (1.13 - xqndmin(ji,jj,jk) )    &
                &          / (xqndmax(ji,jj,jk) - xqndmin(ji,jj,jk) ) )   &
                &          * xqndmax(ji,jj,jk)
@@ -408,12 +415,12 @@
                ! ------------------------------
                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 )
+               ! Computed from Inomura et al. (2020) using Pavlova Lutheri
+               zrpho  = 11.55 * trb(ji,jj,jk,jpnch) / ( trb(ji,jj,jk,jpphy) * 12. + 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.
+               zrassn(ji,jj,jk) = zrass
+               xqpnmin(ji,jj,jk) = ( 0.0 + 0.0078 + 0.62/4. * 0.0783 ) * 16.
+               xqpnmax(ji,jj,jk) = ( zrpho * 0.0089 + zrass * 0.0783 ) * 16.
+               xqpnmax(ji,jj,jk) = xqpnmax(ji,jj,jk) + (0.033 + 0.0078 ) * 16.
                xqpnmax(ji,jj,jk) = MIN( qpnmax, xqpnmax(ji,jj,jk) )
-
 
                ! Size estimation of picophytoplankton based on total biomass
@@ -425,11 +432,12 @@
                ! N/P ratio of picophytoplankton
                ! ------------------------------
-               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 )
+               zfuptk = 0.2 + 0.12 / ( 0.8 * sizep(ji,jj,jk) + rtrn )
+               ! Computed from Inomura et al. (2020) using a synechococcus
+               zrpho = 13.4 * trb(ji,jj,jk,jppch) / ( trb(ji,jj,jk,jppic) * 12. + 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
+               zrassp(ji,jj,jk) = zrass
+               xqppmin(ji,jj,jk) = ( (0.0 + 0.0078 ) + 0.4/4. * 0.0517 ) * 16.
+               xqppmax(ji,jj,jk) = ( zrpho * 0.0076 + zrass * 0.0517 ) * 16.
+               xqppmax(ji,jj,jk) = xqppmax(ji,jj,jk) +  (0.033 + 0.0078 ) * 16
                xqppmax(ji,jj,jk) = MIN( qppmax, xqppmax(ji,jj,jk) )
 
@@ -443,11 +451,11 @@
                ! --------------------
                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 )
+               ! Computed from Inomura et al. (2020) using a synechococcus
+               zrpho = 8.08 * trb(ji,jj,jk,jpdch) / ( trb(ji,jj,jk,jpndi) * 12. + 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.
+               zrassd(ji,jj,jk)=zrass
+               xqpdmin(ji,jj,jk) = ( ( 0.0 + 0.0078 ) + 0.66/4. * 0.0783 ) * 16.
+               xqpdmax(ji,jj,jk) = ( zrpho * 0.0135 + zrass * 0.0783 ) * 16.
+               xqpdmax(ji,jj,jk) = xqpdmax(ji,jj,jk) + ( 0.0078 + 0.033 ) * 16.
                xqpdmax(ji,jj,jk) = MIN(qpdmax, xqpdmax(ji,jj,jk) )
 
@@ -467,11 +475,11 @@
                &        / ( trb(ji,jj,jk,jpnh4) + concnnh4 ) )
                zlim2  = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concnpo4 )
-               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  5.E-11 ) 
-               ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) )
+               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  6.E-11 ) 
+               ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) + 1.8 )
                ztem2  = tsn(ji,jj,jk,jp_tem) - 10.
-               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 )    &
-               &                   * ztem1 / ( 1. + ztem1 ) * MAX( 1., trb(ji,jj,jk,jpphy)*1E6 )   &
+               zetot1 = MAX( 0., etot_ndcy(ji,jj,jk) - 1.) / ( 4. + etot_ndcy(ji,jj,jk) ) * 30. / ( 30. + etot_ndcy(ji,jj,jk) ) 
+
+               xfracal(ji,jj,jk) = caco3r * xlimphy(ji,jj,jk)     &
+               &                   * ztem1 / ( 0.1 + ztem1 ) * MAX( 1., trb(ji,jj,jk,jpphy)*1E6 )   &
                   &                * ( 1. + EXP(-ztem2 * ztem2 / 25. ) )         &
                   &                * zetot1 * MIN( 1., 50. / ( hmld(ji,jj) + rtrn ) )
@@ -508,4 +516,7 @@
         IF( iom_use( "SIZEP"   ) ) CALL iom_put( "SIZEP"  , sizep(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
         IF( iom_use( "SIZED"   ) ) CALL iom_put( "SIZED"  , sized(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "RASSN"   ) ) CALL iom_put( "RASSN"  , zrassn(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "RASSP"   ) ) CALL iom_put( "RASSP"  , zrassp(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "RASSD"   ) ) CALL iom_put( "RASSD"  , zrassd(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
       ENDIF
       !
@@ -640,5 +651,7 @@
          &      fvnuptk (jpi,jpj,jpk), fvduptk (jpi,jpj,jpk),       &
          &      xlimphys(jpi,jpj,jpk), xlimdias(jpi,jpj,jpk),       &
-         &      xlimpics(jpi,jpj,jpk),                              &
+         &      xlimnpp (jpi,jpj,jpk), xlimnpn (jpi,jpj,jpk),       &
+         &      xlimnpd (jpi,jpj,jpk),                              &
+         &      xlimpics(jpi,jpj,jpk), xqfuncfecp(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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 14276)
@@ -88,5 +88,5 @@
       REAL(wp) :: zepsherf, zepshert, zepsherq, zepsherv, zrespirc, zrespirn, zrespirp, zbasresb, zbasresi
       REAL(wp) :: zgraztotc, zgraztotn, zgraztotp, zgraztotf, zbasresn, zbasresp, zbasresf
-      REAL(wp) :: zgratmp, zgradoct, zgradont, zgrareft, zgradopt
+      REAL(wp) :: zgradoct, zgradont, zgrareft, zgradopt
       REAL(wp) :: zprcaca, zmortz, zexcess
       REAL(wp) :: zbeta, zrespz, ztortz, zgrasratp, zgrasratn, zgrasratf
@@ -104,5 +104,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrarep, zgraren, zgrapon, zgrapop
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgradoc, zgradon, zgradop
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigrem, zgramigref, zgramigpoc, zgramigpof, zstrn
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigrem, zgramigref, zgramigpoc, zgramigpof
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigrep, zgramigren, zgramigpop, zgramigpon
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigdoc, zgramigdop, zgramigdon
@@ -190,5 +190,5 @@
                ! to avoid starvation.
                ! ----------------------------------------------------------
-               zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
+               zsigma = 1.0 - zdenom**3/(0.1**3+zdenom**3)
                zsigma = xsigma2 + xsigma2del * zsigma
                ! Nanophytoplankton and diatoms are the only preys considered
@@ -267,5 +267,5 @@
 
                zgraztotc  = zgrazdc + zgrazz + zgraznc + zgrazm + zgrazpoc + zgrazffep + zgrazffeg
-               zgraztotf  = zgrazdf + zgraznf + ( zgrazz + zgrazm ) * ferat3 + zgrazpof &
+               zgraztotf  = zgrazdf + zgraznf + zgrazz * feratz + zgrazm * feratm + zgrazpof &
                &            + zgrazfffp + zgrazfffg
                zgraztotn  = zgrazdn + (zgrazm + zgrazz) * no3rat3 + zgraznn + zgrazpon  &
@@ -291,5 +291,5 @@
                ! Fulton, 2012)
                ! -----------------------------------------------------------------------------------
-               zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / ferat3)
+               zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / feratm)
                zbeta     = MAX(0., (epsher2 - epsher2min) )
                zepsherf  = epsher2min + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
@@ -312,5 +312,5 @@
                zexcess  = ( zgrasratp/ po4rat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
                zbasresp = zbasresi * zexcess * zgrasratp
-               zexcess  = ( zgrasratf/ ferat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
+               zexcess  = ( zgrasratf/ feratm - zepshert ) / ( 1.0 - zepshert + rtrn)
                zbasresf = zbasresi * zexcess * zgrasratf
 
@@ -320,10 +320,9 @@
                zgradont = (1. - unass2n) * zgraztotn - zepsherv * no3rat3 * zgraztotc - zbasresn
                zgradopt = (1. - unass2p) * zgraztotp - zepsherv * po4rat3 * zgraztotc - zbasresp
-               zgrareft = (1. - unass2c) * zgraztotf - zepsherv * ferat3 * zgraztotc - zbasresf
-               ztmp1   = ( 1. - epsher2 - unass2c ) /( 1. - 0.8 * epsher2 ) * ztortz
+               zgrareft = (1. - unass2c) * zgraztotf - zepsherv * feratm * zgraztotc - zbasresf
+               ztmp1    = (1. - epsher2 - unass2c) /( 1. - epsher2 ) * ztortz
                zgradoc(ji,jj,jk) = (zgradoct + ztmp1) * ssigma2
                zgradon(ji,jj,jk) = (zgradont + no3rat3 * ztmp1) * ssigma2
                zgradop(ji,jj,jk) = (zgradopt + po4rat3 * ztmp1) * ssigma2
-               zgratmp = 0.2 * epsher2 /( 1. - 0.8 * epsher2 ) * ztortz
 
                !  Since only semilabile DOM is represented in PISCES
@@ -331,15 +330,15 @@
                !  as dissolved inorganic compounds (ssigma2)
                !  --------------------------------------------------
-               zgrarem(ji,jj,jk) = zgratmp + ( zgradoct + ztmp1 ) * (1.0 - ssigma2)
-               zgraren(ji,jj,jk) = no3rat3 * zgratmp + ( zgradont + no3rat3 * ztmp1 ) * (1.0 - ssigma2)
-               zgrarep(ji,jj,jk) = po4rat3 * zgratmp + ( zgradopt + po4rat3 * ztmp1 ) * (1.0 - ssigma2)
-               zgraref(ji,jj,jk) = zgrareft + ferat3 * ( ztmp1 + zgratmp )
+               zgrarem(ji,jj,jk) = ( zgradoct + ztmp1 ) * (1.0 - ssigma2)
+               zgraren(ji,jj,jk) = ( zgradont + no3rat3 * ztmp1 ) * (1.0 - ssigma2)
+               zgrarep(ji,jj,jk) = ( zgradopt + po4rat3 * ztmp1 ) * (1.0 - ssigma2)
+               zgraref(ji,jj,jk) = zgrareft + feratm * ztmp1
 
                !   Defecation as a result of non assimilated products
                !   --------------------------------------------------
-               zgrapoc(ji,jj,jk)  = zgraztotc * unass2c + unass2c / ( 1. - 0.8 * epsher2 ) * ztortz
-               zgrapon(ji,jj,jk)  = zgraztotn * unass2n + no3rat3 * unass2n / ( 1. - 0.8 * epsher2 ) * ztortz
-               zgrapop(ji,jj,jk)  = zgraztotp * unass2p + po4rat3 * unass2p / ( 1. - 0.8 * epsher2 ) * ztortz
-               zgrapof(ji,jj,jk)  = zgraztotf * unass2c + ferat3  * unass2c / ( 1. - 0.8 * epsher2 ) * ztortz
+               zgrapoc(ji,jj,jk)  = zgraztotc * unass2c + unass2c / ( 1. - epsher2 ) * ztortz
+               zgrapon(ji,jj,jk)  = zgraztotn * unass2n + no3rat3 * unass2n / ( 1. - epsher2 ) * ztortz
+               zgrapop(ji,jj,jk)  = zgraztotp * unass2p + po4rat3 * unass2p / ( 1. - epsher2 ) * ztortz
+               zgrapof(ji,jj,jk)  = zgraztotf * unass2c + feratm  * unass2c / ( 1. - epsher2 ) * ztortz
 
                !  Addition of respiration to the release of inorganic nutrients
@@ -348,5 +347,5 @@
                zgraren(ji,jj,jk) = zgraren(ji,jj,jk) + zbasresn + zrespirc * no3rat3
                zgrarep(ji,jj,jk) = zgrarep(ji,jj,jk) + zbasresp + zrespirc * po4rat3
-               zgraref(ji,jj,jk) = zgraref(ji,jj,jk) + zbasresf + zrespirc * ferat3
+               zgraref(ji,jj,jk) = zgraref(ji,jj,jk) + zbasresf + zrespirc * feratm
 
                !   Update the arrays TRA which contain the biological sources and
@@ -404,5 +403,4 @@
          ALLOCATE( zgramigrep(jpi,jpj), zgramigren(jpi,jpj), zgramigpop(jpi,jpj), zgramigpon(jpi,jpj) )
          ALLOCATE( zgramigdoc(jpi,jpj), zgramigdon(jpi,jpj), zgramigdop(jpi,jpj) )
-         ALLOCATE( zstrn(jpi,jpj) )
          zgramigrem(:,:)  = 0.0   ;   zgramigref(:,:) = 0.0
          zgramigrep(:,:)  = 0.0   ;   zgramigren(:,:) = 0.0
@@ -411,20 +409,4 @@
          zgramigdoc(:,:)  = 0.0   ;   zgramigdon(:,:) = 0.0
          zgramigdop(:,:)  = 0.0   
-
-         ! compute the day length depending on latitude and the day
-         zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
-         zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
-
-         ! day length in hours
-         zstrn(:,:) = 0.
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
-               zargu = MAX( -1., MIN(  1., zargu ) )
-               zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
-               zstrn(ji,jj) = MIN(0.75, MAX( 0.25, zstrn(ji,jj) / 24.) )
-            END DO
-         END DO
-
 
         ! Compute the amount of materials that will go into vertical migration
@@ -435,5 +417,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zmigreltime = (1. - zstrn(ji,jj))
+                  zmigreltime = (1. - strn(ji,jj))
                   IF ( gdept_n(ji,jj,jk) <= heup(ji,jj) ) THEN
                      zgramigrem(ji,jj) = zgramigrem(ji,jj) + xfracmig * zgrarem(ji,jj,jk) * (1. - zmigreltime )    &
@@ -441,5 +423,5 @@
                      zgramigrep(ji,jj) = zgramigrep(ji,jj) + xfracmig * zgrarep(ji,jj,jk) * (1. - zmigreltime )    &
                      &                   * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
-                     zgramigrep(ji,jj) = zgramigren(ji,jj) + xfracmig * zgrarep(ji,jj,jk) * (1. - zmigreltime )    &
+                     zgramigren(ji,jj) = zgramigren(ji,jj) + xfracmig * zgraren(ji,jj,jk) * (1. - zmigreltime )    &
                      &                   * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
                      zgramigref(ji,jj) = zgramigref(ji,jj) + xfracmig * zgraref(ji,jj,jk) * (1. - zmigreltime )   &
@@ -503,5 +485,4 @@
          DEALLOCATE( zgramigrep, zgramigren, zgramigpop, zgramigpon )
          DEALLOCATE( zgramigdoc, zgramigdon, zgramigdop )
-         DEALLOCATE( zstrn )
       ! End of the ln_dvm_meso part
       ENDIF
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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmicro.F90	(revision 14276)
@@ -110,6 +110,7 @@
                ! 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)
+               zproport  = min(sized(ji,jj,jk),1.8)**(-0.48)*min(1.0, exp(-1.1 * MAX(0., ( sized(ji,jj,jk) - 1.8 ))**0.8 ))
+               zproport2 = sizen(ji,jj,jk)**(-0.48)
+               zproport2 = 1.0
                !  linear mortality of mesozooplankton
                !  A michaelis menten modulation term is used to avoid extinction of 
@@ -162,9 +163,9 @@
                ! to avoid starvation.
                ! ----------------------------------------------------------
-               zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
+               zsigma = 1.0 - zdenom**3/(0.1**3+zdenom**3)
                zsigma = xsigma + xsigmadel * zsigma
-               zdiffpn = exp( -ABS(log(0.5 * sizep(ji,jj,jk) / (3.0 * sizen(ji,jj,jk) + rtrn )) )**2 / zsigma**2 )
+               zdiffpn = exp( -ABS(log(0.7 * 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)
+               zdiffdp = exp( -ABS(log(0.7 * 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 )
@@ -211,5 +212,5 @@
                zgraztotn = zgraznn + zgrazpn + zgrazpon + zgrazdn + zgrazz * no3rat3
                zgraztotp = zgraznp + zgrazpp + zgrazpop + zgrazdp + zgrazz * po4rat3
-               zgraztotf = zgraznf + zgrazpf + zgrazpof + zgrazdf + zgrazz * ferat3
+               zgraztotf = zgraznf + zgrazpf + zgrazpof + zgrazdf + zgrazz * feratz
                !
                ! Grazing by microzooplankton
@@ -230,5 +231,5 @@
                ! Fulton, 2012)
                ! -----------------------------------------------------------------------------------
-               zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / ferat3)
+               zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / feratz)
                zbeta     = MAX( 0., (epsher - epshermin) )
                ! Food density deprivation of GGE
@@ -255,5 +256,5 @@
                zexcess  = ( zgrasratp/ po4rat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
                zbasresp = zbasresi * zexcess * zgrasratp
-               zexcess  = ( zgrasratf/ ferat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
+               zexcess  = ( zgrasratf/ feratz - zepshert ) / ( 1.0 - zepshert + rtrn)
                zbasresf = zbasresi * zexcess * zgrasratf
 
@@ -263,5 +264,5 @@
                zgradont   = (1. - unassn) * zgraztotn - zepsherv * no3rat3 * zgraztotc - zbasresn
                zgradopt   = (1. - unassp) * zgraztotp - zepsherv * po4rat3 * zgraztotc - zbasresp
-               zgrareft   = (1. - unassc) * zgraztotf - zepsherv * ferat3 * zgraztotc - zbasresf
+               zgrareft   = (1. - unassc) * zgraztotf - zepsherv * feratz * zgraztotc - zbasresf
 
                ! Since only semilabile DOM is represented in PISCES
@@ -289,5 +290,5 @@
                zgraren = zgraren + zbasresn + zrespirc * no3rat3
                zgrarep = zgrarep + zbasresp + zrespirc * po4rat3
-               zgraref = zgraref + zbasresf + zrespirc * ferat3
+               zgraref = zgraref + zbasresf + zrespirc * feratz
 
                ! Update of the TRA arrays
@@ -330,5 +331,5 @@
                tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + no3rat3 * ztortz + zgrapon - zgrazpon
                tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + po4rat3 * ztortz + zgrapop - zgrazpop
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * ztortz  + zgrapof - zgrazpof
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + feratz * ztortz  + zgrapof - zgrazpof
                !
                ! Calcite production
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmort.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmort.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmort.F90	(revision 14276)
@@ -28,5 +28,4 @@
    REAL(wp), PUBLIC :: wchlp   !: Quadratic mortality rate of picophytoplankton
    REAL(wp), PUBLIC :: wchld   !: Quadratic mortality rate of diatoms
-   REAL(wp), PUBLIC :: wchldm  !: Maximum quadratic mortality rate of diatoms
    REAL(wp), PUBLIC :: mpratn  !: Linear mortality rate of nanophytoplankton
    REAL(wp), PUBLIC :: mpratp  !: Linear mortality rate of picophytoplankton
@@ -69,5 +68,5 @@
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jk
-      REAL(wp) :: zcompaph
+      REAL(wp) :: zcompaph, zlim1, zlim2
       REAL(wp) :: zfactfe, zfactch, zfactn, zfactp, zprcaca
       REAL(wp) :: ztortp , zrespp , zmortp
@@ -86,5 +85,7 @@
                ! blooms (Doney et al. 1996)
                ! -----------------------------------------------------
-               zrespp = wchln * 1.e6 * xstep * xdiss(ji,jj,jk) * zcompaph * trb(ji,jj,jk,jpphy)
+               zlim2   = xlimphy(ji,jj,jk) * xlimphy(ji,jj,jk)
+               zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) * trb(ji,jj,jk,jpphy)
+               zrespp = wchln * 1.e6 * xstep * zlim1 * xdiss(ji,jj,jk) * zcompaph
 
                ! Phytoplankton linear mortality
@@ -142,5 +143,5 @@
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jk
-      REAL(wp) :: zcompaph
+      REAL(wp) :: zcompaph, zlim1, zlim2
       REAL(wp) :: zfactfe, zfactch, zfactn, zfactp
       REAL(wp) :: ztortp , zrespp , zmortp 
@@ -158,5 +159,7 @@
                ! blooms (Doney et al. 1996)
                ! -----------------------------------------------------
-               zrespp = wchlp * 1.e6 * xstep * xdiss(ji,jj,jk) * zcompaph * trb(ji,jj,jk,jppic)
+               zlim2   = xlimpic(ji,jj,jk) * xlimpic(ji,jj,jk)
+               zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) * trb(ji,jj,jk,jppic)
+               zrespp = wchlp * 1.e6 * xstep * zlim1 * xdiss(ji,jj,jk) * zcompaph
 
                ! Phytoplankton linear mortality
@@ -229,5 +232,5 @@
                zlim2   = xlimdia(ji,jj,jk) * xlimdia(ji,jj,jk)
                zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) 
-               zrespp2 = 1.e6 * xstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
+               zrespp2 = 1.e6 * xstep * wchld * zlim1 * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
 
                ! Phytoplankton linear mortality
@@ -291,5 +294,5 @@
       INTEGER :: ios   ! Local integer output status for namelist read
       !!
-      NAMELIST/namp5zmort/ wchln, wchlp, wchld, wchldm, mpratn, mpratp, mpratd
+      NAMELIST/namp5zmort/ wchln, wchlp, wchld, mpratn, mpratp, mpratd
       !!----------------------------------------------------------------------
 
@@ -310,5 +313,4 @@
          WRITE(numout,*) '    quadratic mortality of picophyto.         wchlp     =', wchlp
          WRITE(numout,*) '    quadratic mortality of diatoms            wchld     =', wchld
-         WRITE(numout,*) '    Additional quadratic mortality of diatoms wchldm    =', wchldm
          WRITE(numout,*) '    nanophyto. mortality rate                 mpratn    =', mpratn
          WRITE(numout,*) '    picophyto. mortality rate                 mpratp    =', mpratp
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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zprod.F90	(revision 14276)
@@ -72,15 +72,16 @@
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zsilfac, znanotot, zpicotot, zdiattot, zconctemp, zconctemp2
-      REAL(wp) ::   zration, zratiop, zratiof, zmax, zsilim, ztn, zadap
-      REAL(wp) ::   zpronmax, zpropmax, zprofmax, zrat
+      REAL(wp) ::   zration, zratiop, zratiof, zmax, ztn, zadap
+      REAL(wp) ::   zpronmax, zpropmax, zprofmax, zratio
       REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zprontot, zproptot, zprodtot
       REAL(wp) ::   zprnutmax, zdocprod, zprochln, zprochld, zprochlp
       REAL(wp) ::   zpislopen, zpislopep, zpisloped
       REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
-      REAL(wp) ::   zfact, zrfact2, zmaxsi, zratiosi, zsizetmp, zlimfac
+      REAL(wp) ::   zqfpmax, zqfnmax, zqfdmax
+      REAL(wp) ::   zfact, zrfact2, zmaxsi, zratiosi, zsizetmp, zlimfac, zsilim
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zmixnano, zmixpico, zmixdiat, zstrn
+      REAL(wp), DIMENSION(jpi,jpj    ) :: zmixnano, zmixpico, zmixdiat
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadp, zpislopeadd
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprnut, zprnutp, zprmaxp, zprmaxn, zprmaxd
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprnut, zprmaxp, zprmaxn, zprmaxd
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprbio, zprpic, zprdia, zysopt
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprchln, zprchlp, zprchld
@@ -111,27 +112,13 @@
       zysopt  (:,:,:) = 0._wp
       zrespn  (:,:,:) = 0._wp ; zrespp  (:,:,:) = 0._wp ; zrespd  (:,:,:) = 0._wp 
+      zmxl_fac(:,:,:) = 0._wp ; zmxl_chl(:,:,:) = 0._wp
+      consfe3 (:,:,:) = 0._wp
 
       ! Computation of the optimal production rates and nutrient uptake
       ! rates. Based on a Q10 description of the thermal dependency.
-      zprnut (:,:,:) = 0.6_wp * (1.0 + zpsino3 * qnnmax ) * r1_rday * tgfunc(:,:,:)
-      zprnutp(:,:,:) =  0.6_wp * (1. + zpsino3 * qnpmax ) * r1_rday * tgfunc3(:,:,:)
-      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
-      ! Astronomical parameterization taken from HAMOCC3
-      zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
-      zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
-
-      ! day length in hours
-      zstrn(:,:) = 0.
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
-            zargu = MAX( -1., MIN(  1., zargu ) )
-            zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
-         END DO
-      END DO
+      zprnut (:,:,:) =  0.8_wp * r1_rday * tgfunc(:,:,:)
+      zprmaxn(:,:,:) =  0.8_wp * (1. + zpsino3 * qnnmax ) * r1_rday * tgfunc(:,:,:)
+      zprmaxd(:,:,:) =  0.8_wp * (1. + zpsino3 * qndmax ) * r1_rday * tgfunc(:,:,:)
+      zprmaxp(:,:,:) =  0.6_wp * (1. + zpsino3 * qnpmax ) * r1_rday * tgfunc(:,:,:)
 
       ! Impact of the day duration and light intermittency on phytoplankton growth
@@ -146,5 +133,5 @@
             DO ji = 1, jpi
                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                  zval = MAX( 1., zstrn(ji,jj) )
+                  zval = MAX( 1., strn(ji,jj) )
                   IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
                      zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
@@ -161,7 +148,6 @@
       zprpic(:,:,:) = zprmaxp(:,:,:) * zmxl_fac(:,:,:)
 
-
       ! Maximum light intensity
-      zdaylen(:,:) = MAX(1., zstrn(:,:)) / 24.
+      zdaylen(:,:) = MAX(1., strn(:,:)) / 24.
 
       ! Computation of the P-I slope for nanos, picos and diatoms
@@ -192,7 +178,7 @@
                   ! Actual light levels are used here 
                   !  ---------------------------------------------
-                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) / zmxl_chl(ji,jj,jk) )  )
-                  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))  )
+                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) / zmxl_fac(ji,jj,jk) )  )
+                  zprpic(ji,jj,jk) = zprpic(ji,jj,jk) * ( 1.- EXP( -zpislopep * epico(ji,jj,jk) / zmxl_fac(ji,jj,jk) )  )
+                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) / zmxl_fac(ji,jj,jk) )  )
 
                   !  Computation of production function for Chlorophyll
@@ -227,5 +213,5 @@
                   ! -----------------------------------------------------------------------
                   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) )
+                  zsilim = xlimdia(ji,jj,jk) * zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn )
                   zsiborn = trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil)
                   IF (gphit(ji,jj) < -30 ) THEN
@@ -240,5 +226,5 @@
                      zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zmaxsi
                   ELSE
-                     zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zsilim**0.75 * zmaxsi
+                     zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zsilim**0.7 * zmaxsi
                   ENDIF
                ENDIF
@@ -257,5 +243,4 @@
                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
@@ -284,5 +269,5 @@
                   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 ) )
+                  sizena(ji,jj,jk) = MIN(xsizern, MAX( sizena(ji,jj,jk), zsizetmp ) )
                   ! Maximum potential uptake rate
                   zration = trb(ji,jj,jk,jpnph) / ( trb(ji,jj,jk,jpphy) + rtrn )
@@ -291,20 +276,22 @@
                   zprnutmax = zprnut(ji,jj,jk) * fvnuptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jpphy) * rfact2
                   ! Uptake of nitrogen
-                  zrat = 1.0 - MIN( 1., zration / (xqnnmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
+                  zratio = 1.0 - MIN( 1., zration / (xqnnmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**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) ) )
+                  zpronmax = zpronmax * xqnnmin(ji,jj,jk) / qnnmin
                   zpronewn(ji,jj,jk) = zpronmax * xnanono3(ji,jj,jk)
                   zproregn(ji,jj,jk) = zpronmax * xnanonh4(ji,jj,jk)
                   ! Uptake of phosphorus and DOP
-                  zrat = 1.0 - MIN( 1., zratiop / (xqpnmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
-                  zpropmax = zprnutmax * zmax * xlimnfe(ji,jj,jk) * 16. / 10.
+                  zratio = 1.0 - MIN( 1., zratiop / (xqpnmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**2) ) )
+                  zpropmax = zprnutmax * zmax * xlimnfe(ji,jj,jk)
                   zpropo4n(ji,jj,jk) = zpropmax * xnanopo4(ji,jj,jk)
                   zprodopn(ji,jj,jk) = zpropmax * xnanodop(ji,jj,jk)
                   ! Uptake of iron
-                  zrat = 1.0 - MIN( 1., zratiof / qfnmax )
-                  zmax = MAX(0., MIN(1., zrat**2/ (0.05**2 + zrat**2) ) )
-                  zprofmax = zprnutmax * qfnmax * zmax 
+                  zqfnmax = xqfuncfecn(ji,jj,jk) + ( qfnmax - xqfuncfecn(ji,jj,jk) ) * xlimnpn(ji,jj,jk)
+                  zratio = 1.0 - MIN( 1., zratiof / zqfnmax )
+                  zmax = MAX(0., MIN(1., zratio**2/ (0.05**2 + zratio**2) ) )
+                  zprofmax = zprnutmax * zqfnmax * zmax 
                   zprofen(ji,jj,jk) = zprofmax * xnanofer(ji,jj,jk)    &
                   &          * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn  &
@@ -341,22 +328,24 @@
                   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 = zprnutp(ji,jj,jk) * fvpuptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jppic) * rfact2
+                  zprnutmax = zprnut(ji,jj,jk) * fvpuptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jppic) * rfact2
                   ! Uptake of nitrogen
-                  zrat = 1.0 - MIN( 1., zration / (xqnpmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., zrat**2/ (0.05**2 + zrat**2) ) )
+                  zratio = 1.0 - MIN( 1., zration / (xqnpmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2/ (0.05**2 + zratio**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) ) )
+                  zpronmax = zpronmax * xqnpmin(ji,jj,jk) / qnnmin
                   zpronewp(ji,jj,jk) = zpronmax * xpicono3(ji,jj,jk) 
                   zproregp(ji,jj,jk) = zpronmax * xpiconh4(ji,jj,jk)
                   ! Uptake of phosphorus
-                  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.
+                  zratio = 1.0 - MIN( 1., zratiop / (xqppmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**2) ) )
+                  zpropmax = zprnutmax * zmax * xlimpfe(ji,jj,jk) 
                   zpropo4p(ji,jj,jk) = zpropmax * xpicopo4(ji,jj,jk)
                   zprodopp(ji,jj,jk) = zpropmax * xpicodop(ji,jj,jk)
                   ! Uptake of iron
-                  zrat = 1.0 - MIN( 1., zratiof / qfpmax )
-                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
-                  zprofmax = zprnutmax * qfpmax * zmax
+                  zqfpmax = xqfuncfecp(ji,jj,jk) + ( qfpmax - xqfuncfecp(ji,jj,jk) ) * xlimnpp(ji,jj,jk)
+                  zratio = 1.0 - MIN( 1., zratiof / zqfpmax )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**2) ) )
+                  zprofmax = zprnutmax * zqfpmax * zmax
                   zprofep(ji,jj,jk) = zprofmax * xpicofer(ji,jj,jk)  &
                   &          * (1. + 0.8 * xpicono3(ji,jj,jk) / ( rtrn   &
@@ -395,20 +384,22 @@
                   zprnutmax = zprnut(ji,jj,jk) * fvduptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jpdia) * rfact2
                   ! Uptake of nitrogen
-                  zrat = 1.0 - MIN( 1., zration / (xqndmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
+                  zratio = 1.0 - MIN( 1., zration / (xqndmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**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) ) )
+                  zpronmax = zpronmax * xqndmin(ji,jj,jk) / qnnmin
                   zpronewd(ji,jj,jk) = zpronmax * xdiatno3(ji,jj,jk)
                   zproregd(ji,jj,jk) = zpronmax * xdiatnh4(ji,jj,jk)
                   ! Uptake of phosphorus
-                  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.
+                  zratio = 1.0 - MIN( 1., zratiop / (xqpdmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2/ (0.05**2 + zratio**2) ) )
+                  zpropmax = zprnutmax * zmax * xlimdfe(ji,jj,jk)
                   zpropo4d(ji,jj,jk) = zpropmax * xdiatpo4(ji,jj,jk)
                   zprodopd(ji,jj,jk) = zpropmax * xdiatdop(ji,jj,jk)
                   ! Uptake of iron
-                  zrat = 1.0 - MIN( 1., zratiof / qfdmax )
-                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
-                  zprofmax = zprnutmax * qfdmax * zmax
+                  zqfdmax = xqfuncfecd(ji,jj,jk) + ( qfdmax - xqfuncfecd(ji,jj,jk) ) * xlimnpd(ji,jj,jk)
+                  zratio = 1.0 - MIN( 1., zratiof / zqfdmax )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**2) ) )
+                  zprofmax = zprnutmax * zqfdmax * zmax
                   zprofed(ji,jj,jk) = zprofmax * xdiatfer(ji,jj,jk)    &
                   &          * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn   &
@@ -503,4 +494,6 @@
               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
+              consfe3(ji,jj,jk)   = zfeup * 75.0 / ( rtrn + ( plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) )   &
+              &                     * trb(ji,jj,jk,jpfer) ) / rfact2
               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)  &
@@ -529,8 +522,8 @@
                  zfeup    = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) + texcretp * zprofep(ji,jj,jk)
                  tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp    &
-                 &       - zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
+                 &       - zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * lthet
                  zpligprod1(ji,jj,jk) = zdocprod * ldocp
                  zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) &
-                 &                      + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
+                 &                      + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * lthet
               END DO
            END DO
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/par_sed.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/par_sed.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/par_sed.F90	(revision 14276)
@@ -61,6 +61,6 @@
    INTEGER, PARAMETER ::  &
       jptrased   = jpsol + jpwat , &
-      jpdia3dsed = 2             , &
-      jpdia2dsed = 12
+      jpdia3dsed = 3             , &
+      jpdia2dsed = 20 
 
 END MODULE par_sed
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedarr.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedarr.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedarr.F90	(revision 14276)
@@ -50,4 +50,5 @@
          jid        = MOD( tab_ind(jn) - 1, jpi ) + 1
          jjd        = ( tab_ind(jn) - 1 ) / jpi + 1
+!         IF ( mig(jid) == 112 .and. mjg(jjd) == 25) write(0,*) 'plante ',jn,ndim1d
          tab1d(jn)  = tab2d(jid, jjd)
       END DO 
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/seddta.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/seddta.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/seddta.F90	(revision 14276)
@@ -157,7 +157,7 @@
 !        zf1    = ( 0.02 * (reac_poc - reac_poc * 0.) + zkapp - reac_poc ) / ( reac_poc / 100. - reac_poc )
 !        zf1    = MIN(0.98, MAX(0., zf1 ) )
-         zf1    = 0.48
-         zf0    = 1.0 - 0.02 - zf1
-         zf2    = 0.02
+         zf1    = 0.7
+         zf2    = 0.04
+         zf0    = 1.0 - zf2 - zf1
          rainrm_dta(ji,jspoc) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf0
          rainrm_dta(ji,jspos) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf1
@@ -177,6 +177,6 @@
       rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * conv2 / mol_wgt(jsclay)   &
       &                            + wacc(1:jpoce) * por1(2) * denssol / mol_wgt(jsclay) / ( rsecday * 365.0 )
+      rainrm_dta(1:jpoce,jsfeo)  = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 * 0.5
       rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * 0.965
-      rainrm_dta(1:jpoce,jsfeo)  = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 / 0.965
 !    rainrm_dta(1:jpoce,jsclay) = 1.0E-4 * conv2 / mol_wgt(jsclay)
 
@@ -214,8 +214,8 @@
 
       IF( lk_iomput ) THEN
-          IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", dust(:,:) * conv2 * 1E4 )
-          IF( iom_use("sflxcal" ) )  CALL iom_put( "sflxcal", trc_data(:,:,13) )
-          IF( iom_use("sflxbsi" ) )  CALL iom_put( "sflxbsi", trc_data(:,:,10) )
-          IF( iom_use("sflxpoc" ) )  CALL iom_put( "sflxpoc", trc_data(:,:,11) + trc_data(:,:,12) )
+          IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", dust(:,:) * 1E3 / 1.E4 )
+          IF( iom_use("sflxcal" ) )  CALL iom_put( "sflxcal", trc_data(:,:,14) / 1.E4 )
+          IF( iom_use("sflxbsi" ) )  CALL iom_put( "sflxbsi", trc_data(:,:,11) / 1.E4 )
+          IF( iom_use("sflxpoc" ) )  CALL iom_put( "sflxpoc", ( trc_data(:,:,12) + trc_data(:,:,13) ) / 1.E4 )
       ENDIF
 
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedini.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedini.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedini.F90	(revision 14276)
@@ -371,14 +371,4 @@
       denssol = 2.6
 
-      ! Initialization of diffusion coefficient as function of porosity [cm**2/s]
-      !--------------------------------------------------------------------
-!      DO jn = 1, jpsol
-!         DO jk = 1, jpksed
-!            DO ji = 1, jpoce
-!               diff(ji,jk,jn) = dcoef / ( 1.0 - 2.0 * log(por(jk)) )
-!            END DO
-!         END DO
-!      END DO
-
       ! Accumulation rate from Burwicz et al. (2011). This is used to
       ! compute the flux of clays and minerals
@@ -387,5 +377,6 @@
           ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 )
           ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 )
-          wacc(ji) = ztmp1 + ztmp2
+          wacc(ji) = (ztmp1 + ztmp2)/10.0
+          wacc(ji) = ztmp2 / 10.0
       END DO
 
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedwri.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedwri.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedwri.F90	(revision 14276)
@@ -53,5 +53,5 @@
       IF (lwp) WRITE(numsed,*) ' '
       
-      ALLOCATE( zdta(jpoce,jpksed) )    ;   ALLOCATE( zflx(jpoce,jpwatp1) )
+      ALLOCATE( zdta(jpoce,jpksed) )    ;   ALLOCATE( zflx(jpoce,jptrased+1) )
 
       ! Initialize variables
@@ -87,4 +87,7 @@
       CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,2)  , iarroce(1:jpoce), &
          &                   co3por(1:jpoce,1:jpksed)  )
+
+      CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,3)  , iarroce(1:jpoce), &
+         &                   sedligand(1:jpoce,1:jpksed)  )
       
 !      flxsedi3d = 0.
@@ -94,5 +97,13 @@
          DO ji = 1, jpoce
             zflx(ji,jw) = ( pwcp(ji,1,jw) - pwcp_dta(ji,jw) ) &
-               &         * 1.e3 / 1.e2 * dzkbot(ji) / r2dttrc
+               &         * 1.e3 / ( 1.e2 * dzkbot(ji) ) / 1.E4 / r2dttrc
+         ENDDO
+      ENDDO
+
+      ! Calculation of fluxes g/cm2/s
+      DO js = 1, jpsol
+         zrate =  1.0 / ( por1(jpksed) ) / r2dttrc
+         DO ji = 1, jpoce
+            zflx(ji,jpwat+js) = zflx(ji,jpwat+js) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate
          ENDDO
       ENDDO
@@ -100,7 +111,7 @@
       ! Calculation of accumulation rate per dt
       DO js = 1, jpsol
-         zrate =  1.0 / ( denssol * por1(jpksed) ) / r2dttrc
+         zrate =  1.0 / ( por1(jpksed) ) / r2dttrc
          DO ji = 1, jpoce
-            zflx(ji,jpwatp1) = zflx(ji,jpwatp1) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate
+            zflx(ji,jptrased+1) = zflx(ji,jptrased+1) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate
          ENDDO
       ENDDO
@@ -109,23 +120,24 @@
          CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jn), iarroce(1:jpoce), zflx(1:jpoce,jn)  )
       END DO
+
       zflx(:,1) = dzdep(:) / dtsed
       CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed), iarroce(1:jpoce), zflx(1:jpoce,1) )
 
-       ! Start writing data
-       ! ---------------------
-       DO jn = 1, jptrased
-          cltra = sedtrcd(jn) ! short title for 3D diagnostic
-          CALL iom_put( cltra, trcsedi(:,:,:,jn) )
-       END DO
+      ! Start writing data
+      ! ---------------------
+      DO jn = 1, jptrased
+         cltra = sedtrcd(jn) ! short title for 3D diagnostic
+         CALL iom_put( cltra, trcsedi(:,:,:,jn) )
+      END DO
 
-       DO jn = 1, jpdia3dsed
-          cltra = seddia3d(jn) ! short title for 3D diagnostic
-          CALL iom_put( cltra, flxsedi3d(:,:,:,jn) )
-       END DO
+      DO jn = 1, jpdia3dsed
+         cltra = seddia3d(jn) ! short title for 3D diagnostic
+         CALL iom_put( cltra, flxsedi3d(:,:,:,jn) )
+      END DO
 
-       DO jn = 1, jpdia2dsed
-          cltra = seddia2d(jn) ! short title for 2D diagnostic
-          CALL iom_put( cltra, flxsedi2d(:,:,jn) )
-       END DO
+      DO jn = 1, jpdia2dsed
+         cltra = seddia2d(jn) ! short title for 2D diagnostic
+         CALL iom_put( cltra, flxsedi2d(:,:,jn) )
+      END DO
 
 
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 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/sms_pisces.F90	(revision 14276)
@@ -43,5 +43,6 @@
    REAL(wp) ::   wsbio2scale      !: Length scale for the variations of wsbio2
    REAL(wp) ::   xkmort           !: Mortality half-saturation constant
-   REAL(wp) ::   ferat3           !: Fe/C in heterotrophic organisms
+   REAL(wp) ::   feratz           !: Fe/C in microzooplankton
+   REAL(wp) ::   feratm           !: Fe/C in mesozooplankton
    REAL(wp) ::   ldocp            !: Ligand production ratio during PP
    REAL(wp) ::   ldocz            !: Ligand production ratio by grazing
@@ -68,4 +69,5 @@
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  neln  !: number of T-levels + 1 in the euphotic layer
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  heup  !: euphotic layer depth
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::  strn  !: Day duration in hours
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  etot  !: par (photosynthetic available radiation)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  etot_ndcy      !: PAR over 24h in case of diurnal cycle
@@ -74,5 +76,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  epico          !: PAR for pico
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  epicom         !: mean PAR for pico
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  emoy           !: averaged PAR in the mixed layer
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  emoy, etotm    !: averaged PAR in the mixed layer
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  heup_01 !: Absolute euphotic layer depth
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  xksi  !:  Half-saturation con,stant for diatoms
@@ -98,4 +100,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodgoc    !: GOC production
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   consgoc    !: GOC consumption
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   consfe3    !: GOC consumption
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   blim       !: bacterial production factor
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizen      !: size of nanophyto
@@ -118,5 +121,4 @@
    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
 
    !!----------------------------------------------------------------------
@@ -135,8 +137,7 @@
       !!----------------------------------------------------------------------
       ierr(:) = 0
-      ALLOCATE( etot(jpi,jpj,jpk), neln(jpi,jpj), heup(jpi,jpj),    &
-        &       heup_01(jpi,jpj) , xksi(jpi,jpj)               ,  STAT=ierr(1) )
-      !
-  
+      ALLOCATE( etot(jpi,jpj,jpk), neln(jpi,jpj), heup(jpi,jpj),   &
+        &       heup_01(jpi,jpj) , xksi(jpi,jpj), strn(jpi,jpj),  STAT=ierr(1) )
+
       IF( ln_p4z .OR. ln_p5z ) THEN
 
@@ -144,5 +145,6 @@
          ALLOCATE(  enano(jpi,jpj,jpk)    , ediat(jpi,jpj,jpk) ,   &
            &        enanom(jpi,jpj,jpk)   , ediatm(jpi,jpj,jpk),   &
-           &        etot_ndcy(jpi,jpj,jpk), emoy(jpi,jpj,jpk)  ,  STAT=ierr(2) ) 
+           &        etot_ndcy(jpi,jpj,jpk), emoy(jpi,jpj,jpk)  ,   &
+           &        etotm(jpi,jpj,jpk),                           STAT=ierr(2) )
 
          !* Biological SMS
@@ -155,5 +157,5 @@
             &      prodpoc(jpi,jpj,jpk) , conspoc(jpi,jpj,jpk) ,    &
             &      prodgoc(jpi,jpj,jpk) , consgoc(jpi,jpj,jpk) ,    &
-            &      blim   (jpi,jpj,jpk) ,                         STAT=ierr(4) )
+            &      blim   (jpi,jpj,jpk) , consfe3(jpi,jpj,jpk) ,  STAT=ierr(4) )
 
          !* Carbonate chemistry
@@ -164,6 +166,5 @@
          !
          !* Temperature dependency of SMS terms
-         ALLOCATE( tgfunc (jpi,jpj,jpk) , tgfunc2(jpi,jpj,jpk),     &
-            &      tgfunc3(jpi,jpj,jpk) ,                         STAT=ierr(6) )
+         ALLOCATE( tgfunc (jpi,jpj,jpk) , tgfunc2(jpi,jpj,jpk),   STAT=ierr(6) )
          !
          !* Sinking speed
@@ -174,8 +175,5 @@
            &       sizena(jpi,jpj,jpk), sizeda(jpi,jpj,jpk),      STAT=ierr(8) )
          ! 
-         !* Prognostic ligand
-         IF( ln_ligand ) THEN
-           ALLOCATE( plig(jpi,jpj,jpk)  ,                         STAT=ierr(9) )
-         ENDIF
+         ALLOCATE( plig(jpi,jpj,jpk)  ,                           STAT=ierr(9) )
       ENDIF
       !
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/TRP/trcrad.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/TRP/trcrad.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/TRP/trcrad.F90	(revision 14276)
@@ -107,5 +107,5 @@
       ENDIF
       !
-      ALLOCATE( gainmass(jptra,2) )
+      ALLOCATE( gainmass(jptra,4) )
       gainmass(:,:) = 0.
       !
@@ -169,7 +169,8 @@
                         IF( ptrb(ji,jj,jk,jn) > 0. ) THEN                         ! use positive values to compensate mass gain
                            zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptrb > 0
+                           gainmass(jn,1) = gainmass(jn,1) + ptrb(ji,jj,jk,jn) * cvol(ji,jj,jk) * ( zcoef - 1.0 )  ! we are correcting mass...
                            ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * zcoef
                            IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value
-                              gainmass(jn,1) = gainmass(jn,1) - ptrb(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass...
+                              gainmass(jn,2) = gainmass(jn,2) - ptrb(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass...
                               ptrb(ji,jj,jk,jn) = 0.                              ! limit the compensation to keep positive value
                            ENDIF
@@ -190,9 +191,11 @@
          IF( kt == nitend ) THEN
             CALL mpp_sum( 'trcrad', gainmass(:,1) )
+            CALL mpp_sum( 'trcrad', gainmass(:,2) )
             DO jn = jp_sms0, jp_sms1
-               IF( gainmass(jn,1) > 0. ) THEN
+               IF( gainmass(jn,1) < 0. ) THEN
                   ztotmass = glob_sum( 'trcrad', ptrb(:,:,:,jn) * cvol(:,:,:) )
-                  IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn  &
-                     &        , ' total mass : ', ztotmass, ', mass gain : ',  gainmass(jn,1)
+                  IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn  &
+                     &        , ' total mass : ', ztotmass, ', mass correction : ',  gainmass(jn,1)   &
+                     &        , ' mass gain : ', gainmass(jn,2)
                END IF
             END DO
@@ -219,7 +222,8 @@
                         IF( ptrn(ji,jj,jk,jn) > 0. ) THEN                         ! use positive values to compensate mass gain
                            zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn)       ! ztrpos > 0 as ptrb > 0
+                           gainmass(jn,3) = gainmass(jn,3) + ptrn(ji,jj,jk,jn) * cvol(ji,jj,jk) * ( zcoef - 1.0 )  ! we are correcting mass...
                            ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * zcoef
                            IF( zcoef < 0. ) THEN                                  ! if the compensation exceed the positive value
-                              gainmass(jn,2) = gainmass(jn,2) - ptrn(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass...
+                              gainmass(jn,4) = gainmass(jn,4) - ptrn(ji,jj,jk,jn) * cvol(ji,jj,jk)   ! we are adding mass...
                               ptrn(ji,jj,jk,jn) = 0.                              ! limit the compensation to keep positive value
                            ENDIF
@@ -239,10 +243,12 @@
  
          IF( kt == nitend ) THEN
-            CALL mpp_sum( 'trcrad', gainmass(:,2) )
+            CALL mpp_sum( 'trcrad', gainmass(:,3) )
+            CALL mpp_sum( 'trcrad', gainmass(:,4) )
             DO jn = jp_sms0, jp_sms1
-               IF( gainmass(jn,2) > 0. ) THEN
+               IF( gainmass(jn,3) < 0. ) THEN
                   ztotmass = glob_sum( 'trcrad', ptrn(:,:,:,jn) * cvol(:,:,:) )
-                  WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrn, traceur ', jn  &
-                     &        , ' total mass : ', ztotmass, ', mass gain : ',  gainmass(jn,1)
+                  IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16, a, D23.16)') 'trcrad ptrn, traceur ', jn  &
+                     &        , ' total mass : ', ztotmass, ', mass correction : ',  gainmass(jn,3)   &
+                     &        , ' mass gain : ', gainmass(jn,4)
                END IF
             END DO
Index: /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcnam.F90
===================================================================
--- /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcnam.F90	(revision 14275)
+++ /NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/trcnam.F90	(revision 14276)
@@ -248,5 +248,5 @@
       !!---------------------------------------------------------------------
 #if defined key_trdmxl_trc  || defined key_trdtrc
-      INTEGER  ::   ios, ierr                 ! Local integer
+      INTEGER  ::   jn, ios, ierr                 ! Local integer
       !!
       NAMELIST/namtrc_trd/ nn_trd_trc, nn_ctls_trc, rn_ucf_trc, &
@@ -260,4 +260,5 @@
       !
       ALLOCATE( ln_trdtrc(jptra) ) 
+      ln_trdtrc(:) = .false.
       !
       REWIND( numnat_ref )              ! Namelist namtrc_trd in reference namelist : Passive tracer trends
@@ -267,5 +268,5 @@
       READ  ( numnat_cfg, namtrc_trd, IOSTAT = ios, ERR = 906 )
 906   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtrc_trd in configuration namelist' )
-      IF(lwm) WRITE( numont, namtrc_trd )
+     IF(lwm) WRITE( numont, namtrc_trd )
 
       IF(lwp) THEN
