diff --git a/cfgs/SHARED/field_def_nemo-pisces.xml b/cfgs/SHARED/field_def_nemo-pisces.xml
index d73853f06c74506f33d5eb346555ee0588a91d03..17028a230b458794ec9601aebc3a264ac5602042 100644
--- a/cfgs/SHARED/field_def_nemo-pisces.xml
+++ b/cfgs/SHARED/field_def_nemo-pisces.xml
@@ -183,6 +183,9 @@
     <field id="PPNEWN"      long_name="New Primary production of nanophyto"     unit="molC/m3/s"  grid_ref="grid_T_3D" />
     <field id="PPNEWP"      long_name="New Primary production of picophyto"     unit="molC/m3/s"  grid_ref="grid_T_3D" />
     <field id="PPNEWD"      long_name="New Primary production of diatoms"       unit="molC/m3/s"  grid_ref="grid_T_3D" />
+    <field id="GPPHYN"      long_name="Gross Primary production of nanophyto"   unit="molC/m3/s"  grid_ref="grid_T_3D" />
+    <field id="GPPHYP"      long_name="Gross Primary production of picophyto"   unit="molC/m3/s"  grid_ref="grid_T_3D" />
+    <field id="GPPHYD"      long_name="Gross Primary production of diatoms"     unit="molC/m3/s"  grid_ref="grid_T_3D" />
     <field id="PBSi"        long_name="Primary production of Si diatoms"        unit="molC/m3/s"  grid_ref="grid_T_3D" />
     <field id="PFeN"        long_name="Primary production of nano iron"         unit="molC/m3/s"  grid_ref="grid_T_3D" />
     <field id="PFeP"        long_name="Primary production of pico iron"         unit="molC/m3/s"  grid_ref="grid_T_3D" />
diff --git a/cfgs/SHARED/namelist_pisces_ref b/cfgs/SHARED/namelist_pisces_ref
index 10aff870eab7ba858ad399fb796d591282357b67..25d0e9f9849b09cdc4794d8131f3a64017eba415 100644
--- a/cfgs/SHARED/namelist_pisces_ref
+++ b/cfgs/SHARED/namelist_pisces_ref
@@ -152,7 +152,6 @@
 !-----------------------------------------------------------------------
    pislopen   =  2.       ! P-I slope
    pisloped   =  2.       ! P-I slope  for diatoms
-   xadap      =  0.       ! Adaptation factor to low light
    excretn    =  0.05     ! excretion ratio of phytoplankton
    excretd    =  0.05     ! excretion ratio of diatoms
    bresp      =  0.033    ! Basal respiration rate
@@ -172,7 +171,6 @@
    excretn    =  0.05     ! excretion ratio of phytoplankton
    excretp    =  0.05     ! excretion ratio of picophytoplankton
    excretd    =  0.05     ! excretion ratio of diatoms
-   xadap      =  0.       ! Adaptation factor to low light
    bresp      =  0.02     ! Basal respiration rate
    thetannm   =  0.3      ! Maximum Chl/N in nanophytoplankton
    thetanpm   =  0.3      ! Maximum Chl/N in picophytoplankton
diff --git a/src/TOP/PISCES/P4Z/p4zfechem.F90 b/src/TOP/PISCES/P4Z/p4zfechem.F90
index 2598766485bdad722d29f17aaab666676789cfb7..b5d56529a0d1e1af6373a495192becc37722b72d 100644
--- a/src/TOP/PISCES/P4Z/p4zfechem.F90
+++ b/src/TOP/PISCES/P4Z/p4zfechem.F90
@@ -82,13 +82,16 @@ CONTAINS
       ! Parameterization from Pham and Ito (2018)
       ! -------------------------------------------------
       xfecolagg(:,:,:) = ligand * 1E9 + MAX(0., chemo2(:,:,:) - tr(:,:,:,jpoxy,Kbb) ) / 400.E-6
-      IF( ln_ligvar ) THEN
-         ztotlig(:,:,:) =  0.09 * 0.667 * tr(:,:,:,jpdoc,Kbb) * 1E6 + xfecolagg(:,:,:)
-         ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
-      ELSE
-        IF( ln_ligand ) THEN  ;   ztotlig(:,:,:) = tr(:,:,:,jplgw,Kbb) * 1E9
-        ELSE                  ;   ztotlig(:,:,:) = ligand * 1E9 
-        ENDIF
+      !
+      IF( ln_ligand ) THEN  
+         ztotlig(:,:,:) = tr(:,:,:,jplgw,Kbb) * 1E9
+      ELSE 
+         IF( ln_ligvar ) THEN
+            ztotlig(:,:,:) =  0.09 * 0.667 * tr(:,:,:,jpdoc,Kbb) * 1E6 + xfecolagg(:,:,:)
+            ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
+         ELSE
+            ztotlig(:,:,:) = ligand * 1E9 
+         ENDIF
       ENDIF
 
       ! ------------------------------------------------------------
@@ -124,13 +127,22 @@ CONTAINS
       ! prognostic or non variable, all the colloidal fraction is supposed
       ! to coagulate
       ! ----------------------------------------------------------------------
-      IF (ln_ligand) THEN
-         zfecoll(:,:,:) = 0.5 * zFeL1(:,:,:) * MAX(0., tr(:,:,:,jplgw,Kbb) - xfecolagg(:,:,:) * 1.0E-9 ) / ( tr(:,:,:,jplgw,Kbb) + rtrn ) 
+      IF( ln_ligand ) THEN
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+            zfecoll(ji,jj,jk) = 0.5 * zFeL1(ji,jj,jk) * MAX(0., ztotlig(ji,jj,jk) - xfecolagg(ji,jj,jk) ) &
+                  &              / ( ztotlig(ji,jj,jk) + rtrn )
+         END_3D
       ELSE
-         IF (ln_ligvar) THEN
-            zfecoll(:,:,:) = 0.5 * zFeL1(:,:,:) * MAX(0., tr(:,:,:,jplgw,Kbb) - xfecolagg(:,:,:) * 1.0E-9 ) / ( tr(:,:,:,jplgw,Kbb) + rtrn )   
+         IF( ln_ligvar ) THEN
+            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+               zfecoll(ji,jj,jk) = 0.5 * zFeL1(ji,jj,jk) * MAX(0., ztotlig(ji,jj,jk) - xfecolagg(ji,jj,jk) ) &
+                  &              / ( ztotlig(ji,jj,jk) + rtrn )
+            END_3D
          ELSE
-            zfecoll(:,:,:) = 0.5 * zFeL1(:,:,:)
+            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+               zfecoll(ji,jj,jk) = 0.5 * zFeL1(ji,jj,jk) * MAX(0., 0.09 * 0.667 * tr(ji,jj,jk,jpdoc,Kbb) * 1E6 ) &
+                  &             / ( ztotlig(ji,jj,jk) + rtrn )
+            END_3D
          ENDIF
       ENDIF
 
diff --git a/src/TOP/PISCES/P4Z/p4zlim.F90 b/src/TOP/PISCES/P4Z/p4zlim.F90
index 4e227808aad4e067b9348d662bcf0637977d7b93..6ea372375746cc5bd9ab004f9234536153168436 100644
--- a/src/TOP/PISCES/P4Z/p4zlim.F90
+++ b/src/TOP/PISCES/P4Z/p4zlim.F90
@@ -242,7 +242,7 @@ CONTAINS
 
          xfracal(ji,jj,jk) = caco3r * MIN( zlim1, zlim2, zlim3 )                  &
             &                       * ztem1 / ( 0.1 + ztem1 )                     &
-            &                       * MAX( 1., tr(ji,jj,jk,jpphy,Kbb) * 1.e6 / 2. )  &
+            &                       * MAX( 1., tr(ji,jj,jk,jpphy,Kbb) * xsizephy )  &
             &                       * zetot1 * zetot2               &
             &                       * ( 1. + EXP(-ztem2 * ztem2 / 25. ) )         &
             &                       * MIN( 1., 50. / ( hmld(ji,jj) + rtrn ) )
@@ -250,18 +250,6 @@ CONTAINS
          xfracal(ji,jj,jk) = MAX( 0.02, xfracal(ji,jj,jk) )
       END_3D
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-         ! denitrification factor computed from O2 levels
-         nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - tr(ji,jj,jk,jpoxy,Kbb) )    &
-            &                                / ( oxymin + tr(ji,jj,jk,jpoxy,Kbb) )  )
-         nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
-         !
-         ! denitrification factor computed from NO3 levels
-         nitrfac2(ji,jj,jk) = MAX( 0.e0,       ( 1.E-6 - tr(ji,jj,jk,jpno3,Kbb) )  &
-            &                                / ( 1.E-6 + tr(ji,jj,jk,jpno3,Kbb) ) )
-         nitrfac2(ji,jj,jk) = MIN( 1., nitrfac2(ji,jj,jk) )
-      END_3D
-      !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics
         CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! euphotic layer deptht
         CALL iom_put( "LNnut"  , xlimphy(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
diff --git a/src/TOP/PISCES/P4Z/p4zpoc.F90 b/src/TOP/PISCES/P4Z/p4zpoc.F90
index 54b8f202e626b8367f6da45d6a2181337bf2f339..bc0ae47026fdc95c0c9b5b7472811e2ca93c8232 100644
--- a/src/TOP/PISCES/P4Z/p4zpoc.F90
+++ b/src/TOP/PISCES/P4Z/p4zpoc.F90
@@ -118,7 +118,7 @@ CONTAINS
      ! a standard parameterisation with a constant lability
      ! -----------------------------------------------------------------------
      ztremint(:,:,:) = zremigoc(:,:,:)
-     DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+     DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1)
         IF (tmask(ji,jj,jk) == 1.) THEN
           zdep = hmld(ji,jj)
           !
diff --git a/src/TOP/PISCES/P4Z/p4zprod.F90 b/src/TOP/PISCES/P4Z/p4zprod.F90
index 59ca90ccc62b5933ddad31ffff1460965cfc528d..a0c736a2f49b2e8c7446dd2339b6ef37c6d1c6a4 100644
--- a/src/TOP/PISCES/P4Z/p4zprod.F90
+++ b/src/TOP/PISCES/P4Z/p4zprod.F90
@@ -27,7 +27,6 @@ MODULE p4zprod
 
    REAL(wp), PUBLIC ::   pislopen     !:  P-I slope of nanophytoplankton
    REAL(wp), PUBLIC ::   pisloped     !:  P-I slope of diatoms
-   REAL(wp), PUBLIC ::   xadap        !:  Adaptation factor to low light 
    REAL(wp), PUBLIC ::   excretn      !:  Excretion ratio of nanophyto
    REAL(wp), PUBLIC ::   excretd      !:  Excretion ratio of diatoms
    REAL(wp), PUBLIC ::   bresp        !:  Basal respiration rate
@@ -69,8 +68,8 @@ CONTAINS
       INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices
       !
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zsilfac, znanotot, zdiattot, zconctemp, zconctemp2
-      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap, zlim, zsiborn
+      REAL(wp) ::   zsilfac, znanotot, zdiattot
+      REAL(wp) ::   zratio, zmax, zsilim, ztn, zlim, zsiborn
       REAL(wp) ::   zpptot, zpnewtot, zpregtot, zprochln, zprochld
       REAL(wp) ::   zproddoc, zprodsil, zprodfer, zprodlig
       REAL(wp) ::   zpislopen, zpisloped, zfact
@@ -112,8 +111,9 @@ CONTAINS
 
          DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
             IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+               zval = 24.0
                IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
-                  zval = MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
+                  zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
                ENDIF
                zmxl_chl(ji,jj,jk) = zval / 24.
                zmxl_fac(ji,jj,jk) = 1.0 - exp( -0.26 * zval )
@@ -143,23 +143,19 @@ CONTAINS
       ! -----------------------------------------------------------------------
       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-            ztn         = MAX( 0., ts(ji,jj,jk,jp_tem,Kmm) - 15. )
-            zadap       = xadap * ztn / ( 2.+ ztn )
-            zconctemp   = MAX( 0.e0 , tr(ji,jj,jk,jpdia,Kbb) - xsizedia )
-            zconctemp2  = tr(ji,jj,jk,jpdia,Kbb) - zconctemp
             !
             ! The initial slope of the PI curve can be increased for nano
             ! to account for photadaptation, for instance in the DCM
             ! This parameterization is adhoc and should be either 
             ! improved or removed in future versions of the model
-
             ! Nanophytoplankton
-            zpislopeadn(ji,jj,jk) = pislopen * ( 1.+ zadap  * EXP( -0.25 * enano(ji,jj,jk) ) )  &
-            &                   * tr(ji,jj,jk,jpnch,Kbb) /( tr(ji,jj,jk,jpphy,Kbb) * 12. + rtrn)
+            zpislopeadn(ji,jj,jk) = pislopen * tr(ji,jj,jk,jpnch,Kbb)  &
+            &                   /( tr(ji,jj,jk,jpphy,Kbb) * 12. + rtrn)
 
             ! Diatoms
-            zpislopeadd(ji,jj,jk) = (pislopen * zconctemp2 + pisloped * zconctemp) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )   &
-            &                   * tr(ji,jj,jk,jpdch,Kbb) /( tr(ji,jj,jk,jpdia,Kbb) * 12. + rtrn)
+            zpislopeadd(ji,jj,jk) = pisloped * tr(ji,jj,jk,jpdch,Kbb)   &
+            &                   /( tr(ji,jj,jk,jpdia,Kbb) * 12. + rtrn)
+
          ENDIF
       END_3D
 
@@ -435,7 +431,7 @@ CONTAINS
       INTEGER ::   ios   ! Local integer
       !
       ! Namelist block
-      NAMELIST/namp4zprod/ pislopen, pisloped, xadap, bresp, excretn, excretd,  &
+      NAMELIST/namp4zprod/ pislopen, pisloped, bresp, excretn, excretd,  &
          &                 chlcnm, chlcdm, chlcmin, fecnm, fecdm, grosip
       !!----------------------------------------------------------------------
       !
@@ -456,7 +452,6 @@ CONTAINS
          WRITE(numout,*) '   Namelist : namp4zprod'
          WRITE(numout,*) '      mean Si/C ratio                           grosip       =', grosip
          WRITE(numout,*) '      P-I slope                                 pislopen     =', pislopen
-         WRITE(numout,*) '      Acclimation factor to low light           xadap        =', xadap
          WRITE(numout,*) '      excretion ratio of nanophytoplankton      excretn      =', excretn
          WRITE(numout,*) '      excretion ratio of diatoms                excretd      =', excretd
          WRITE(numout,*) '      basal respiration in phytoplankton        bresp        =', bresp
diff --git a/src/TOP/PISCES/P4Z/p4zrem.F90 b/src/TOP/PISCES/P4Z/p4zrem.F90
index c63a6e63f59e3fd8d803181902ac7a046a6950e1..126f1e57f16c23f30a7fb992e304c933681071f7 100644
--- a/src/TOP/PISCES/P4Z/p4zrem.F90
+++ b/src/TOP/PISCES/P4Z/p4zrem.F90
@@ -95,15 +95,13 @@ CONTAINS
       ! recent version of PISCES with bacteria.
       ! ----------------------------------------------------------------
       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-         zdep = MAX( hmld(ji,jj), heup_01(ji,jj) )
+         zdep = MAX( hmld(ji,jj), heup_01(ji,jj), gdept(ji,jj,1,Kmm) )
          IF ( gdept(ji,jj,jk,Kmm) < zdep ) THEN
             zdepbac(ji,jj,jk) = 0.6 * ( MAX(0.0, tr(ji,jj,jk,jpzoo,Kbb) + tr(ji,jj,jk,jpmes,Kbb) ) * 1.0E6 )**0.6 * 1.E-6
             ztempbac(ji,jj)   = zdepbac(ji,jj,jk)
-!         IF( gdept(ji,jj,jk,Kmm) >= zdep ) THEN
          ELSE
             zdepmin = MIN( 1., zdep / gdept(ji,jj,jk,Kmm) )
             zdepbac (ji,jj,jk) = zdepmin**0.683 * ztempbac(ji,jj)
-!            zdepeff(ji,jj,jk) = zdepeff(ji,jj,jk) * zdepmin**0.3
             zdepeff(ji,jj,jk) = zdepeff(ji,jj,jk) * zdepmin**0.6
          ENDIF
       END_3D
@@ -180,8 +178,8 @@ CONTAINS
          ! Bacteries are obliged to take up iron from the water. Some
          ! studies (especially at Papa) have shown this uptake to be significant
          ! ----------------------------------------------------------
-         zbactfer = feratb * 0.6_wp * xstep * tgfunc(ji,jj,jk) * xlimbacl(ji,jj,jk) * tr(ji,jj,jk,jpfer,Kbb)    &
-           &       / ( xkferb + tr(ji,jj,jk,jpfer,Kbb) ) * zdepeff(ji,jj,jk) * zdepbac(ji,jj,jk)
+         zbactfer = feratb * 0.6_wp * xstep * tgfunc(ji,jj,jk) * xlimbacl(ji,jj,jk) * biron(ji,jj,jk)    &
+           &       / ( xkferb + biron(ji,jj,jk) ) * 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 
diff --git a/src/TOP/PISCES/P4Z/p4zsed.F90 b/src/TOP/PISCES/P4Z/p4zsed.F90
index ca6bdb281c10b0e6c24143440ce85618e10042f0..f6ad15c38ac8ec1b2af53b71403259111b5d9721 100644
--- a/src/TOP/PISCES/P4Z/p4zsed.F90
+++ b/src/TOP/PISCES/P4Z/p4zsed.F90
@@ -118,7 +118,7 @@ CONTAINS
                 !
               zflx = (  tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj)   &
                 &     + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) ) * 1E6
-              zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2
+              zbureff = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2 * MIN(gdepw(ji,jj,ikt+1,Kmm) / 1000.00, 1.0)
            ENDIF
          END_2D
          !
@@ -251,7 +251,7 @@ CONTAINS
             zfact = zlim * rfact2
             ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
             ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) )
-            ztrdop(ji,jj,jk) = tr(ji,jj,jk,jpdop,Kbb) / ( 1E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (1. - ztrpo4(ji,jj,jk))
+            ztrdop(ji,jj,jk) = tr(ji,jj,jk,jpdop,Kbb) / ( 10E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (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)
          END_3D
@@ -288,7 +288,7 @@ CONTAINS
             tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zfact * 1.0 / 3.0
             tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zfact * 1.0 / 3.0
             tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + 16.0 / 46.0 * zfact / 3.0  &
-            &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   &
+            &                     - 16.0 / 46.0 * zfact  * 2.0 / 3.0  * ztrdop(ji,jj,jk)   &
             &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
             tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zfact * 1.0 / 3.0 * 2.0 / 3.0
             tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) + zfact * 1.0 / 3.0 * 2.0 /3.0
diff --git a/src/TOP/PISCES/P4Z/p5zlim.F90 b/src/TOP/PISCES/P4Z/p5zlim.F90
index 91e94183a3e6f79b83a9db8ac4c3129fc5be99be..f4528230a8e46e713e1b795a1714286cd06639ab 100644
--- a/src/TOP/PISCES/P4Z/p5zlim.F90
+++ b/src/TOP/PISCES/P4Z/p5zlim.F90
@@ -482,13 +482,6 @@ CONTAINS
          xfracal(ji,jj,jk) = MAX( 0.02, MIN( 0.8 , xfracal(ji,jj,jk) ) )
       END_3D
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-         ! denitrification factor computed from O2 levels
-         nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - tr(ji,jj,jk,jpoxy,Kbb) )    &
-            &                                / ( oxymin + tr(ji,jj,jk,jpoxy,Kbb) )  )
-         nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
-      END_3D
-      !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics
         CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! euphotic layer deptht
         CALL iom_put( "LNnut"  , xlimphy(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
@@ -610,17 +603,18 @@ CONTAINS
       xpsinh4  = 1.8 * rno3
       xpsiuptk = 1.0 / 6.625
       !
-      nitrfac(:,:,jpk) = 0._wp
-      xfracal(:,:,jpk) = 0._wp
-      xlimphy(:,:,jpk) = 0._wp
-      xlimpic(:,:,jpk) = 0._wp
-      xlimdia(:,:,jpk) = 0._wp
-      xlimnfe(:,:,jpk) = 0._wp
-      xlimpfe(:,:,jpk) = 0._wp
-      xlimdfe(:,:,jpk) = 0._wp
-      sizen  (:,:,jpk) = 0._wp
-      sizep  (:,:,jpk) = 0._wp
-      sized  (:,:,jpk) = 0._wp
+      nitrfac (:,:,jpk) = 0._wp
+      nitrfac2(:,:,jpk) = 0._wp
+      xfracal (:,:,jpk) = 0._wp
+      xlimphy (:,:,jpk) = 0._wp
+      xlimpic (:,:,jpk) = 0._wp
+      xlimdia (:,:,jpk) = 0._wp
+      xlimnfe (:,:,jpk) = 0._wp
+      xlimpfe (:,:,jpk) = 0._wp
+      xlimdfe (:,:,jpk) = 0._wp
+      sizen   (:,:,jpk) = 0._wp
+      sizep   (:,:,jpk) = 0._wp
+      sized   (:,:,jpk) = 0._wp
       !
    END SUBROUTINE p5z_lim_init
 
diff --git a/src/TOP/PISCES/P4Z/p5zmeso.F90 b/src/TOP/PISCES/P4Z/p5zmeso.F90
index 24ca0260d8b5f0d17aebf1da57293231fb8eebf1..1506928c5f115818c63fec208a33799df8fa87a5 100644
--- a/src/TOP/PISCES/P4Z/p5zmeso.F90
+++ b/src/TOP/PISCES/P4Z/p5zmeso.F90
@@ -106,10 +106,10 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing, zfezoo2
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrarem, zgraref, zgrapoc, zgrapof
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrarep, zgraren, zgrapon, zgrapop
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgradoc, zgradon, zgradop
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgradoc, zgradon, zgradop, zgrabsi
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigrem, zgramigref, zgramigpoc, zgramigpof
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigrep, zgramigren, zgramigpop, zgramigpon
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigdoc, zgramigdop, zgramigdon
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigdoc, zgramigdop, zgramigdon, zgramigbsi
       
 
       !!---------------------------------------------------------------------
@@ -187,7 +187,7 @@ CONTAINS
          ! have low abundance, .i.e. zooplankton become less specific 
          ! to avoid starvation.
          ! ----------------------------------------------------------
-         zsigma = 1.0 - zdenom**3/(0.1**3+zdenom**3)
+         zsigma = 1.0 - zdenom**3/(0.05**3+zdenom**3)
          zsigma = xsigma2 + xsigma2del * zsigma
          ! Nanophytoplankton and diatoms are the only preys considered
          ! to be close enough to have potential interference
@@ -366,7 +366,7 @@ CONTAINS
          !   sinks
          !   --------------------------------------------------------------
          tr(ji,jj,jk,jpmes,Krhs) = tr(ji,jj,jk,jpmes,Krhs) + zepsherv * zgraztotc - zrespirc   &
-         &                     - ztortz - zgrazm
+         &                       - ztortz - zgrazm
          tr(ji,jj,jk,jpdia,Krhs) = tr(ji,jj,jk,jpdia,Krhs) - zgrazdc
          tr(ji,jj,jk,jpndi,Krhs) = tr(ji,jj,jk,jpndi,Krhs) - zgrazdn
          tr(ji,jj,jk,jppdi,Krhs) = tr(ji,jj,jk,jppdi,Krhs) - zgrazdp
@@ -379,7 +379,7 @@ CONTAINS
          tr(ji,jj,jk,jpnch,Krhs) = tr(ji,jj,jk,jpnch,Krhs) - zgraznc * tr(ji,jj,jk,jpnch,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) + rtrn )
          tr(ji,jj,jk,jpdch,Krhs) = tr(ji,jj,jk,jpdch,Krhs) - zgrazdc * tr(ji,jj,jk,jpdch,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
          tr(ji,jj,jk,jpdsi,Krhs) = tr(ji,jj,jk,jpdsi,Krhs) - zgrazdc * tr(ji,jj,jk,jpdsi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
-         tr(ji,jj,jk,jpgsi,Krhs) = tr(ji,jj,jk,jpgsi,Krhs) + zgrazdc * tr(ji,jj,jk,jpdsi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
+         zgrabsi(ji,jj,jk)       = zgrazdc * tr(ji,jj,jk,jpdsi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
 
          tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) - zgrazpoc - zgrazffep + zfracc
          prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zfracc
@@ -415,13 +415,13 @@ CONTAINS
       IF( ln_dvm_meso ) THEN
           ALLOCATE( zgramigrem(jpi,jpj), zgramigref(jpi,jpj), zgramigpoc(jpi,jpj), zgramigpof(jpi,jpj) )
           ALLOCATE( zgramigrep(jpi,jpj), zgramigren(jpi,jpj), zgramigpop(jpi,jpj), zgramigpon(jpi,jpj) )
-          ALLOCATE( zgramigdoc(jpi,jpj), zgramigdon(jpi,jpj), zgramigdop(jpi,jpj) )
+          ALLOCATE( zgramigdoc(jpi,jpj), zgramigdon(jpi,jpj), zgramigdop(jpi,jpj), zgramigbsi(jpi,jpj) )
           zgramigrem(:,:)  = 0.0   ;   zgramigref(:,:) = 0.0
           zgramigrep(:,:)  = 0.0   ;   zgramigren(:,:) = 0.0
           zgramigpoc(:,:)  = 0.0   ;   zgramigpof(:,:) = 0.0
           zgramigpop(:,:)  = 0.0   ;   zgramigpon(:,:) = 0.0
           zgramigdoc(:,:)  = 0.0   ;   zgramigdon(:,:) = 0.0
-          zgramigdop(:,:)  = 0.0   
+          zgramigdop(:,:)  = 0.0   ;   zgramigbsi(:,:) = 0.0 
                  
           ! Compute the amount of materials that will go into vertical migration
           ! This fraction is sumed over the euphotic zone and is removed from 
@@ -442,6 +442,8 @@ CONTAINS
                 zgramigdoc(ji,jj) = zgramigdoc(ji,jj) + xfracmig * zgradoc(ji,jj,jk) * zmigthick
                 zgramigdop(ji,jj) = zgramigdop(ji,jj) + xfracmig * zgradop(ji,jj,jk) * zmigthick
                 zgramigdon(ji,jj) = zgramigdon(ji,jj) + xfracmig * zgradon(ji,jj,jk) * zmigthick
+                zgramigbsi(ji,jj) = zgramigbsi(ji,jj) + xfracmig * zgrabsi(ji,jj,jk) * zmigthick
+
 
                 zgrarem(ji,jj,jk)  = zgrarem(ji,jj,jk) * ( (1.0 - xfracmig) + xfracmig * zmigreltime )
                 zgrarep(ji,jj,jk)  = zgrarep(ji,jj,jk) * ( (1.0 - xfracmig) + xfracmig * zmigreltime )
@@ -454,6 +456,7 @@ CONTAINS
                 zgradoc(ji,jj,jk)  = zgradoc(ji,jj,jk) * ( (1.0 - xfracmig) + xfracmig * zmigreltime )
                 zgradop(ji,jj,jk)  = zgradop(ji,jj,jk) * ( (1.0 - xfracmig) + xfracmig * zmigreltime )
                 zgradon(ji,jj,jk)  = zgradon(ji,jj,jk) * ( (1.0 - xfracmig) + xfracmig * zmigreltime )
+                zgrabsi(ji,jj,jk)  = zgrabsi(ji,jj,jk) * ( (1.0 - xfracmig) + xfracmig * zmigreltime )
              ENDIF
           END_3D
 
@@ -475,6 +478,7 @@ CONTAINS
                  zgradoc(ji,jj,jkt) = zgradoc(ji,jj,jkt) + zgramigdoc(ji,jj) * zdep
                  zgradop(ji,jj,jkt) = zgradop(ji,jj,jkt) + zgramigdop(ji,jj) * zdep 
                  zgradon(ji,jj,jkt) = zgradon(ji,jj,jkt) + zgramigdon(ji,jj) * zdep 
+                 zgrabsi(ji,jj,jkt) = zgrabsi(ji,jj,jkt) + zgramigbsi(ji,jj) * zdep
               ENDIF
           END_2D
                    !
@@ -482,7 +486,7 @@ CONTAINS
           ! ------------------------------
           DEALLOCATE( zgramigrem, zgramigref, zgramigpoc, zgramigpof )
           DEALLOCATE( zgramigrep, zgramigren, zgramigpop, zgramigpon )
-          DEALLOCATE( zgramigdoc, zgramigdon, zgramigdop )
+          DEALLOCATE( zgramigdoc, zgramigdon, zgramigdop, zgramigbsi )
          ! End of the ln_dvm_meso part
       ENDIF
 
@@ -510,6 +514,7 @@ CONTAINS
          tr(ji,jj,jk,jpgon,Krhs) = tr(ji,jj,jk,jpgon,Krhs) + zgrapon(ji,jj,jk)
          tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) + zgrapop(ji,jj,jk)
          tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + zgrapof(ji,jj,jk)
+         tr(ji,jj,jk,jpgsi,Krhs) = tr(ji,jj,jk,jpgsi,Krhs) + zgrabsi(ji,jj,jk)
       END_3D
       !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
diff --git a/src/TOP/PISCES/P4Z/p5zmicro.F90 b/src/TOP/PISCES/P4Z/p5zmicro.F90
index 6edfde55e504a8916bb3098a293078c1da86c3c0..97459e7375f8a658178ce812f63f23969af20cac 100644
--- a/src/TOP/PISCES/P4Z/p5zmicro.F90
+++ b/src/TOP/PISCES/P4Z/p5zmicro.F90
@@ -242,7 +242,7 @@ CONTAINS
          ! Excess carbon in the food is used preferentially
          ! when activated by zmetexcess
          ! ------------------------------------------------
-         zexcess  = zgraztotc * zepsherf * (1.0 - zepshert) * zmetexcess
+         zexcess  = zgraztotc * zepsherq * zepsherf * (1.0 - zepshert) * zmetexcess
          zbasresb = MAX(0., zrespz - zexcess)
          zbasresi = zexcess + MIN(0., zrespz - zexcess)  
          zrespirc = srespir * zepsherv * zgraztotc + zbasresb
diff --git a/src/TOP/PISCES/P4Z/p5zprod.F90 b/src/TOP/PISCES/P4Z/p5zprod.F90
index 848f8550539968bc10d6efcfcfc6ae007fc215d2..5dbc757b6de30f9c098e527aa0c2bec6dc5e733c 100644
--- a/src/TOP/PISCES/P4Z/p5zprod.F90
+++ b/src/TOP/PISCES/P4Z/p5zprod.F90
@@ -77,11 +77,11 @@ CONTAINS
       REAL(wp) ::   zsilfac, znanotot, zpicotot, zdiattot, zconctemp, zconctemp2
       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) ::   zlim, zsiborn, zprod, zprontot, zproptot, zprodtot
       REAL(wp) ::   zproddoc, zproddon, zproddop, zprodsil, zprodfer, zprodlig, zresptot     
       REAL(wp) ::   zprnutmax, zprochln, zprochld, zprochlp
       REAL(wp) ::   zpislopen, zpislopep, zpisloped
-      REAL(wp) ::   zval, zpptot, zpnewtot, zpregtot
+      REAL(wp) ::   zval, zpptot, zpnewtot, zpregtot, zpo4tot
       REAL(wp) ::   zqfpmax, zqfnmax, zqfdmax
       REAL(wp) ::   zfact, zrfact2, zmaxsi, zratiosi, zsizetmp, zlimfac, zsilim
       CHARACTER (len=25) :: charout
@@ -96,7 +96,6 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zproregn, zproregp, zproregd
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpropo4n, zpropo4p, zpropo4d
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprodopn, zprodopp, zprodopd
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrespn, zrespp, zrespd
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmxl_fac, zmxl_chl
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zpligprod
       !!---------------------------------------------------------------------
@@ -112,7 +111,6 @@ CONTAINS
       zprdia  (:,:,:) = 0._wp ; zprpic  (:,:,:) = 0._wp ; zprbio  (:,:,:) = 0._wp
       zprodopn(:,:,:) = 0._wp ; zprodopp(:,:,:) = 0._wp ; zprodopd(:,:,:) = 0._wp
       zysopt  (:,:,:) = 0._wp 
-      zrespn  (:,:,:) = 0._wp ; zrespp  (:,:,:) = 0._wp ; zrespd  (:,:,:) = 0._wp 
       zmxl_fac(:,:,:) = 0._wp ; zmxl_chl(:,:,:) = 0._wp
 
       ! Computation of the optimal production rates and nutrient uptake
@@ -134,8 +132,9 @@ CONTAINS
 
          DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
             IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+               zval = 24.
                IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
-                  zval = MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
+                  zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
                ENDIF
                zmxl_chl(ji,jj,jk) = zval / 24.
                zmxl_fac(ji,jj,jk) = 1.0 - exp( -0.26 * zval )
@@ -168,17 +167,14 @@ CONTAINS
       ! The formulation proposed by Geider et al. (1997) has been used.
       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-            ! Computation of the P-I slope for nanos and diatoms
-            ztn         = MAX( 0., ts(ji,jj,jk,jp_tem,Kmm) - 15. )
-            zadap       = xadap * ztn / ( 2.+ ztn )
             !
             ! Nanophytoplankton
             zpislopeadn(ji,jj,jk) = pislopen * tr(ji,jj,jk,jpnch,Kbb)    &
             &                       /( tr(ji,jj,jk,jpphy,Kbb) * 12. + rtrn)
 
             ! Picophytoplankton
-            zpislopeadp(ji,jj,jk) = pislopep * ( 1. + zadap * EXP( -0.25 * epico(ji,jj,jk) ) )   &
-            &                       * tr(ji,jj,jk,jppch,Kbb) /( tr(ji,jj,jk,jppic,Kbb) * 12. + rtrn)
+            zpislopeadp(ji,jj,jk) = pislopep * tr(ji,jj,jk,jppch,Kbb)    &
+            &                       /( tr(ji,jj,jk,jppic,Kbb) * 12. + rtrn)
 
             ! Diatoms
             zpislopeadd(ji,jj,jk) = pisloped * tr(ji,jj,jk,jpdch,Kbb)    &
@@ -191,9 +187,9 @@ CONTAINS
             ! Computation of production function for Carbon
             ! Actual light levels are used here 
             !  ---------------------------------------------
-            zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) / zmxl_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) )  )
+            zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) )
+            zprpic(ji,jj,jk) = zprpic(ji,jj,jk) * ( 1.- EXP( -zpislopep * epico(ji,jj,jk) ) )
+            zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) ) )
 
             !  Computation of production function for Chlorophyll
             !  Mean light level in the mixed layer (when appropriate)
@@ -203,6 +199,7 @@ CONTAINS
             zpislopen = zpislopen * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
             zpisloped = zpisloped * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
             zpislopep = zpislopep * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
+
             zprchln(ji,jj,jk) = zprmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) )  )
             zprchlp(ji,jj,jk) = zprmaxp(ji,jj,jk) * ( 1.- EXP( -zpislopep * epicom(ji,jj,jk) )  )
             zprchld(ji,jj,jk) = zprmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) )  )
@@ -215,7 +212,7 @@ CONTAINS
             ! ------------------------
             ! Si/C increases with iron stress and silicate availability (zsilfac)
             ! Si/C is arbitrariliy increased for very high Si concentrations
-            ! to mimic the very high ratios observed in the Southern Ocean (zsilfac2)
+            ! 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
             ! proposed by Gurney and Davidson (1999).
@@ -224,29 +221,31 @@ CONTAINS
             zsilim = xlimdia(ji,jj,jk) * zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn )
             zsiborn = tr(ji,jj,jk,jpsil,Kbb) * tr(ji,jj,jk,jpsil,Kbb) * tr(ji,jj,jk,jpsil,Kbb)
             IF (gphit(ji,jj) < -30 ) THEN
-              zsilfac2 = 1. + 1. * zsiborn / ( zsiborn + xksi2**3 )
+              zsilfac = 1. + 1. * zsiborn / ( zsiborn + xksi2**3 )
             ELSE
-              zsilfac2 = 1.
+              zsilfac = 1.
             ENDIF
-            zratiosi = 1.0 - tr(ji,jj,jk,jpdsi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn ) / ( zsilfac2 * grosip * 3.0 + rtrn )
+            zratiosi = 1.0 - tr(ji,jj,jk,jpdsi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + 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
       END_3D
 
-      !  Sea-ice effect on production
+      ! Sea-ice effect on production
       ! No production is assumed below sea ice
-      ! -------------------------------------- 
+      ! --------------------------------------
       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-         zprbio(ji,jj,jk)  = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
-         zprpic(ji,jj,jk)  = zprpic(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
-         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) )
+         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+            zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * (1.0 - fr_i(ji,jj) )
+            zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * (1.0 - fr_i(ji,jj) )
+            zprpic(ji,jj,jk) = zprpic(ji,jj,jk) * (1.0 - fr_i(ji,jj) )
+            zprnut(ji,jj,jk) = zprnut(ji,jj,jk) * (1.0 - fr_i(ji,jj) )
+         ENDIF
       END_3D
 
       ! Computation of the various production and uptake terms of nanophytoplankton 
@@ -429,10 +428,12 @@ CONTAINS
 
       !   Update the arrays TRA which contain the biological sources and sinks
       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-        zpptot   = zpropo4n(ji,jj,jk) + zpropo4d(ji,jj,jk) + zpropo4p(ji,jj,jk)
+        zpo4tot  = zpropo4n(ji,jj,jk) + zpropo4d(ji,jj,jk) + zpropo4p(ji,jj,jk)
         zpnewtot = zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk) + zpronewp(ji,jj,jk)
         zpregtot = zproregn(ji,jj,jk) + zproregd(ji,jj,jk) + zproregp(ji,jj,jk)
 
+        zpptot   = zprorcap(ji,jj,jk) + zprorcan(ji,jj,jk) + zprorcad(ji,jj,jk)
+
         zprontot = zpronewn(ji,jj,jk) + zproregn(ji,jj,jk)
         zproptot = zpronewp(ji,jj,jk) + zproregp(ji,jj,jk)
         zprodtot = zpronewd(ji,jj,jk) + zproregd(ji,jj,jk)
@@ -448,17 +449,15 @@ CONTAINS
         zproddon =  excretd * zprodtot + excretn * zprontot + excretp * zproptot
 
         zprodfer = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) + texcretp * zprofep(ji,jj,jk)
-        zresptot = zrespn(ji,jj,jk) + zrespp(ji,jj,jk) + zrespd(ji,jj,jk) 
         !
-        tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - zpptot
+        tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - zpo4tot
         tr(ji,jj,jk,jpno3,Krhs) = tr(ji,jj,jk,jpno3,Krhs) - zpnewtot
         tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) - zpregtot  
         !
         tr(ji,jj,jk,jpphy,Krhs) = tr(ji,jj,jk,jpphy,Krhs)         &
            &                     + zprorcan(ji,jj,jk) * texcretn  &
            &                     - xpsino3 * zpronewn(ji,jj,jk)   &
-           &                     - xpsinh4 * zproregn(ji,jj,jk)   &
-           &                     - zrespn(ji,jj,jk) 
+           &                     - xpsinh4 * zproregn(ji,jj,jk)   
 
         tr(ji,jj,jk,jpnph,Krhs) = tr(ji,jj,jk,jpnph,Krhs) + zprontot * texcretn
         tr(ji,jj,jk,jppph,Krhs) = tr(ji,jj,jk,jppph,Krhs) + ( zpropo4n(ji,jj,jk) + zprodopn(ji,jj,jk) ) * texcretn
@@ -468,8 +467,7 @@ CONTAINS
         tr(ji,jj,jk,jppic,Krhs) = tr(ji,jj,jk,jppic,Krhs)         &
            &                     + zprorcap(ji,jj,jk) * texcretp  &
            &                     - xpsino3 * zpronewp(ji,jj,jk)   &
-           &                     - xpsinh4 * zproregp(ji,jj,jk)   &
-           &                     - zrespp(ji,jj,jk) 
+           &                     - xpsinh4 * zproregp(ji,jj,jk)   
 
         tr(ji,jj,jk,jpnpi,Krhs) = tr(ji,jj,jk,jpnpi,Krhs) + zproptot * texcretp
         tr(ji,jj,jk,jpppi,Krhs) = tr(ji,jj,jk,jpppi,Krhs) + ( zpropo4p(ji,jj,jk) + zprodopp(ji,jj,jk) ) * texcretp
@@ -479,8 +477,7 @@ CONTAINS
         tr(ji,jj,jk,jpdia,Krhs) = tr(ji,jj,jk,jpdia,Krhs)         &
            &                    + zprorcad(ji,jj,jk) * texcretd   &
            &                    - xpsino3 * zpronewd(ji,jj,jk)    &
-           &                    - xpsinh4 * zproregd(ji,jj,jk)    &
-           &                    - zrespd(ji,jj,jk) 
+           &                    - xpsinh4 * zproregd(ji,jj,jk)    
 
         tr(ji,jj,jk,jpndi,Krhs) = tr(ji,jj,jk,jpndi,Krhs) + zprodtot * texcretd
         tr(ji,jj,jk,jppdi,Krhs) = tr(ji,jj,jk,jppdi,Krhs) + ( zpropo4d(ji,jj,jk) + zprodopd(ji,jj,jk) ) * texcretd
@@ -491,7 +488,7 @@ CONTAINS
         tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + zproddop
   
         tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) &
-           &                     + o2ut * zpregtot + ( o2ut + o2nit ) * zpnewtot - o2ut * zresptot 
+           &                     + o2ut * zpregtot + ( o2ut + o2nit ) * zpnewtot
 
         tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - zprodfer
         consfe3(ji,jj,jk)       = zprodfer * 75.0 / ( rtrn + ( plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) )   &
@@ -524,14 +521,19 @@ CONTAINS
 
     ! Total primary production per year
     IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
-      & tpp = glob_sum( 'p5zprod', ( zprorcan(:,:,:) + zprorcad(:,:,:) + zprorcap(:,:,:) ) * cvol(:,:,:) )
+      & tpp = glob_sum( 'p5zprod', ( ( zprorcan(:,:,:) - xpsino3 * zpronewn(:,:,:) - xpsinh4 * zproregn(:,:,:)  &
+              &                      + zprorcad(:,:,:) - xpsino3 * zpronewd(:,:,:) - xpsinh4 * zproregd(:,:,:)  &
+              &                      + zprorcap(:,:,:) - xpsino3 * zpronewp(:,:,:) - xpsinh4 * zproregp(:,:,:) ) * cvol(:,:,:) ) )
 
     IF( lk_iomput .AND.  knt == nrdttrc ) THEN
        zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
        !
-       CALL iom_put( "PPPHYP"  , zprorcap(:,:,:) * zfact * tmask(:,:,:)   ) ! primary production by picophyto
-       CALL iom_put( "PPPHYN"  , zprorcan(:,:,:) * zfact * tmask(:,:,:) )  ! primary production by nanophyto
-       CALL iom_put( "PPPHYD"  , zprorcad(:,:,:) * zfact * tmask(:,:,:)   ) ! primary production by diatomes
+       CALL iom_put( "GPPHYP"  , zprorcap(:,:,:) * zfact * tmask(:,:,:) ) ! primary production by picophyto
+       CALL iom_put( "GPPHYN"  , zprorcan(:,:,:) * zfact * tmask(:,:,:) )  ! primary production by nanophyto
+       CALL iom_put( "GPPHYD"  , zprorcad(:,:,:) * zfact * tmask(:,:,:) ) ! primary production by diatomes
+       CALL iom_put( "PPPHYP"  , ( zprorcap(:,:,:) - xpsino3 * zpronewp(:,:,:) - xpsinh4 * zproregp(:,:,:) ) * zfact * tmask(:,:,:) ) ! primary production by picophyto
+       CALL iom_put( "PPPHYN"  , ( zprorcan(:,:,:) - xpsino3 * zpronewn(:,:,:) - xpsinh4 * zproregn(:,:,:) ) * zfact * tmask(:,:,:) )  ! primary production by nanophyto
+       CALL iom_put( "PPPHYD"  , ( zprorcad(:,:,:) - xpsino3 * zpronewd(:,:,:) - xpsinh4 * zproregd(:,:,:) ) * zfact * tmask(:,:,:) ) ! primary production by diatomes
        CALL iom_put( "PPNEWN"  , zpronewp(:,:,:) * zfact * tmask(:,:,:)    ) ! new primary production by picophyto
        CALL iom_put( "PPNEWN"  , zpronewn(:,:,:) * zfact * tmask(:,:,:)    ) ! new primary production by nanophyto
        CALL iom_put( "PPNEWD"  , zpronewd(:,:,:) * zfact * tmask(:,:,:)   ) ! new primary production by diatomes
@@ -590,7 +592,7 @@ CONTAINS
       INTEGER :: ios    ! Local integer output status for namelist read
       !!
       NAMELIST/namp5zprod/ pislopen, pislopep, pisloped, excretn, excretp, excretd,     &
-         &                 thetannm, thetanpm, thetandm, chlcmin, grosip, bresp, xadap
+         &                 thetannm, thetanpm, thetandm, chlcmin, grosip, bresp
       !!----------------------------------------------------------------------
 
       READ  ( numnatp_ref, namp5zprod, IOSTAT = ios, ERR = 901)
@@ -608,7 +610,6 @@ CONTAINS
          WRITE(numout,*) '    P-I slope                                 pislopen     =', pislopen
          WRITE(numout,*) '    P-I slope  for diatoms                    pisloped     =', pisloped
          WRITE(numout,*) '    P-I slope  for picophytoplankton          pislopep     =', pislopep
-         WRITE(numout,*) '    Acclimation factor to low light           xadap        =', xadap
          WRITE(numout,*) '    excretion ratio of nanophytoplankton      excretn      =', excretn
          WRITE(numout,*) '    excretion ratio of picophytoplankton      excretp      =', excretp
          WRITE(numout,*) '    excretion ratio of diatoms                excretd      =', excretd