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 13233)
+++ 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) )
          !
