diff --git a/src/TOP/PISCES/P4Z/p4zlys.F90 b/src/TOP/PISCES/P4Z/p4zlys.F90
index 007835e179b9aed11fd54f31bedf6e3bc42d167e..cc0f80ff519095c8e27752d39de1e53ed3c00b73 100644
--- a/src/TOP/PISCES/P4Z/p4zlys.F90
+++ b/src/TOP/PISCES/P4Z/p4zlys.F90
@@ -64,41 +64,43 @@ CONTAINS
       !
       INTEGER  ::   ji, jj, jk, jn
       REAL(wp) ::   zdispot, zrhd, zcalcon
-      REAL(wp) ::   zomegaca, zexcess, zexcess0, zkd
+      REAL(wp) ::   zomegaca, zexcess, zexcess0, zkd, zco3, ztra
       CHARACTER (len=25) ::   charout
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zco3, zcaldiss, zhinit, zhi, zco3sat
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zhinit, zhi
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zcaldiss, zbicarb, zco3sat
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )  CALL timing_start('p4z_lys')
       !
+      IF( iom_use( "CO3"    ) ) ALLOCATE( zbicarb (A2D(0),jpk) )  
+      IF( iom_use( "CO3sat" ) ) ALLOCATE( zco3sat (A2D(0),jpk) ) 
+      IF( iom_use( "DCAL"   ) ) ALLOCATE( zcaldiss(A2D(0),jpk) ) 
 
-      zhinit  (:,:,:) = hi(:,:,:) / ( rhd(:,:,:) + 1._wp )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+        zhinit(ji,jj,jk) = hi(ji,jj,jk) / ( rhd(ji,jj,jk) + 1._wp )
+      END_3D
       !
       !     -------------------------------------------
       !     COMPUTE [CO3--] and [H+] CONCENTRATIONS
       !     -------------------------------------------
 
       CALL solve_at_general( zhinit, zhi, Kbb )
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-         zco3(ji,jj,jk) = tr(ji,jj,jk,jpdic,Kbb) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   &
+
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+         !
+         zco3           = tr(ji,jj,jk,jpdic,Kbb) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   &
             &             + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn )
          hi  (ji,jj,jk) = zhi(ji,jj,jk) * ( rhd(ji,jj,jk) + 1._wp )
-      END_3D
 
-      !     ---------------------------------------------------------
-      !        CALCULATE DEGREE OF CACO3 SATURATION AND CORRESPONDING
-      !        DISSOLOUTION AND PRECIPITATION OF CACO3 (BE AWARE OF
-      !        MGCO3)
-      !     ---------------------------------------------------------
+         ! CALCULATE DEGREE OF CACO3 SATURATION AND CORRESPONDING
+         ! DISSOLOUTION AND PRECIPITATION OF CACO3 (BE AWARE OF MGCO3)
 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
 
          ! DEVIATION OF [CO3--] FROM SATURATION VALUE
          ! Salinity dependance in zomegaca and divide by rhd to have good units
          zcalcon  = calcon * ( salinprac(ji,jj,jk) / 35._wp )
          zrhd    = rhd(ji,jj,jk) + 1._wp
-         zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zrhd + rtrn )
-         zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zrhd / ( zcalcon + rtrn )
+         zomegaca = ( zcalcon * zco3 ) / ( aksp(ji,jj,jk) * zrhd + rtrn )
 
          ! SET DEGREE OF UNDER-/SUPERSATURATION
          excess(ji,jj,jk) = 1._wp - zomegaca
@@ -116,24 +118,35 @@ CONTAINS
 
         !  CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
         !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
-        zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss ! calcite dissolution
+        ztra =  zdispot * rfact2 / rmtss ! calcite dissolution
+        !
+        tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + 2. * ztra
+        tr(ji,jj,jk,jpcal,Krhs) = tr(ji,jj,jk,jpcal,Krhs) -      ztra
+        tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) +      ztra
+        ! 
+        IF( iom_use( "CO3"    ) ) zbicarb (ji,jj,jk) = zco3                       !  bicarbonate
+        IF( iom_use( "CO3sat" ) ) zco3sat (ji,jj,jk) = zco3 / ( zomegaca + rtrn )  ! calcite saturation
+        IF( iom_use( "DCAL"   ) ) zcaldiss(ji,jj,jk) = ztra                        ! calcite dissolution
         !
-        tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + 2. * zcaldiss(ji,jj,jk)
-        tr(ji,jj,jk,jpcal,Krhs) = tr(ji,jj,jk,jpcal,Krhs) -      zcaldiss(ji,jj,jk)
-        tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) +      zcaldiss(ji,jj,jk)
       END_3D
       !
 
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-         CALL iom_put( "PH" , -1. * LOG10( MAX( hi(:,:,:), rtrn ) ) * tmask(:,:,:) )
-         IF( iom_use( "CO3" ) ) THEN
-            zco3(:,:,jpk) = 0.    ; CALL iom_put( "CO3"   , zco3(:,:,:)     * 1.e+3           * tmask(:,:,:) )
+         IF( iom_use( "PH" ) )  CALL iom_put( "PH" , -1. * LOG10( MAX( hi(A2D(0),1:jpk), rtrn ) ) * tmask(A2D(0),1:jpk) )
+         IF( iom_use( "CO3"    ) ) THEN  ! bicarbonate 
+            zbicarb(A2D(0),jpk)  = 0. 
+            CALL iom_put( "CO3"   , zbicarb(A2D(0),1:jpk) * 1.e+3 * tmask(A2D(0),1:jpk) ) 
+            DEALLOCATE( zbicarb )
          ENDIF
-         IF( iom_use( "CO3sat" ) ) THEN
-           zco3sat(:,:,jpk) = 0.  ; CALL iom_put( "CO3sat", zco3sat(:,:,:)  * 1.e+3           * tmask(:,:,:) )
+         IF( iom_use( "CO3sat" ) ) THEN   ! calcite saturation
+            zco3sat(A2D(0),jpk)  = 0. 
+            CALL iom_put( "CO3sat", zco3sat(A2D(0),1:jpk) * 1.e+3 * tmask(A2D(0),1:jpk) ) 
+            DEALLOCATE( zco3sat )
          ENDIF
-         IF( iom_use( "DCAL" ) ) THEN
-           zcaldiss(:,:,jpk) = 0. ; CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) )
+         IF( iom_use( "DCAL" ) ) THEN  ! calcite dissolution
+            zcaldiss(A2D(0),jpk)  = 0. 
+            CALL iom_put( "DCAL", zcaldiss(A2D(0),1:jpk) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpk) ) 
+            DEALLOCATE( zcaldiss )
          ENDIF              
       ENDIF
       !
diff --git a/src/TOP/PISCES/P4Z/p4zprod.F90 b/src/TOP/PISCES/P4Z/p4zprod.F90
index 59ca90ccc62b5933ddad31ffff1460965cfc528d..357611278a1c4d7e49cead02c90f1a584c39989a 100644
--- a/src/TOP/PISCES/P4Z/p4zprod.F90
+++ b/src/TOP/PISCES/P4Z/p4zprod.F90
@@ -75,73 +75,64 @@ CONTAINS
       REAL(wp) ::   zproddoc, zprodsil, zprodfer, zprodlig
       REAL(wp) ::   zpislopen, zpisloped, zfact
       REAL(wp) ::   zratiosi, zmaxsi, zlimfac, zsizetmp, zfecnm, zfecdm
-      REAL(wp) ::   zprod, zval
+      REAL(wp) ::   zprod, zval, zmxl_fac, zmxl_chl
+      REAL(wp) ::   zprmax, zpislopeadn,zpislopeadd,zprchld, zprchln   
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprmaxn,zprmaxd
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadd, zysopt  
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprdia, zprbio, zprchld, zprchln   
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprorcan, zprorcad, zprofed, zprofen
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpronewn, zpronewd
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmxl_fac, zmxl_chl
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zpligprod
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zprdia, zprbio, zysopt, zmxl    
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zprorcan, zprorcad, zprofed, zprofen
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zpronewn, zpronewd
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_prod')
       !
-      !  Allocate temporary workspace
-      !
-      zprorcan(:,:,:) = 0._wp ; zprorcad(:,:,:) = 0._wp ; zprofed(:,:,:) = 0._wp
-      zprofen (:,:,:) = 0._wp ; zysopt  (:,:,:) = 0._wp
-      zpronewn(:,:,:) = 0._wp ; zpronewd(:,:,:) = 0._wp ; zprdia(:,:,:) = 0._wp
-      zprbio  (:,:,:) = 0._wp ; zprchld (:,:,:) = 0._wp ; zprchln(:,:,:) = 0._wp 
-      zmxl_fac(:,:,:) = 0._wp ; zmxl_chl(:,:,:) = 0._wp 
-
-      ! Computation of the maximimum production. Based on a Q10 description
-      ! of the thermal dependency
-      ! Parameters are taken from Bissinger et al. (2008)
-      zprmaxn(:,:,:) = 0.65_wp * r1_rday * tgfunc(:,:,:)
-      zprmaxd(:,:,:) = zprmaxn(:,:,:)
 
       ! Intermittency is supposed to have a similar effect on production as 
-      ! day length (Shatwell et al., 2012). The correcting factor is zmxl_fac. 
-      ! zmxl_chl is the fractional day length and is used to compute the mean
-      ! PAR during daytime. The effect of mixing is computed using the 
-      ! absolute light level definition of the euphotic zone
+      ! day length (Shatwell et al., 2012).  
       ! ------------------------------------------------------------------------- 
       IF ( ln_p4z_dcyc ) THEN    ! Diurnal cycle in PISCES
-
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
             IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
                   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 )
+               zmxl(ji,jj,jk) = zval
             ENDIF
          END_3D
  
       ELSE ! No diurnal cycle in PISCES
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
             IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                zval = MAX( 1., strn(ji,jj) )
                IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
                   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 )
+               zmxl(ji,jj,jk) = zval
             ENDIF
          END_3D
 
       ENDIF
 
-      zprbio(:,:,:) = zprmaxn(:,:,:) * zmxl_fac(:,:,:)
-      zprdia(:,:,:) = zprmaxd(:,:,:) * zmxl_fac(:,:,:)
-
       ! The formulation proposed by Geider et al. (1997) has been modified 
       ! to exclude the effect of nutrient limitation and temperature in the PI
       ! curve following Vichi et al. (2007)
       ! -----------------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+         !
+         ! Computation of the maximimum production. Based on a Q10 description
+         ! of the thermal dependency. Parameters are taken from Bissinger et al. (2008)
+         zprmax = 0.65_wp * r1_rday * tgfunc(ji,jj,jk)
+         !
+         ! The correcting factor of Intermittency is zmxl_fac. 
+         ! zmxl_chl is the fractional day length and is used to compute the mean
+         ! PAR during daytime. The effect of mixing is computed using the 
+         ! absolute light level definition of the euphotic zone
+         zmxl_fac = 1.0 - EXP( -0.26 * zmxl(ji,jj,jk) )
+         zmxl_chl = zmxl(ji,jj,jk) / 24.
+         !
+         zprbio(ji,jj,jk) = zprmax * zmxl_fac 
+         zprdia(ji,jj,jk) = zprmax * zmxl_fac
+         !
          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 )
@@ -154,66 +145,59 @@ CONTAINS
             ! improved or removed in future versions of the model
 
             ! Nanophytoplankton
-            zpislopeadn(ji,jj,jk) = pislopen * ( 1.+ zadap  * EXP( -0.25 * enano(ji,jj,jk) ) )  &
+            zpislopeadn = pislopen * ( 1.+ zadap  * EXP( -0.25 * enano(ji,jj,jk) ) )  &
             &                   * 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 )   &
+            zpislopeadd = (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)
+            !
+            ! Computation of production function for Carbon ; Actual light levels are used here 
+            zfact = 1._wp / ( ( r1_rday + bresp * r1_rday ) * zmxl_fac * rday + rtrn)
+            zpislopen = zpislopeadn * zfact
+            zpisloped = zpislopeadd * zfact
+            !
+            zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
+            zprdia(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
          ENDIF
-      END_3D
+         !  Computation of a proxy of the N/C quota from nutrient limitation 
+         !  and light limitation. Steady state is assumed to allow the computation
+         !  ----------------------------------------------------------------------
+         zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) )   &
+         &      * zprmax / ( zprbio(ji,jj,jk) + rtrn )
+         quotan(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
+         !
+         zval = MIN( xdiatpo4(ji,jj,jk), ( xdiatnh4(ji,jj,jk) + xdiatno3(ji,jj,jk) ) )   &
+         &      * zprmax / ( zprdia(ji,jj,jk) + rtrn )
+         quotad(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
+
+         ! Sea-ice effect on production : No production is assumed below sea ice
+         zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
+         zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-             ! Computation of production function for Carbon
-             ! Actual light levels are used here 
-             ! ----------------------------------------------
-             zpislopen = zpislopeadn(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
-             &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
-             zpisloped = zpislopeadd(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
-             &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
-             zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(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)
-             !  is used here (acclimation is in general slower than 
-             !  the characteristic time scales of vertical mixing)
-             !  ------------------------------------------------------
-             zpislopen = zpislopeadn(ji,jj,jk) / ( zprmaxn(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
-             zpisloped = zpislopeadd(ji,jj,jk) / ( zprmaxd(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
-             zprchln(ji,jj,jk) = zprmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) ) )
-             zprchld(ji,jj,jk) = zprmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) ) )
-         ENDIF
-      END_3D
-
-      !  Computation of a proxy of the N/C quota from nutrient limitation 
-      !  and light limitation. Steady state is assumed to allow the computation
-      !  ----------------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-          zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) )   &
-          &      * zprmaxn(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn )
-          quotan(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
-          zval = MIN( xdiatpo4(ji,jj,jk), ( xdiatnh4(ji,jj,jk) + xdiatno3(ji,jj,jk) ) )   &
-          &      * zprmaxd(ji,jj,jk) / ( zprdia(ji,jj,jk) + rtrn )
-          quotad(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
-      END_3D
-
-
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-
-          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-             ! Si/C of diatoms
-             ! ------------------------
-             ! 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 (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).
-             ! -----------------------------------------------------------------------
+            !  Computation of production function for Chlorophyll
+            !  Mean light level in the mixed layer (when appropriate)
+            !  is used here (acclimation is in general slower than 
+            !  the characteristic time scales of vertical mixing)
+            !  ------------------------------------------------------
+            zfact = 1._wp / ( zprmax * zmxl_chl * rday + rtrn )
+            zpislopen = zpislopeadn * zfact 
+            zpisloped = zpislopeadd * zfact
+            zprchln = zprmax * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) ) )
+            zprchld = zprmax * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) ) )
+            ! Si/C of diatoms
+            ! ------------------------
+            ! 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 (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).
+            ! -----------------------------------------------------------------------
             zlim  = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi1 )
-            zsilim = xlimdia(ji,jj,jk) * zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn )
+            !
+            zsilim = xlimdia(ji,jj,jk) * zprdia(ji,jj,jk) / ( zprmax + 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
               zsilfac = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
@@ -228,21 +212,9 @@ CONTAINS
             ELSE
                zysopt(ji,jj,jk) = zlim * zsilfac * grosip * 1.0 * zsilim**0.7 * zmaxsi
             ENDIF
-        ENDIF
-      END_3D
-
-      ! 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) )
-         zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
-      END_3D
 
-      ! Computation of the various production  and nutrient uptake terms
-      ! ---------------------------------------------------------------
-      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 various production  and nutrient uptake terms
+            ! ---------------------------------------------------------------
             !  production terms for nanophyto. (C)
             zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * tr(ji,jj,jk,jpphy,Kbb) * rfact2
 
@@ -255,7 +227,7 @@ CONTAINS
             ! size at time step t+1 and is thus updated at the end of the 
             ! current time step
             ! --------------------------------------------------------------------
-            zlimfac = xlimphy(ji,jj,jk) * zprchln(ji,jj,jk) / ( zprmaxn(ji,jj,jk) + rtrn )
+            zlimfac = xlimphy(ji,jj,jk) * zprchln / ( zprmax + 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 ) )
 
@@ -266,7 +238,7 @@ CONTAINS
             zfecnm = xqfuncfecn(ji,jj,jk) + ( fecnm - xqfuncfecn(ji,jj,jk) ) * ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) )
             zratio = 1.0 - MIN(1.0,tr(ji,jj,jk,jpnfe,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) * zfecnm + rtrn ) )
             zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) ) 
-            zprofen(ji,jj,jk) = zfecnm * zprmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
+            zprofen(ji,jj,jk) = zfecnm * zprmax * ( 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) ) )   &
             &          * xnanofer(ji,jj,jk) * zmax * tr(ji,jj,jk,jpphy,Kbb) * rfact2
@@ -282,7 +254,7 @@ CONTAINS
             ! size at time step t+1 and is thus updated at the end of the 
             ! current time step. 
             ! --------------------------------------------------------------------
-            zlimfac = zprchld(ji,jj,jk) * xlimdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn )
+            zlimfac = zprchld * xlimdia(ji,jj,jk) / ( zprmax + rtrn )
             zsizetmp = 1.0 + 1.3 * ( xsizerd - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3)
             sizeda(ji,jj,jk) = min(xsizerd, max( sizeda(ji,jj,jk), zsizetmp ) )
 
@@ -293,45 +265,37 @@ CONTAINS
             zfecdm = xqfuncfecd(ji,jj,jk) + ( fecdm - xqfuncfecd(ji,jj,jk) ) * ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) )
             zratio = 1.0 - MIN(1.0, tr(ji,jj,jk,jpdfe,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) * zfecdm + rtrn ) )
             zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) ) 
-            zprofed(ji,jj,jk) = zfecdm * zprmaxd(ji,jj,jk) * (1.0 - fr_i(ji,jj) )  &
+            zprofed(ji,jj,jk) = zfecdm * zprmax * (1.0 - fr_i(ji,jj) )  &
             &          * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn + xdiatno3(ji,jj,jk)  &
             &          + xdiatnh4(ji,jj,jk) ) * (1. - xdiatfer(ji,jj,jk) ) )   &
             &          * xdiatfer(ji,jj,jk) * zmax * tr(ji,jj,jk,jpdia,Kbb) * rfact2
-         ENDIF
-      END_3D
 
-      ! Computation of the chlorophyll production terms
-      ! The parameterization is taken from Geider et al. (1997)
-      ! -------------------------------------------------------
-      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 chlorophyll production terms
+            ! The parameterization is taken from Geider et al. (1997)
+
             !  production terms for nanophyto. ( chlorophyll )
-            znanotot = enanom(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
-            zprod    = rday * zprorcan(ji,jj,jk) * zprchln(ji,jj,jk) * xlimphy(ji,jj,jk)
+            znanotot = enanom(ji,jj,jk) / ( zmxl_chl + rtrn )
+            zprod    = rday * zprorcan(ji,jj,jk) * zprchln * xlimphy(ji,jj,jk)
             zprochln = chlcmin * 12. * zprorcan (ji,jj,jk)
             zprochln = zprochln + (chlcnm - chlcmin) * 12. * zprod / &
-                                  & (  zpislopeadn(ji,jj,jk) * znanotot +rtrn)
+                                  & (  zpislopeadn * znanotot +rtrn)
 
             !  production terms for diatoms ( chlorophyll )
-            zdiattot = ediatm(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
-            zprod    = rday * zprorcad(ji,jj,jk) * zprchld(ji,jj,jk) * xlimdia(ji,jj,jk)
+            zdiattot = ediatm(ji,jj,jk) / ( zmxl_chl + rtrn )
+            zprod    = rday * zprorcad(ji,jj,jk) * zprchld * xlimdia(ji,jj,jk)
             zprochld = chlcmin * 12. * zprorcad(ji,jj,jk)
             zprochld = zprochld + (chlcdm - chlcmin) * 12. * zprod / &
-                                  & ( zpislopeadd(ji,jj,jk) * zdiattot +rtrn )
+                                  & ( zpislopeadd * zdiattot +rtrn )
 
             !   Update the arrays TRA which contain the Chla sources and sinks
             tr(ji,jj,jk,jpnch,Krhs) = tr(ji,jj,jk,jpnch,Krhs) + zprochln * texcretn
             tr(ji,jj,jk,jpdch,Krhs) = tr(ji,jj,jk,jpdch,Krhs) + zprochld * texcretd
-         ENDIF
-      END_3D
 
-      !   Update the arrays TRA which contain the biological sources and sinks
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
-        IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+            !   Update the arrays TRA which contain the biological sources and sinks
            zpptot   = zprorcan(ji,jj,jk) + zprorcad(ji,jj,jk)
            zpnewtot = zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk)
            zpregtot = zpptot - zpnewtot
-           zprodsil  = zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk) * rfact2 * tr(ji,jj,jk,jpdia,Kbb)
+           zprodsil  = zprmax * zysopt(ji,jj,jk) * rfact2 * tr(ji,jj,jk,jpdia,Kbb)
            zproddoc  = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)
            zprodfer  = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)
            !
@@ -353,7 +317,7 @@ CONTAINS
            !
            tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) - zpptot
            tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * ( zpnewtot - zpregtot )
-        ENDIF
+         ENDIF
       END_3D
 
      ! Production and uptake of ligands by phytoplankton. This part is activated 
@@ -362,7 +326,7 @@ CONTAINS
      ! Shaked et al. (2020)
      ! -------------------------------------------------------------------------
      IF( ln_ligand ) THEN
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
            IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
               zproddoc = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)
               zprodfer = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)
@@ -373,41 +337,109 @@ CONTAINS
          END_3D
      ENDIF
 
-
     ! Output of the diagnostics
     ! Total primary production per year
-    IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
-         & tpp = glob_sum( 'p4zprod', ( zprorcan(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
+    IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  THEN
+       ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp  
+       zw3d(A2D(0),1:jpkm1) = ( zprorcan(A2D(0),1:jpkm1) + zprorcad(A2D(0),1:jpkm1) ) * cvol(A2D(0),1:jpkm1) 
+       tpp = glob_sum( 'p4zprod', zw3d )
+       DEALLOCATE ( zw3d ) 
+    ENDIF
 
     IF( lk_iomput .AND.  knt == nrdttrc ) THEN
        zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
        !
-       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( "PPNEWN"  , zpronewn(:,:,:) * zfact * tmask(:,:,:)    ) ! new primary production by nanophyto
-       CALL iom_put( "PPNEWD"  , zpronewd(:,:,:) * zfact * tmask(:,:,:)   ) ! new primary production by diatomes
-       CALL iom_put( "PBSi"    , zprorcad(:,:,:) * zfact * tmask(:,:,:) * zysopt(:,:,:)  ) ! biogenic silica production
-       CALL iom_put( "PFeN"    , zprofen(:,:,:)  * zfact * tmask(:,:,:)  ) ! biogenic iron production by nanophyto
-       CALL iom_put( "PFeD"    , zprofed(:,:,:)  * zfact * tmask(:,:,:)  ) ! biogenic iron production by  diatomes
+       ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp  
+       !
+       IF( iom_use ( "PPPHYN" ) ) THEN   ! primary production by nanophyto
+         zw3d(A2D(0),1:jpkm1) = zprorcan(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "PPPHYN", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "PPPHYD" ) ) THEN    ! primary production by diatomes
+         zw3d(A2D(0),1:jpkm1) = zprorcad(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "PPPHYD", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "PPNEWN" ) ) THEN    ! new primary production by nano
+         zw3d(A2D(0),1:jpkm1) = zpronewn(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "PPNEWN", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "PPNEWD" ) ) THEN    ! new primary production by diatomes
+         zw3d(A2D(0),1:jpkm1) = zpronewd(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "PPNEWD", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "PBSi" ) ) THEN    !  biogenic silica production
+         zw3d(A2D(0),1:jpkm1) = zprorcad(A2D(0),1:jpkm1) * zysopt(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "PBSi", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "PFeN" ) ) THEN    ! biogenic iron production by nanophyto
+         zw3d(A2D(0),1:jpkm1) = zprofen(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "PFeN", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "PFeD" ) ) THEN    ! biogenic iron production by nanophyto
+         zw3d(A2D(0),1:jpkm1) = zprofed(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "PFeD", zw3d )  
+       ENDIF
+       !
        IF( ln_ligand .AND. ( iom_use( "LPRODP" ) .OR. iom_use( "LDETP" ) ) ) THEN
-           ALLOCATE(  zpligprod(jpi,jpj,jpk) )
-           zpligprod(:,:,:) = excretd * zprorcad(:,:,:) + excretn * zprorcan(:,:,:)
-           CALL iom_put( "LPRODP"  , zpligprod(:,:,:) * ldocp * 1e9 * zfact * tmask(:,:,:) )
+           zw3d(A2D(0),1:jpkm1) = excretd * zprorcad(A2D(0),1:jpkm1) + excretn * zprorcan(A2D(0),1:jpkm1) &
+             &                   * zfact * tmask(A2D(0),1:jpkm1) 
+           CALL iom_put( "LPRODP"  , zw3d * ldocp * 1e9 )
            !
-           zpligprod(:,:,:) = ( texcretn * zprofen(:,:,:) + texcretd * zprofed(:,:,:) ) & 
-             &                  * plig(:,:,:) / ( rtrn + plig(:,:,:) + 75.0 * (1.0 - plig(:,:,:) ) )
-           CALL iom_put( "LDETP"   , zpligprod(:,:,:)  * lthet * 1e9 * zfact * tmask(:,:,:) )
-           DEALLOCATE(  zpligprod )
+           zw3d(A2D(0),1:jpkm1) = ( texcretn * zprofen(A2D(0),1:jpkm1) + texcretd * zprofed(A2D(0),1:jpkm1) ) & 
+             &                  * plig(A2D(0),1:jpkm1) / ( rtrn + plig(A2D(0),1:jpkm1) + 75.0 * (1.0 - plig(A2D(0),1:jpkm1) ) )  &
+             &                  * zfact * tmask(A2D(0),1:jpkm1) 
+           CALL iom_put( "LDETP"   , zw3d * lthet * 1e9 )
+       ENDIF
+       !
+       IF( iom_use ( "Mumax" ) ) THEN    ! Maximum growth rate
+         zw3d(A2D(0),1:jpkm1) = 0.65_wp * r1_rday * tgfunc(A2D(0),1:jpkm1)  * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "Mumax", zw3d )  
        ENDIF
-       CALL iom_put( "Mumax"   , zprmaxn(:,:,:) * tmask(:,:,:)  ) ! Maximum growth rate
-       CALL iom_put( "MuN"     , zprbio(:,:,:) * xlimphy(:,:,:) * tmask(:,:,:) ) ! Realized growth rate for nanophyto
-       CALL iom_put( "MuD"     , zprdia(:,:,:) * xlimdia(:,:,:) * tmask(:,:,:) ) ! Realized growth rate for diatoms
-       CALL iom_put( "LNlight" , zprbio (:,:,:) / (zprmaxn(:,:,:) + rtrn) * tmask(:,:,:)  )  ! light limitation term
-       CALL iom_put( "LDlight" , zprdia (:,:,:) / (zprmaxd(:,:,:) + rtrn) * tmask(:,:,:)   )
-       CALL iom_put( "TPP"     , ( zprorcan(:,:,:) + zprorcad(:,:,:) ) * zfact * tmask(:,:,:)  )  ! total primary production
-       CALL iom_put( "TPNEW"   , ( zpronewn(:,:,:) + zpronewd(:,:,:) ) * zfact * tmask(:,:,:)  ) ! total new production
-       CALL iom_put( "TPBFE"   , ( zprofen(:,:,:) + zprofed(:,:,:) ) * zfact * tmask(:,:,:)  )  ! total biogenic iron production
-       CALL iom_put( "tintpp"  , tpp * zfact )  !  global total integrated primary production molC/s
+       !
+       IF( iom_use ( "MuN" ) ) THEN    ! Realized growth rate for nanophyto
+         zw3d(A2D(0),1:jpkm1) = zprbio(A2D(0),1:jpkm1) * xlimphy(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "MuN", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "MuD" ) ) THEN    ! Realized growth rate for diatoms
+         zw3d(A2D(0),1:jpkm1) = zprdia(A2D(0),1:jpkm1) * xlimdia(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "MuD", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "LNlight" ) ) THEN    ! light limitation term for nano
+         zw3d(A2D(0),1:jpkm1) = zprbio(A2D(0),1:jpkm1) / ( 0.65_wp * r1_rday * tgfunc(A2D(0),1:jpkm1) + rtrn ) * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "LNlight", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "LDlight" ) ) THEN    ! light limitation term for diatomes
+         zw3d(A2D(0),1:jpkm1) = zprdia(A2D(0),1:jpkm1) / ( 0.65_wp * r1_rday * tgfunc(A2D(0),1:jpkm1) + rtrn ) * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "LDlight", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "TPP" ) ) THEN   ! total primary production
+         zw3d(A2D(0),1:jpkm1) = ( zprorcan(A2D(0),1:jpkm1) +  zprorcad(A2D(0),1:jpkm1) ) * zfact * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "TPP", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "TPNEW" ) ) THEN   ! total new production
+         zw3d(A2D(0),1:jpkm1) = ( zpronewn(A2D(0),1:jpkm1) +  zpronewd(A2D(0),1:jpkm1) ) * zfact * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "TPNEW", zw3d )  
+       ENDIF
+       !
+       IF( iom_use ( "TPBFE" ) ) THEN   ! total biogenic iron production
+         zw3d(A2D(0),1:jpkm1) = ( zprofen(A2D(0),1:jpkm1) +  zprofed(A2D(0),1:jpkm1) ) * zfact * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "TPBFE", zw3d )  
+       ENDIF
+       !
+       IF( iom_use( "tintpp") )  CALL iom_put( "tintpp"  , tpp * zfact )  !  global total integrated primary production molC/s
+       !
+       DEALLOCATE( zw3d )
      ENDIF
 
      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
@@ -480,7 +512,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_prod_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( quotan(jpi,jpj,jpk), quotad(jpi,jpj,jpk), STAT = p4z_prod_alloc )
+      ALLOCATE( quotan(A2D(0),jpk), quotad(A2D(0),jpk), STAT = p4z_prod_alloc )
       !
       IF( p4z_prod_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_prod_alloc : failed to allocate arrays.' )
       !
diff --git a/src/TOP/PISCES/sms_pisces.F90 b/src/TOP/PISCES/sms_pisces.F90
index ba890a006293caee267e405e6f29a9ed15445d83..58679aaeb53ac924ff2f50e78fe31b2fc4e1fd4d 100644
--- a/src/TOP/PISCES/sms_pisces.F90
+++ b/src/TOP/PISCES/sms_pisces.F90
@@ -124,6 +124,8 @@ MODULE sms_pisces
 
    LOGICAL, SAVE :: lk_sed
 
+      !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
    !! $Id: sms_pisces.F90 15459 2021-10-29 08:19:18Z cetlod $ 
@@ -145,47 +147,47 @@ CONTAINS
       IF( ln_p4z .OR. ln_p5z ) THEN
 
          !* Optics
-         ALLOCATE(  enano(jpi,jpj,jpk) , ediat(jpi,jpj,jpk) ,   &
-           &        enanom(jpi,jpj,jpk), ediatm(jpi,jpj,jpk),   &
-           &        emoy(jpi,jpj,jpk)  , etotm(jpi,jpj,jpk),   STAT=ierr(2) )
+         ALLOCATE(  enano(A2D(0),jpk) , ediat(A2D(0),jpk) ,   &
+           &        enanom(A2D(0),jpk), ediatm(A2D(0),jpk),   &
+           &        emoy(A2D(0),jpk)  , etotm(A2D(0),jpk),   STAT=ierr(2) )
 
          !* Biological SMS
-         ALLOCATE( xksimax(jpi,jpj)  , biron(jpi,jpj,jpk)      ,  STAT=ierr(3) )
+         ALLOCATE( xksimax(jpi,jpj)  , biron(A2D(0),jpk)      ,  STAT=ierr(3) )
 
          ! Biological SMS
-         ALLOCATE( xfracal  (jpi,jpj,jpk), orem    (jpi,jpj,jpk),  &
-            &      nitrfac  (jpi,jpj,jpk), nitrfac2(jpi,jpj,jpk),  &
-            &      prodcal  (jpi,jpj,jpk), xdiss   (jpi,jpj,jpk),  &
-            &      prodpoc  (jpi,jpj,jpk), conspoc (jpi,jpj,jpk),  &
-            &      prodgoc  (jpi,jpj,jpk), consgoc (jpi,jpj,jpk),  &
-            &      blim     (jpi,jpj,jpk), consfe3 (jpi,jpj,jpk),  &
-            &      xfecolagg(jpi,jpj,jpk), xcoagfe (jpi,jpj,jpk), STAT=ierr(4) )
+         ALLOCATE( xfracal  (A2D(0),jpk), orem    (A2D(0),jpk),  &
+            &      nitrfac  (A2D(0),jpk), nitrfac2(A2D(0),jpk),  &
+            &      prodcal  (A2D(0),jpk), xdiss   (A2D(0),jpk),  &
+            &      prodpoc  (A2D(0),jpk), conspoc (A2D(0),jpk),  &
+            &      prodgoc  (A2D(0),jpk), consgoc (A2D(0),jpk),  &
+            &      blim     (A2D(0),jpk), consfe3 (A2D(0),jpk),  &
+            &      xfecolagg(A2D(0),jpk), xcoagfe (A2D(0),jpk), STAT=ierr(4) )
 
          !* Carbonate chemistry
-         ALLOCATE( ak13  (jpi,jpj,jpk)  ,                          &
-            &      ak23(jpi,jpj,jpk)    , aksp  (jpi,jpj,jpk) ,    &
-            &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,    &
-            &      aphscale(jpi,jpj,jpk),                         STAT=ierr(5) )
+         ALLOCATE( ak13  (A2D(0),jpk)  ,                          &
+            &      ak23(A2D(0),jpk)    , aksp  (A2D(0),jpk) ,    &
+            &      hi  (A2D(0),jpk)    , excess(A2D(0),jpk) ,    &
+            &      aphscale(A2D(0),jpk),                         STAT=ierr(5) )
          !
          !* Temperature dependency of SMS terms
-         ALLOCATE( tgfunc (jpi,jpj,jpk) , tgfunc2(jpi,jpj,jpk),   STAT=ierr(6) )
+         ALLOCATE( tgfunc (A2D(0),jpk) , tgfunc2(A2D(0),jpk),   STAT=ierr(6) )
          !
          !* Sinking speed
-         ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4 (jpi,jpj,jpk),   STAT=ierr(7) )
+         ALLOCATE( wsbio3 (A2D(0),jpk) , wsbio4 (A2D(0),jpk),   STAT=ierr(7) )
 
          !*  Size of phytoplankton cells
-         ALLOCATE( sizen (jpi,jpj,jpk), sized (jpi,jpj,jpk),        &
-           &       sizena(jpi,jpj,jpk), sizeda(jpi,jpj,jpk),      STAT=ierr(8) )
+         ALLOCATE( sizen (A2D(0),jpk), sized (A2D(0),jpk),        &
+           &       sizena(A2D(0),jpk), sizeda(A2D(0),jpk),      STAT=ierr(8) )
          ! 
-         ALLOCATE( plig(jpi,jpj,jpk)  ,                           STAT=ierr(9) )
+         ALLOCATE( plig(A2D(0),jpk)  ,                           STAT=ierr(9) )
       ENDIF
       !
       IF( ln_p5z ) THEN
          ! PISCES-QUOTA specific part      
-         ALLOCATE( epico(jpi,jpj,jpk)   , epicom(jpi,jpj,jpk) ,   STAT=ierr(10) ) 
+         ALLOCATE( epico(A2D(0),jpk)   , epicom(A2D(0),jpk) ,   STAT=ierr(10) ) 
 
          !*  Size of phytoplankton cells
-         ALLOCATE( sizep(jpi,jpj,jpk), sizepa(jpi,jpj,jpk),       STAT=ierr(11) )
+         ALLOCATE( sizep(A2D(0),jpk), sizepa(A2D(0),jpk),       STAT=ierr(11) )
       ENDIF
       !
       sms_pisces_alloc = MAXVAL( ierr )
diff --git a/src/TOP/PISCES/trcice_pisces.F90 b/src/TOP/PISCES/trcice_pisces.F90
index 8ce76bc76b9740e84bcd3bac29fea74ad22daae7..06e458a54852ae81a83dc6da4f26c0d7d9c97b60 100644
--- a/src/TOP/PISCES/trcice_pisces.F90
+++ b/src/TOP/PISCES/trcice_pisces.F90
@@ -19,6 +19,8 @@ MODULE trcice_pisces
 
    PUBLIC   trc_ice_ini_pisces ! called by trcini.F90 module
 
+      !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
    !! $Id: trcice_pisces.F90 10794 2019-03-22 09:25:28Z cetlod $ 
@@ -283,15 +285,15 @@ CONTAINS
       ENDIF
 ! 
       DO jn = jp_pcs0, jp_pcs1
-         IF( cn_trc_o(jn) == 'GL ' ) trc_o(:,:,jn) = zpisc(jn,1)  ! Global case
+         IF( cn_trc_o(jn) == 'GL ' ) trc_o(A2D(0),jn) = zpisc(jn,1)  ! Global case
          IF( cn_trc_o(jn) == 'AA ' ) THEN 
-            WHERE( gphit(:,:) >= 0._wp ) ; trc_o(:,:,jn) = zpisc(jn,2) ; END WHERE ! Arctic 
-            WHERE( gphit(:,:) <  0._wp ) ; trc_o(:,:,jn) = zpisc(jn,3) ; END WHERE ! Antarctic 
+            WHERE( gphit(A2D(0)) >= 0._wp ) ; trc_o(A2D(0),jn) = zpisc(jn,2) ; END WHERE ! Arctic 
+            WHERE( gphit(A2D(0)) <  0._wp ) ; trc_o(A2D(0),jn) = zpisc(jn,3) ; END WHERE ! Antarctic 
          ENDIF
          IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN     !  Baltic Sea particular case for ORCA configurations
-             WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.    &
-                    54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp )
-                    trc_o(:,:,jn) = zpisc(jn,4)
+             WHERE( 14._wp <= glamt(A2D(0)) .AND. glamt(A2D(0)) <= 32._wp .AND.    &
+                    54._wp <= gphit(A2D(0)) .AND. gphit(A2D(0)) <= 66._wp )
+                    trc_o(A2D(0),jn) = zpisc(jn,4)
             END WHERE
          ENDIF 
       ENDDO
@@ -321,16 +323,16 @@ CONTAINS
       DO jn = jp_pcs0, jp_pcs1
          !-- Everywhere but in the Baltic
          IF ( trc_ice_ratio(jn) >= -1._wp ) THEN ! no prescribed conc. ; typically everything but iron) 
-            trc_i(:,:,jn) = zratio(jn,1) * trc_o(:,:,jn) 
+            trc_i(A2D(0),jn) = zratio(jn,1) * trc_o(A2D(0),jn) 
          ELSE                                    ! prescribed concentration
-            trc_i(:,:,jn) = trc_ice_prescr(jn)
+            trc_i(A2D(0),jn) = trc_ice_prescr(jn)
          ENDIF
          !-- Baltic
          IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN     
             IF ( trc_ice_ratio(jn) >= - 1._wp ) THEN ! no prescribed conc. ; typically everything but iron) 
-               WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.    &
-                      54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp )
-                     trc_i(:,:,jn) = zratio(jn,2) * trc_o(:,:,jn) 
+               WHERE( 14._wp <= glamt(A2D(0)) .AND. glamt(A2D(0)) <= 32._wp .AND.    &
+                      54._wp <= gphit(A2D(0)) .AND. gphit(A2D(0)) <= 66._wp )
+                     trc_i(A2D(0),jn) = zratio(jn,2) * trc_o(A2D(0),jn) 
                END WHERE
             ENDIF
          ENDIF
diff --git a/src/TOP/PISCES/trcwri_pisces.F90 b/src/TOP/PISCES/trcwri_pisces.F90
index e5cce04d5a93d31375531ac4d92fb9e5e00a2520..d44b7b11eaa800b2cb374084b3cd80f5dae4fb93 100644
--- a/src/TOP/PISCES/trcwri_pisces.F90
+++ b/src/TOP/PISCES/trcwri_pisces.F90
@@ -38,7 +38,7 @@ CONTAINS
       CHARACTER (len=20)           :: cltra
       REAL(wp)                     :: zfact
       INTEGER                      :: ji, jj, jk, jn
-      REAL(wp), DIMENSION(jpi,jpj) :: zdic, zo2min, zdepo2min
+      REAL(wp), DIMENSION(A2D(0)) :: zdic, zo2min, zdepo2min
       !!---------------------------------------------------------------------
  
       ! write the tracer concentrations in the file
@@ -60,15 +60,19 @@ CONTAINS
          IF( iom_use( "INTDIC" ) ) THEN                     !   DIC content in kg/m2
             zdic(:,:) = 0.
             DO jk = 1, jpkm1
-               zdic(:,:) = zdic(:,:) + tr(:,:,jk,jpdic,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) * 12.
+               DO_2D( 0, 0, 0, 0 )
+                 zdic(ji,jj) = zdic(ji,jj) + tr(ji,jj,jk,jpdic,Kmm) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * 12.
+               END_2D
             ENDDO
-            CALL iom_put( 'INTDIC', zdic )     
+            CALL iom_put( 'INTDIC', zdic )
          ENDIF
          !
-         IF( iom_use( "O2MIN" ) .OR. iom_use ( "ZO2MIN" ) ) THEN  ! Oxygen minimum concentration and depth 
-            zo2min   (:,:) = tr(:,:,1,jpoxy,Kmm) * tmask(:,:,1)
-            zdepo2min(:,:) = gdepw(:,:,1,Kmm)   * tmask(:,:,1)
-            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 ) 
+         IF( iom_use( "O2MIN" ) .OR. iom_use ( "ZO2MIN" ) ) THEN  ! Oxygen minimum concentration and depth
+            DO_2D( 0, 0, 0, 0 )
+               zo2min   (ji,jj) = tr(ji,jj,1,jpoxy,Kmm) * tmask(ji,jj,1)
+               zdepo2min(ji,jj) = gdepw(ji,jj,1,Kmm)    * tmask(ji,jj,1)
+            END_2D
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                IF( tmask(ji,jj,jk) == 1 ) then
                   IF( tr(ji,jj,jk,jpoxy,Kmm) < zo2min(ji,jj) ) then
                      zo2min   (ji,jj) = tr(ji,jj,jk,jpoxy,Kmm)