diff --git a/src/OCE/trc_oce.F90 b/src/OCE/trc_oce.F90
index 6810cb875c9a6faf919f653d4e94082114a8617a..9f1bb10e85c338f584f512cd26875cbcaa8f5b3f 100644
--- a/src/OCE/trc_oce.F90
+++ b/src/OCE/trc_oce.F90
@@ -43,6 +43,9 @@ MODULE trc_oce
    !!----------------------------------------------------------------------
    LOGICAL, PUBLIC, PARAMETER ::   lk_top     = .FALSE.   !: TOP model
 #endif
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
    !! $Id: trc_oce.F90 13286 2020-07-09 15:48:29Z smasson $ 
diff --git a/src/TOP/PISCES/P2Z/p2zbio.F90 b/src/TOP/PISCES/P2Z/p2zbio.F90
index babf325b7bef50d389d5d63aa87cd471fe5efb31..4100534f899bafbc79999a5531bc38bcd0f3594b 100644
--- a/src/TOP/PISCES/P2Z/p2zbio.F90
+++ b/src/TOP/PISCES/P2Z/p2zbio.F90
@@ -113,7 +113,7 @@ CONTAINS
 
       IF( lk_iomput ) THEN
          ALLOCATE( zw3d(A2D(0),jpk,3) )   ;   zw3d(:,:,jpk,:) = 0._wp
-         ALLOCATE( zw2d(A2D(0),17) )
+         ALLOCATE( zw2d(A2D(0),17) )      ;   zw2d(:,:,:) = 0._wp
       ENDIF
       !
       xksi(:,:) = 0.e0        ! zooplakton closure ( fbod)
@@ -366,8 +366,6 @@ CONTAINS
          CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
-      IF( lk_iomput )   DEALLOCATE( zw2d, zw3d )
-      !
       IF( ln_timing )  CALL timing_stop('p2z_bio')
       !
    END SUBROUTINE p2z_bio
diff --git a/src/TOP/PISCES/P4Z/p4zpoc.F90 b/src/TOP/PISCES/P4Z/p4zpoc.F90
index 37347fce0c8b851ffb541f9c2913b87fb4b2af96..823d416a856829fef852d7616d780b6d50aeaf6c 100644
--- a/src/TOP/PISCES/P4Z/p4zpoc.F90
+++ b/src/TOP/PISCES/P4Z/p4zpoc.F90
@@ -111,7 +111,7 @@ CONTAINS
      ! lability class is specified in the namelist, this is equivalent to 
      ! a standard parameterisation with a constant lability
      ! -----------------------------------------------------------------------
-     DO_3D( 0, 0, 0, 0, 1, jpkm1)
+     DO_3D( 0, 0, 0, 0, 2, jpkm1)
         IF (tmask(ji,jj,jk) == 1.) THEN
           zdep = hmld(ji,jj)
           !