diff --git a/cfgs/SHARED/field_def_nemo-pisces.xml b/cfgs/SHARED/field_def_nemo-pisces.xml
index 41e65d48fce0445af8ba6fbe5d1a0429c49dd9e0..23ec65218916b654e53034d2d242bafc562ae0a2 100644
--- a/cfgs/SHARED/field_def_nemo-pisces.xml
+++ b/cfgs/SHARED/field_def_nemo-pisces.xml
@@ -189,10 +189,10 @@
     <field id="PFeP"        long_name="Primary production of pico iron"         unit="molC/m3/s"  grid_ref="grid_T_3D" />
     <field id="PFeD"        long_name="Primary production of diatoms iron"      unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
     <field id="xfracal"     long_name="Calcifying fraction"                     unit="1"          grid_ref="grid_T_3D_inner" />
-    <field id="PCAL"        long_name="Calcite production"                      unit="mol/m3/s"   grid_ref="grid_T_3D" />
-    <field id="GRAZ1"       long_name="Grazing by microzooplankton"             unit="mol/m3/s"   grid_ref="grid_T_3D" />
-    <field id="GRAZ2"       long_name="Grazing by mesozooplankton"              unit="mol/m3/s"   grid_ref="grid_T_3D" />
-    <field id="REMIN"       long_name="Oxic remineralization of OM"             unit="mol/m3/s"   grid_ref="grid_T_3D" />
+    <field id="PCAL"        long_name="Calcite production"                      unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
+    <field id="GRAZ1"       long_name="Grazing by microzooplankton"             unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
+    <field id="GRAZ2"       long_name="Grazing by mesozooplankton"              unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
+    <field id="REMIN"       long_name="Oxic remineralization of OM"             unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
     <field id="DENIT"       long_name="Anoxic remineralization of OM"           unit="mol/m3/s"   grid_ref="grid_T_3D" />
     <field id="REMINP"       long_name="Oxic remineralization rate of POC"      unit="d-1"        grid_ref="grid_T_3D" />
     <field id="REMING"       long_name="Oxic remineralization rate of GOC"      unit="d-1"        grid_ref="grid_T_3D" />
@@ -253,18 +253,18 @@
     <field id="FECOLL"      long_name="Colloidal Pumping of FeL"                unit="mmol-FeL/m3/s"  grid_ref="grid_T_3D" />
     <field id="LGWCOLL"     long_name="Coagulation loss of ligands"             unit="mmol-L/m3/s"  grid_ref="grid_T_3D" />
     <field id="REMINF"      long_name="Oxic remineralization suppy of Fe"       unit="mmol-Fe/m3/s"  grid_ref="grid_T_3D" />
-    <field id="BACT"        long_name="Bacterial Biomass"                       unit="mmol/m3"  grid_ref="grid_T_3D" />
-    <field id="FEBACT"      long_name="Bacterial uptake of Fe"                  unit="molFe/m3/s"  grid_ref="grid_T_3D" />
+    <field id="BACT"        long_name="Bacterial Biomass"                       unit="mmol/m3"  grid_ref="grid_T_3D_inner" />
+    <field id="FEBACT"      long_name="Bacterial uptake of Fe"                  unit="molFe/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="FEPREC"      long_name="Precipitation of Fe"                     unit="molFe/m3/s"  grid_ref="grid_T_3D" />
     <field id="LPRODR"      long_name="OM remineralisation ligand production rate" unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
     <field id="LPRODP"      long_name="phytoplankton ligand production rate"    unit="nmol-L/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="LIGREM"      long_name="Remineralisation loss of ligands"        unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
     <field id="LIGPR"       long_name="Photochemical loss of ligands"           unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
     <field id="LDETP"       long_name="Ligand destruction during phytoplankton uptake" unit="nmol-L/m3/s"  grid_ref="grid_T_3D_inner" />
-    <field id="LPRODZ2"     long_name="mesozooplankton ligand production rate"  unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
+    <field id="LPRODZ2"     long_name="mesozooplankton ligand production rate"  unit="nmol-L/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="LPRODZ"      long_name="microzooplankton ligand production rate" unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
-    <field id="FEZOO"       long_name="microzooplankton iron recycling rate"    unit="nmol-FeL/m3/s"  grid_ref="grid_T_3D" />
-    <field id="FEZOO2"      long_name="mesozooplankton iron recycling rate"     unit="nmol-FeL/m3/s"  grid_ref="grid_T_3D" />
+    <field id="FEZOO"       long_name="microzooplankton iron recycling rate"    unit="nmol-FeL/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="FEZOO2"      long_name="mesozooplankton iron recycling rate"     unit="nmol-FeL/m3/s"  grid_ref="grid_T_3D_inner" />
 
     <!-- PISCES tracers trends -->
     <field id="INTdtAlk"    long_name="Vertically int. of change of alkalinity"             unit="mol/m2/s"                       />
diff --git a/src/TOP/PISCES/P4Z/p4zbio.F90 b/src/TOP/PISCES/P4Z/p4zbio.F90
index d2a21ea1cd049af3a7c8ae96ce5feadaf2e14ff0..65f61f4f90552a637b0f6c8f53b2058bc9b708fa 100644
--- a/src/TOP/PISCES/P4Z/p4zbio.F90
+++ b/src/TOP/PISCES/P4Z/p4zbio.F90
@@ -72,7 +72,7 @@ CONTAINS
       ! OF PHYTOPLANKTON AND DETRITUS. Shear rate is supposed to equal 1
       ! in the mixed layer and 0.1 below the mixed layer.
       xdiss(:,:,:) = 1.
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          IF( gdepw(ji,jj,jk+1,Kmm) > hmld(ji,jj) )   xdiss(ji,jj,jk) = 0.01
       END_3D
 
diff --git a/src/TOP/PISCES/P4Z/p4zflx.F90 b/src/TOP/PISCES/P4Z/p4zflx.F90
index 051985ba2037053eb13b4094b8a8aff532c7f2de..1d67f902a0a8868b43261222cb1f3cb1668d6f16 100644
--- a/src/TOP/PISCES/P4Z/p4zflx.F90
+++ b/src/TOP/PISCES/P4Z/p4zflx.F90
@@ -83,7 +83,7 @@ CONTAINS
       REAL(wp) ::   zph, zdic, zsch_o2, zsch_co2
       REAL(wp) ::   zyr_dec, zdco2dt
       CHARACTER (len=25) ::   charout
-      REAL(wp), DIMENSION(jpi,jpj) ::   zkgco2, zkgo2, zh2co3, zoflx,  zpco2atm, zpco2oce  
+      REAL(wp), DIMENSION(A2D(0)) ::   zkgco2, zkgo2, zh2co3, zoflx,  zpco2atm, zpco2oce  
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_flx')
@@ -110,7 +110,7 @@ CONTAINS
 
       IF( l_co2cpl )   satmco2(:,:) = atm_co2(:,:)
 
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
          zrhd = rhd(ji,jj,1) + 1._wp
          zdic  = tr(ji,jj,1,jpdic,Kbb)
@@ -126,7 +126,7 @@ CONTAINS
       ! FIRST COMPUTE GAS EXCHANGE COEFFICIENTS
       ! -------------------------------------------
 
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          ztc  = MIN( 35., ts(ji,jj,1,jp_tem,Kmm) )
          ztc2 = ztc * ztc
          ztc3 = ztc * ztc2 
@@ -145,7 +145,7 @@ CONTAINS
       END_2D
 
 
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          ztkel = tempis(ji,jj,1) + 273.15
          zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
          zvapsw    = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal)
diff --git a/src/TOP/PISCES/P4Z/p4zmeso.F90 b/src/TOP/PISCES/P4Z/p4zmeso.F90
index 0033e0b8939a11324152aeec6fe0a82eeb086f71..f9f1f149d835f927c045e76ca8fbf3f7158d3d1e 100644
--- a/src/TOP/PISCES/P4Z/p4zmeso.F90
+++ b/src/TOP/PISCES/P4Z/p4zmeso.F90
@@ -89,10 +89,11 @@ CONTAINS
       REAL(wp) :: zgrazfffp, zgrazfffg, zgrazffep, zgrazffeg, zrum, zcodel, zargu, zval, zdep
       REAL(wp) :: zsigma, zdiffdn, ztmp1, ztmp2, ztmp3, ztmp4, ztmptot, zmigthick 
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing, zfezoo2
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrarem, zgraref, zgrapoc, zgrapof, zgrabsi
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zgrazing, zfezoo2
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zgrarem, zgraref, zgrapoc, zgrapof, zgrabsi
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigrem, zgramigref, zgramigpoc, zgramigpof
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigbsi
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_meso')
@@ -108,7 +109,7 @@ CONTAINS
       ! ---------------------------------------------
       IF (ln_dvm_meso) CALL p4z_meso_depmig( Kbb, Kmm )
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          zcompam   = MAX( ( tr(ji,jj,jk,jpmes,Kbb) - 1.e-9 ), 0.e0 )
          zfact     = xstep * tgfunc2(ji,jj,jk) * zcompam
 
@@ -345,7 +346,7 @@ CONTAINS
         ! This fraction is sumed over the euphotic zone and is removed from 
         ! the fluxes driven by mesozooplankton in the euphotic zone.
         ! --------------------------------------------------------------------
-        DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk)
+        DO_3D( 0, 0, 0, 0, 1, jpk)
             zmigreltime = (1. - strn(ji,jj))
             zmigthick   = (1. - zmigreltime ) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
             IF ( gdept(ji,jj,jk,Kmm) <= heup(ji,jj) ) THEN
@@ -366,7 +367,7 @@ CONTAINS
          ! The inorganic and organic fluxes induced by migrating organisms are added at the 
          ! the migration depth (corresponding indice is set by kmig)
          ! --------------------------------------------------------------------------------
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
             IF( tmask(ji,jj,1) == 1.) THEN
                jkt = kmig(ji,jj)
                zdep = 1. / e3t(ji,jj,jkt,Kmm)
@@ -387,7 +388,7 @@ CONTAINS
       !   Update the arrays TRA which contain the biological sources and sinks
       !   This only concerns the variables which are affected by DVM (inorganic 
       !   nutrients, DOC agands, and particulate organic carbon). 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk)
+      DO_3D( 0, 0, 0, 0, 1, jpk)
          tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zgrarem(ji,jj,jk) * sigma2
          tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zgrarem(ji,jj,jk) * sigma2
          tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zgrarem(ji,jj,jk) * ( 1. - sigma2 )
@@ -408,11 +409,30 @@ CONTAINS
       !
       ! Write the output
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-        CALL iom_put( "PCAL"  , prodcal(:,:,:) * 1.e+3  * rfact2r * tmask(:,:,:) )  !  Calcite production 
-        CALL iom_put( "GRAZ2" , zgrazing(:,:,:) * 1.e+3  * rfact2r * tmask(:,:,:) ) ! Total grazing of phyto by zoo
-        CALL iom_put( "FEZOO2", zfezoo2(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
-        IF( ln_ligand ) &
-         & CALL iom_put( "LPRODZ2", zgrarem(ji,jj,jk) * ( 1. - sigma2 ) * ldocz * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:)  )
+        !
+        ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+        !
+        IF( iom_use ( "PCAL" ) ) THEN   ! Calcite production
+          zw3d(A2D(0),1:jpkm1) = prodcal(A2D(0),1:jpkm1) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "PCAL", zw3d )
+        ENDIF
+        !
+        IF( iom_use ( "GRAZ2" ) ) THEN   ! Total grazing of phyto by zoo
+          zw3d(A2D(0),1:jpkm1) = zgrazing(A2D(0),1:jpkm1) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "GRAZ2", zw3d )
+        ENDIF
+        !
+        IF( iom_use ( "FEZOO2" ) ) THEN   
+          zw3d(A2D(0),1:jpkm1) = zfezoo2(A2D(0),1:jpkm1) * 1e9 * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "FEZOO2", zw3d )
+        ENDIF
+        !
+        IF( ln_ligand .AND. iom_use( "LPRODZ2" )  ) THEN
+            zw3d(A2D(0),1:jpkm1) = zgrarem(A2D(0),1:jpkm1) * ( 1. - sigma2 ) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
+            CALL iom_put( "LPRODZ2"  , zw3d * ldocz * 1e9 )
+        ENDIF
+            !
+        DEALLOCATE( zw3d )
       ENDIF
       !
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
@@ -517,7 +537,7 @@ CONTAINS
       ! Compute the averaged values of oxygen, temperature over the domain 
       ! 150m to 500 m depth.
       ! ------------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk)
+      DO_3D( 0, 0, 0, 0, 1, jpk)
          IF( tmask(ji,jj,jk) == 1.) THEN
             IF( gdept(ji,jj,jk,Kmm) >= 150. .AND. gdept(ji,jj,jk,kmm) <= 500.) THEN
                oxymoy(ji,jj)  = oxymoy(ji,jj)  + tr(ji,jj,jk,jpoxy,Kbb) * 1E6 * e3t(ji,jj,jk,Kmm)
@@ -530,7 +550,7 @@ CONTAINS
       ! Compute the difference between surface values and the mean values in the mesopelagic
       ! domain
       ! ------------------------------------------------------------------------------------
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          z1dep = 1. / ( zdepmoy(ji,jj) + rtrn )
          oxymoy(ji,jj)  = tr(ji,jj,1,jpoxy,Kbb) * 1E6 - oxymoy(ji,jj)  * z1dep
          tempmoy(ji,jj) = ts(ji,jj,1,jp_tem,Kmm)      - tempmoy(ji,jj) * z1dep
@@ -539,7 +559,7 @@ CONTAINS
       ! Computation of the migration depth based on the parameterization of 
       ! Bianchi et al. (2013)
       ! -------------------------------------------------------------------
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          IF( tmask(ji,jj,1) == 1. ) THEN
             ztotchl = ( tr(ji,jj,1,jpnch,Kbb) + tr(ji,jj,1,jpdch,Kbb) ) * 1E6
             depmig(ji,jj) = 398. - 0.56 * oxymoy(ji,jj) -115. * log10(ztotchl) + 0.36 * hmld(ji,jj) -2.4 * tempmoy(ji,jj)
@@ -548,7 +568,7 @@ CONTAINS
       ! 
       ! Computation of the corresponding jk indice 
       ! ------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          IF( depmig(ji,jj) >= gdepw(ji,jj,jk,Kmm) .AND. depmig(ji,jj) < gdepw(ji,jj,jk+1,Kmm) ) THEN
              kmig(ji,jj) = jk
           ENDIF
@@ -560,7 +580,7 @@ CONTAINS
       ! to 0. Thus, to avoid that problem, the migration depth is adjusted so
       ! that it falls above the OMZ
       ! -----------------------------------------------------------------------
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          IF( tr(ji,jj,kmig(ji,jj),jpoxy,Kbb) < 5E-6 ) THEN
             DO jk = kmig(ji,jj),1,-1
                IF( tr(ji,jj,jk,jpoxy,Kbb) >= 5E-6 .AND. tr(ji,jj,jk+1,jpoxy,Kbb)  < 5E-6) THEN
diff --git a/src/TOP/PISCES/P4Z/p4zmicro.F90 b/src/TOP/PISCES/P4Z/p4zmicro.F90
index f91b04d7f185c17ab3e6145823c730a9b39a334e..71255455359da6daf82531df804363767aca82c3 100644
--- a/src/TOP/PISCES/P4Z/p4zmicro.F90
+++ b/src/TOP/PISCES/P4Z/p4zmicro.F90
@@ -78,8 +78,9 @@ CONTAINS
       REAL(wp) :: zrespz, ztortz, zgrasratf, zgrasratn
       REAL(wp) :: zgraznc, zgrazpoc, zgrazdc, zgrazpof, zgrazdf, zgraznf
       REAL(wp) :: zsigma, zdiffdn, ztmp1, ztmp2, ztmp3, ztmptot, zproport
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing, zfezoo
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zgrazing, zfezoo
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zzligprod
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zw3d
       CHARACTER (len=25) :: charout
 
       !!---------------------------------------------------------------------
@@ -91,7 +92,7 @@ CONTAINS
          zzligprod(:,:,:) = 0._wp
       ENDIF
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          zcompaz = MAX( ( tr(ji,jj,jk,jpzoo,Kbb) - 1.e-9 ), 0.e0 )
          zfact   = xstep * tgfunc2(ji,jj,jk) * zcompaz
 
@@ -257,15 +258,24 @@ CONTAINS
       END_3D
       !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-        IF( iom_use("GRAZ1") ) THEN  !   Total grazing of phyto by zooplankton
-           zgrazing(:,:,jpk) = 0._wp   ; CALL iom_put( "GRAZ1" , zgrazing(:,:,:) * 1.e+3  * rfact2r * tmask(:,:,:) ) 
-         ENDIF
-         IF( iom_use("FEZOO") ) THEN  
-           zfezoo (:,:,jpk) = 0._wp    ; CALL iom_put( "FEZOO", zfezoo(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
-         ENDIF
-         IF( ln_ligand ) THEN
-            zzligprod(:,:,jpk) = 0._wp ; CALL iom_put( "LPRODZ", zzligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:))
-         ENDIF
+        !
+        ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+        !
+        IF( iom_use ( "GRAZ1" ) ) THEN   ! Total grazing of phyto by zooplankton
+          zw3d(A2D(0),1:jpkm1) = zgrazing(A2D(0),1:jpkm1) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "GRAZ1", zw3d )
+        ENDIF
+        !
+        IF( iom_use ( "FEZOO" ) ) THEN 
+          zw3d(A2D(0),1:jpkm1) = zfezoo(A2D(0),1:jpkm1)  * 1e9 * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "FEZOO", zw3d )
+        ENDIF
+        !
+        IF( ln_ligand .AND. iom_use( "LPRODZ2" ) ) THEN
+           zzligprod(:,:,jpk) = 0._wp ; CALL iom_put( "LPRODZ", zzligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:))
+        ENDIF
+        !
+        DEALLOCATE( zw3d )
       ENDIF
       !
       IF(sn_cfctl%l_prttrc) THEN      ! print mean trends (used for debugging)
diff --git a/src/TOP/PISCES/P4Z/p4zrem.F90 b/src/TOP/PISCES/P4Z/p4zrem.F90
index c63a6e63f59e3fd8d803181902ac7a046a6950e1..13d813117057f67577890fad94f12a49a0bcb282 100644
--- a/src/TOP/PISCES/P4Z/p4zrem.F90
+++ b/src/TOP/PISCES/P4Z/p4zrem.F90
@@ -76,8 +76,9 @@ CONTAINS
       REAL(wp) ::   zbactfer, zonitr, zrfact2
       REAL(wp) ::   zammonic, zoxyremc, zosil, ztem, zdenitnh4, zolimic
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepbac, zolimi, zfacsi, zfacsib, zdepeff, zfebact
-      REAL(wp), DIMENSION(jpi,jpj    ) :: ztempbac
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zdepbac, zolimi, zfacsi, zfacsib, zdepeff, zfebact
+      REAL(wp), DIMENSION(A2D(0)    ) :: ztempbac
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_rem')
@@ -90,11 +91,11 @@ CONTAINS
 
       ! Computation of the mean bacterial concentration
       ! this parameterization has been deduced from a model version
-      ! that was modeling explicitely bacteria. This is a very old param 
+      ! that was modeling explicitely bacteria. This is a very old parame
       ! that will be very soon updated based on results from a much more
       ! recent version of PISCES with bacteria.
       ! ----------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          zdep = MAX( hmld(ji,jj), heup_01(ji,jj) )
          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
@@ -108,7 +109,7 @@ CONTAINS
          ENDIF
       END_3D
 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          ! DOC ammonification. Depends on depth, phytoplankton biomass
          ! and a limitation term which is supposed to be a parameterization of the bacterial activity. 
          ! --------------------------------------------------------------------------
@@ -152,7 +153,7 @@ CONTAINS
          ENDIF
       END_3D
 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          ! NH4 nitrification to NO3. Ceased for oxygen concentrations
          ! below 2 umol/L. Inhibited at strong light 
          ! ----------------------------------------------------------
@@ -174,7 +175,7 @@ CONTAINS
          CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
 
          ! Bacterial uptake of iron. No iron is available in DOC. So
          ! Bacteries are obliged to take up iron from the water. Some
@@ -202,7 +203,7 @@ CONTAINS
       ! Initialization of the array which contains the labile fraction
       ! of bSi. Set to a constant in the upper ocean
       ! ---------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          ! Remineralization rate of BSi dependent on T and saturation
          ! The parameterization is taken from Ridgwell et al. (2002) 
          ! ---------------------------------------------------------
@@ -236,20 +237,35 @@ CONTAINS
          WRITE(charout, FMT="('rem3')")
          CALL prt_ctl_info( charout, cdcomp = 'top' )
          CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
-       ENDIF
+      ENDIF
 
       IF( knt == nrdttrc ) THEN
           zrfact2 = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
           !
-          IF( iom_use( "REMIN" ) )  THEN !  Remineralisation rate
-             zolimi(:,:,jpk) = 0. ; CALL iom_put( "REMIN"  , zolimi(:,:,:) * tmask(:,:,:) * zrfact2  )
+          ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+          !
+          IF( iom_use ( "REMIN" ) ) THEN   ! Remineralisation rate
+            zw3d(A2D(0),1:jpkm1) = zolimi(A2D(0),1:jpkm1) * rfact2r * tmask(A2D(0),1:jpkm1)
+            CALL iom_put( "REMIN", zw3d )
           ENDIF
-          CALL iom_put( "DENIT"  , denitr(:,:,:) * rdenit * rno3 * tmask(:,:,:) * zrfact2 ) ! Denitrification 
-          IF( iom_use( "BACT" ) )  THEN ! Bacterial biomass
-             zdepbac(:,:,jpk) = 0.  ;   CALL iom_put( "BACT", zdepbac(:,:,:) * 1.E6 * tmask(:,:,:) )
+          !
+          IF( iom_use ( "BACT" ) ) THEN   ! Bacterial biomass
+            zw3d(A2D(0),1:jpkm1) = zdepbac(A2D(0),1:jpkm1) * 1.E6 * tmask(A2D(0),1:jpkm1)
+            CALL iom_put( "BACT", zw3d )
           ENDIF
-          CALL iom_put( "FEBACT" , zfebact(:,:,:) * 1E9 * tmask(:,:,:) * zrfact2  )
-       ENDIF
+          !
+          IF( iom_use ( "FEBACT" ) ) THEN 
+            zw3d(A2D(0),1:jpkm1) = zfebact(A2D(0),1:jpkm1) * 1E9 * rfact2r * tmask(A2D(0),1:jpkm1)
+            CALL iom_put( "FEBACT", zw3d )
+          ENDIF
+          !
+          IF( iom_use( "DENIT" ) )  THEN ! Denitrification
+             denitr(:,:,jpk) = 0._wp  
+             CALL iom_put( "DENIT", denitr(:,:,:) * rdenit * rno3 * tmask(:,:,:) * zrfact2 )
+          ENDIF
+          !
+          DEALLOCATE( zw3d )
+      ENDIF
       !
       IF( ln_timing )   CALL timing_stop('p4z_rem')
       !
diff --git a/src/TOP/PISCES/P4Z/p4zsink.F90 b/src/TOP/PISCES/P4Z/p4zsink.F90
index 53db2671e02ac76bedea89eb4c7cfd259172fad5..e900c82cafe7c5fcd0f3e59122a6144fd52557e4 100644
--- a/src/TOP/PISCES/P4Z/p4zsink.F90
+++ b/src/TOP/PISCES/P4Z/p4zsink.F90
@@ -86,7 +86,7 @@ CONTAINS
       ! CaCO3 and bSi are supposed to sink at the big particles speed 
       ! due to their high density
       ! ---------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          zmax  = MAX( heup_01(ji,jj), hmld(ji,jj) )
          zfact = MAX( 0., gdepw(ji,jj,jk+1,Kmm) - zmax ) / wsbio2scale
          wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact