From f24a87f04eeb8513f57737a382e40060b67cdfbf Mon Sep 17 00:00:00 2001
From: Renaud Person <renaud.person@locean.ipsl.fr>
Date: Mon, 25 Jul 2022 17:35:34 +0200
Subject: [PATCH] Continuation of PISCES-p4z/p5z routines halo cleanup

---
 cfgs/SHARED/field_def_nemo-pisces.xml |  40 ++---
 src/TOP/PISCES/P4Z/p4zfechem.F90      |  44 ++++--
 src/TOP/PISCES/P4Z/p4zflx.F90         |  54 +++++--
 src/TOP/PISCES/P4Z/p4zligand.F90      |  12 +-
 src/TOP/PISCES/P4Z/p4zmeso.F90        |   4 +-
 src/TOP/PISCES/P4Z/p4zmicro.F90       |  12 +-
 src/TOP/PISCES/P4Z/p4zopt.F90         |  20 +--
 src/TOP/PISCES/P4Z/p4zpoc.F90         |  23 ++-
 src/TOP/PISCES/P4Z/p4zprod.F90        |   4 +-
 src/TOP/PISCES/P4Z/p4zsed.F90         |  57 ++++---
 src/TOP/PISCES/P4Z/p4zsink.F90        |  44 ++++--
 src/TOP/PISCES/P4Z/p4zsms.F90         |   2 +-
 src/TOP/PISCES/P4Z/p5zlim.F90         |  93 +++++++++---
 src/TOP/PISCES/P4Z/p5zmeso.F90        |  57 ++++---
 src/TOP/PISCES/P4Z/p5zmicro.F90       |  34 +++--
 src/TOP/PISCES/P4Z/p5zprod.F90        | 207 +++++++++++++++++++-------
 src/TOP/PISCES/SED/sedchem.F90        |   2 +-
 src/TOP/PISCES/SED/sedsfc.F90         |   2 +-
 src/TOP/PISCES/SED/trcdmp_sed.F90     |   2 +-
 19 files changed, 507 insertions(+), 206 deletions(-)

diff --git a/cfgs/SHARED/field_def_nemo-pisces.xml b/cfgs/SHARED/field_def_nemo-pisces.xml
index 23ec65218..a94dec2ab 100644
--- a/cfgs/SHARED/field_def_nemo-pisces.xml
+++ b/cfgs/SHARED/field_def_nemo-pisces.xml
@@ -179,14 +179,14 @@
     <field id="DCAL"        long_name="Calcite dissolution"                     unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
     <field id="PAR"         long_name="Photosynthetically Available Radiation"  unit="W/m2"       grid_ref="grid_T_3D" />
     <field id="PPPHYN"      long_name="Primary production of nanophyto"         unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
-    <field id="PPPHYP"      long_name="Primary production of picophyto"         unit="molC/m3/s"  grid_ref="grid_T_3D" />
+    <field id="PPPHYP"      long_name="Primary production of picophyto"         unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="PPPHYD"      long_name="Primary production of diatoms"           unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="PPNEWN"      long_name="New Primary production of nanophyto"     unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
-    <field id="PPNEWP"      long_name="New Primary production of picophyto"     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_inner" />
     <field id="PPNEWD"      long_name="New Primary production of diatoms"       unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="PBSi"        long_name="Primary production of Si diatoms"        unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="PFeN"        long_name="Primary production of nano iron"         unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
-    <field id="PFeP"        long_name="Primary production of pico 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_inner" />
     <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_inner" />
@@ -194,12 +194,12 @@
     <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" />
+    <field id="REMINP"       long_name="Oxic remineralization rate of POC"      unit="d-1"        grid_ref="grid_T_3D_inner" />
+    <field id="REMING"       long_name="Oxic remineralization rate of GOC"      unit="d-1"        grid_ref="grid_T_3D_inner" />
     <field id="Nfix"        long_name="Nitrogen fixation"                       unit="mol/m3/s"   grid_ref="grid_T_3D" />
     <field id="Mumax"       long_name="Maximum growth rate"                     unit="s-1"        grid_ref="grid_T_3D_inner" />
     <field id="MuN"         long_name="Realized growth rate for nanophyto"      unit="s-1"        grid_ref="grid_T_3D_inner" />
-    <field id="MuP"         long_name="Realized growth rate for picophyto"      unit="s-1"        grid_ref="grid_T_3D" />
+    <field id="MuP"         long_name="Realized growth rate for picophyto"      unit="s-1"        grid_ref="grid_T_3D_inner" />
     <field id="MuD"         long_name="Realized growth rate for diatomes"       unit="s-1"        grid_ref="grid_T_3D_inner" />
     <field id="MunetN"      long_name="Net growth rate for nanophyto"           unit="s-1"        grid_ref="grid_T_3D" />
     <field id="MunetP"      long_name="Net growth rate for picophyto"           unit="s-1"        grid_ref="grid_T_3D" />
@@ -211,7 +211,7 @@
     <field id="LPFe"        long_name="Iron limitation term in Picophyto"       unit="-"          grid_ref="grid_T_3D_inner" />
     <field id="LDFe"        long_name="Iron limitation term in Diatoms"         unit=""           grid_ref="grid_T_3D_inner" />
     <field id="LNlight"     long_name="Light limitation term in Nanophyto"      unit=""           grid_ref="grid_T_3D_inner" />
-    <field id="LPlight"     long_name="Light limitation term in Picophyto"      unit="-"          grid_ref="grid_T_3D" />
+    <field id="LPlight"     long_name="Light limitation term in Picophyto"      unit="-"          grid_ref="grid_T_3D_inner" />
     <field id="LDlight"     long_name="Light limitation term in Diatoms"        unit=""           grid_ref="grid_T_3D_inner" />
     <field id="SIZEN"       long_name="Mean relative size of nanophyto."        unit="-"          grid_ref="grid_T_3D_inner" />
     <field id="SIZEP"       long_name="Mean relative size of picophyto."        unit="-"          grid_ref="grid_T_3D_inner" />
@@ -227,9 +227,9 @@
     <field id="Biron"       long_name="Bioavailable iron"                       unit="nmol/m3"    grid_ref="grid_T_3D" />
     <field id="Sdenit"      long_name="Nitrate reduction in the sediments"      unit="mol/m2/s"                        />
     <field id="Ironice"     long_name="Iron input/uptake due to sea ice"        unit="mol/m2/s"                        />
-    <field id="SedCal"      long_name="Calcite burial in the sediments"         unit="molC/m2/s"                       />
-    <field id="SedSi"       long_name="Silicon burial in the sediments"         unit="molSi/m2/s"                      />
-    <field id="SedC"        long_name="Organic C burial in the sediments"       unit="molC/m2/s"                       />
+    <field id="SedCal"      long_name="Calcite burial in the sediments"         unit="molC/m2/s"  grid_ref="grid_T_2D_inner" />
+    <field id="SedSi"       long_name="Silicon burial in the sediments"         unit="molSi/m2/s" grid_ref="grid_T_2D_inner" />
+    <field id="SedC"        long_name="Organic C burial in the sediments"       unit="molC/m2/s"  grid_ref="grid_T_2D_inner" />
     <field id="HYDR"        long_name="Iron input from hydrothemal vents"       unit="mol/m2/s"   grid_ref="grid_T_3D" />
     <field id="EPC100"      long_name="Export of carbon particles at 100 m"     unit="mol/m2/s"                        />
     <field id="EPFE100"     long_name="Export of biogenic iron at 100 m"        unit="mol/m2/s"                        />
@@ -240,10 +240,10 @@
     <field id="EXPSI"       long_name="Export of Silicate"                      unit="mol/m2/s"   grid_ref="grid_T_3D" />
     <field id="EXPCAL"      long_name="Export of Calcite"                       unit="mol/m2/s"   grid_ref="grid_T_3D" />
     <field id="Cflx"        long_name="DIC flux"                                unit="mol/m2/s"                        />
-    <field id="Oflx"        long_name="Oxygen flux"                             unit="mol/m2/s"                        />
-    <field id="Kg"          long_name="Gas transfer"                            unit="mol/m2/s/uatm"                   />
-    <field id="Dpco2"       long_name="Delta CO2"                               unit="uatm"                            />
-    <field id="pCO2sea"     long_name="surface ocean pCO2"                      unit="uatm"                            />
+    <field id="Oflx"        long_name="Oxygen flux"                             unit="mol/m2/s"   grid_ref="grid_T_2D_inner" />
+    <field id="Kg"          long_name="Gas transfer"                            unit="mol/m2/s/uatm"  grid_ref="grid_T_2D_inner" />
+    <field id="Dpco2"       long_name="Delta CO2"                               unit="uatm"       grid_ref="grid_T_2D_inner" />
+    <field id="pCO2sea"     long_name="surface ocean pCO2"                      unit="uatm"       grid_ref="grid_T_2D_inner" />
     <field id="Dpo2"        long_name="Delta O2"                                unit="uatm"                            />
     <field id="Heup"        long_name="Euphotic layer depth"                    unit="m"                               />
     <field id="AtmCo2"      long_name="Atmospheric CO2 concentration"           unit="ppm"                               />
@@ -251,18 +251,18 @@
     <field id="Ironsed"     long_name="Iron deposition from sediment"           unit="mol/m2/s"   grid_ref="grid_T_3D" />
     <field id="FESCAV"      long_name="Scavenging of Iron"                      unit="mmol-Fe/m3/s"   grid_ref="grid_T_3D" />
     <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="LGWCOLL"     long_name="Coagulation loss of ligands"             unit="mmol-L/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="REMINF"      long_name="Oxic remineralization suppy of Fe"       unit="mmol-Fe/m3/s"  grid_ref="grid_T_3D_inner" />
     <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="LPRODR"      long_name="OM remineralisation ligand production rate" unit="nmol-L/m3/s"  grid_ref="grid_T_3D_inner" />
     <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="LIGREM"      long_name="Remineralisation loss of ligands"        unit="nmol-L/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="LIGPR"       long_name="Photochemical loss of ligands"           unit="nmol-L/m3/s"  grid_ref="grid_T_3D_inner" />
     <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_inner" />
-    <field id="LPRODZ"      long_name="microzooplankton ligand production rate" unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
+    <field id="LPRODZ"      long_name="microzooplankton ligand production rate" unit="nmol-L/m3/s"  grid_ref="grid_T_3D_inner" />
     <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" />
 
diff --git a/src/TOP/PISCES/P4Z/p4zfechem.F90 b/src/TOP/PISCES/P4Z/p4zfechem.F90
index 259876648..f98370b0a 100644
--- a/src/TOP/PISCES/P4Z/p4zfechem.F90
+++ b/src/TOP/PISCES/P4Z/p4zfechem.F90
@@ -96,7 +96,7 @@ CONTAINS
       ! This model is based on one ligand, Fe2+ and Fe3+ 
       ! Chemistry is supposed to be fast enough to be at equilibrium
       ! ------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
           zTL1(ji,jj,jk)  = ztotlig(ji,jj,jk)
           zkeq            = fekeq(ji,jj,jk)
           zklight         = 4.77E-7 * etot(ji,jj,jk) * 0.5 / ( 10**(-6.3) )
@@ -134,7 +134,7 @@ CONTAINS
          ENDIF
       ENDIF
 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          ! Scavenging rate of iron. This scavenging rate depends on the load of particles of sea water. 
          ! This parameterization assumes a simple second order kinetics (k[Particles][Fe]).
          ! Scavenging onto dust is also included as evidenced from the DUNE experiments.
@@ -205,14 +205,38 @@ CONTAINS
       !     ---------------------------------
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
          zrfact2 = 1.e3 * rfact2r  ! conversion from mol/L/timestep into mol/m3/s
-         IF( iom_use("Fe3")    )  CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
-         IF( iom_use("FeL1")   )  CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
-         IF( iom_use("TL1")    )  CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
-         IF( iom_use("Totlig") )  CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
-         IF( iom_use("Biron")  )  CALL iom_put("Biron"  , biron  (:,:,:)  * 1e9 * tmask(:,:,:) )   ! biron
-         IF( iom_use("FESCAV") )  CALL iom_put("FESCAV" , zscav3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
-         IF( iom_use("FECOLL") )  CALL iom_put("FECOLL" , zcoll3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
-         IF( iom_use("FEPREC") )  CALL iom_put("FEPREC" , zfeprecip(:,:,:) *1e9*tmask(:,:,:)*zrfact2 )
+         !
+         IF( iom_use ( "Fe3" ) ) THEN   ! Fe3+
+           CALL iom_put( "Fe3",  zFe3(:,:,:) * tmask(:,:,:))
+         ENDIF
+         !
+         IF( iom_use ( "FeL1" ) ) THEN   ! FeL1
+           CALL iom_put( "FeL1",  zFeL1(:,:,:) * tmask(:,:,:))
+         ENDIF
+         !
+         IF( iom_use ( "TL1" ) ) THEN   ! TL1
+           CALL iom_put( "TL1",  zTL1(:,:,:) * tmask(:,:,:))
+         ENDIF
+         !
+         IF( iom_use ( "Totlig" ) ) THEN   ! Totlig
+           CALL iom_put( "Totlig",  ztotlig(:,:,:) * tmask(:,:,:))
+         ENDIF
+         !
+         IF( iom_use ( "Biron" ) ) THEN   ! biron
+           CALL iom_put( "Biron",  biron(:,:,:) * tmask(:,:,:))
+         ENDIF
+         !
+         IF( iom_use ( "FESCAV" ) ) THEN   ! FESCAV
+           CALL iom_put( "FESCAV",  zscav3d(:,:,:) * 1e9 * tmask(:,:,:)  * zrfact2)
+         ENDIF
+         !
+         IF( iom_use ( "FECOLL" ) ) THEN   ! FECOLL
+           CALL iom_put( "FECOLL",  zcoll3d(:,:,:) * 1e9 * tmask(:,:,:)  * zrfact2)
+         ENDIF
+         !
+         IF( iom_use ( "FEPREC" ) ) THEN   ! FEPREC
+           CALL iom_put( "FEPREC",  zfeprecip(:,:,:) * 1e9 * tmask(:,:,:)  * zrfact2)
+         ENDIF
       ENDIF
 
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
diff --git a/src/TOP/PISCES/P4Z/p4zflx.F90 b/src/TOP/PISCES/P4Z/p4zflx.F90
index 1d67f902a..bd632bab8 100644
--- a/src/TOP/PISCES/P4Z/p4zflx.F90
+++ b/src/TOP/PISCES/P4Z/p4zflx.F90
@@ -84,6 +84,7 @@ CONTAINS
       REAL(wp) ::   zyr_dec, zdco2dt
       CHARACTER (len=25) ::   charout
       REAL(wp), DIMENSION(A2D(0)) ::   zkgco2, zkgo2, zh2co3, zoflx,  zpco2atm, zpco2oce  
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zw2d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_flx')
@@ -184,15 +185,50 @@ CONTAINS
       ENDIF
 
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-         CALL iom_put( "AtmCo2"  , satmco2(:,:) * tmask(:,:,1) )   ! Atmospheric CO2 concentration
-         CALL iom_put( "Cflx"    , oce_co2(:,:) * 1000. ) 
-         CALL iom_put( "Oflx"    , zoflx(:,:) * 1000.  )
-         CALL iom_put( "Kg"      , zkgco2(:,:) * tmask(:,:,1)  )
-         CALL iom_put( "Dpco2"   , ( zpco2atm(:,:) - zpco2oce(:,:) ) * tmask(:,:,1) )
-         CALL iom_put( "pCO2sea" , zpco2oce(:,:) * tmask(:,:,1) )
-         CALL iom_put( "Dpo2"    , ( atcox * patm(:,:) - atcox * tr(:,:,1,jpoxy,Kbb) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) )
-         CALL iom_put( "tcflx"   , t_oce_co2_flx     )   ! molC/s
-         CALL iom_put( "tcflxcum", t_oce_co2_flx_cum )   ! molC
+        !
+        ALLOCATE( zw2d(A2D(0)) )  ;  zw2d(A2D(0)) = 0._wp
+        !
+        IF( iom_use ( "AtmCo2" ) )  THEN  ! Atmospheric CO2 concentration
+          CALL iom_put( "AtmCo2"  , satmco2(:,:) * tmask(:,:,1) ) 
+        ENDIF
+        !
+        IF( iom_use ( "Cflx" ) ) THEN 
+          CALL iom_put( "Cflx"    , oce_co2(:,:) * 1000. ) 
+        ENDIF
+        !
+        IF( iom_use ( "Dpo2" ) ) THEN
+          CALL iom_put( "Dpo2"    , ( atcox * patm(:,:) - atcox * tr(:,:,1,jpoxy,Kbb) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) )
+        ENDIF
+        !
+        IF( iom_use ( "tcflx" ) ) THEN  ! molC/s
+          CALL iom_put( "tcflx"   , t_oce_co2_flx ) 
+        ENDIF
+        !
+        IF( iom_use ( "tcflxcum" ) ) THEN ! molC
+          CALL iom_put( "tcflxcum", t_oce_co2_flx_cum )   ! molC
+        ENDIF
+        !
+        IF( iom_use ( "Oflx" ) ) THEN 
+          zw2d(A2D(0)) = zoflx(A2D(0)) * 1000. 
+          CALL iom_put( "Oflx", zw2d )
+        ENDIF
+        !
+        IF( iom_use ( "Kg" ) ) THEN
+          zw2d(A2D(0)) = zkgco2(A2D(0)) * tmask(A2D(0),1)
+          CALL iom_put( "Kg", zw2d )
+        ENDIF
+        !
+        IF( iom_use ( "Dpco2" ) ) THEN
+          zw2d(A2D(0)) =  ( zpco2atm(A2D(0)) - zpco2oce(A2D(0)) ) * tmask(A2D(0),1)
+          CALL iom_put( "Dpco2", zw2d )
+        ENDIF
+        !
+        IF( iom_use ( "pCO2sea" ) ) THEN
+          zw2d(A2D(0)) =  zpco2oce(A2D(0)) * tmask(A2D(0),1)
+          CALL iom_put( "pCO2sea", zw2d )
+        ENDIF
+        !
+        DEALLOCATE( zw2d )
       ENDIF
       !
       IF( ln_timing )   CALL timing_stop('p4z_flx')
diff --git a/src/TOP/PISCES/P4Z/p4zligand.F90 b/src/TOP/PISCES/P4Z/p4zligand.F90
index 07ce54c43..77a6d54d0 100644
--- a/src/TOP/PISCES/P4Z/p4zligand.F90
+++ b/src/TOP/PISCES/P4Z/p4zligand.F90
@@ -47,13 +47,13 @@ CONTAINS
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zlgwp, zlgwpr, zlgwr, zlablgw 
       REAL(wp) ::   zlam1a, zlam1b, zaggliga, zligco
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zligrem, zligpr, zligprod, zlcoll3d
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zligrem, zligpr, zligprod, zlcoll3d
       CHARACTER (len=25) ::   charout
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_ligand')
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          !
          ! ------------------------------------------------------------------
          ! Remineralization of iron ligands
@@ -89,16 +89,16 @@ CONTAINS
       !     ---------------------------------
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
          IF( iom_use( "LIGREM" ) ) THEN
-           zligrem(:,:,jpk) = 0.  ; CALL iom_put( "LIGREM", zligrem(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
+           zligrem(A2D(0),jpk) = 0.  ; CALL iom_put( "LIGREM", zligrem(A2D(0),:) * 1e9 * 1.e+3 * rfact2r * tmask(A2D(0),:) )
          ENDIF
          IF( iom_use( "LIGPR" ) ) THEN
-           zligpr(:,:,jpk) = 0.   ; CALL iom_put( "LIGPR" , zligpr(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
+           zligpr(A2D(0),jpk) = 0.   ; CALL iom_put( "LIGPR" , zligpr(A2D(0),:) * 1e9 * 1.e+3 * rfact2r * tmask(A2D(0),:) )
          ENDIF
          IF( iom_use( "LPRODR" ) ) THEN
-           zligprod(:,:,jpk) = 0. ; CALL iom_put( "LPRODR", zligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
+           zligprod(A2D(0),jpk) = 0. ; CALL iom_put( "LPRODR", zligprod(A2D(0),:) * 1e9 * 1.e+3 * rfact2r * tmask(A2D(0),:) )
          ENDIF
          IF( iom_use( "LGWCOLL" ) ) THEN
-            zlcoll3d(:,:,jpk) = 0. ; CALL iom_put( "LGWCOLL", zlcoll3d(:,:,:) * 1.e9 *  1.e+3 * rfact2r * tmask(:,:,:) )
+            zlcoll3d(A2D(0),jpk) = 0. ; CALL iom_put( "LGWCOLL", zlcoll3d(A2D(0),:) * 1.e9 *  1.e+3 * rfact2r * tmask(A2D(0),:) )
          ENDIF
       ENDIF
       !
diff --git a/src/TOP/PISCES/P4Z/p4zmeso.F90 b/src/TOP/PISCES/P4Z/p4zmeso.F90
index f9f1f149d..8dc46de78 100644
--- a/src/TOP/PISCES/P4Z/p4zmeso.F90
+++ b/src/TOP/PISCES/P4Z/p4zmeso.F90
@@ -428,8 +428,8 @@ CONTAINS
         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 )
+           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 )
diff --git a/src/TOP/PISCES/P4Z/p4zmicro.F90 b/src/TOP/PISCES/P4Z/p4zmicro.F90
index 712554553..9b07f7453 100644
--- a/src/TOP/PISCES/P4Z/p4zmicro.F90
+++ b/src/TOP/PISCES/P4Z/p4zmicro.F90
@@ -79,7 +79,6 @@ CONTAINS
       REAL(wp) :: zgraznc, zgrazpoc, zgrazdc, zgrazpof, zgrazdf, zgraznf
       REAL(wp) :: zsigma, zdiffdn, ztmp1, ztmp2, ztmp3, ztmptot, zproport
       REAL(wp), DIMENSION(A2D(0),jpk) :: zgrazing, zfezoo
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zzligprod
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zw3d
       CHARACTER (len=25) :: charout
 
@@ -87,11 +86,6 @@ CONTAINS
       !
       IF( ln_timing )   CALL timing_start('p4z_micro')
       !
-      IF (ln_ligand) THEN
-         ALLOCATE( zzligprod(jpi,jpj,jpk) )
-         zzligprod(:,:,:) = 0._wp
-      ENDIF
-      !
       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
@@ -216,7 +210,6 @@ CONTAINS
          !
          IF( ln_ligand ) THEN
             tr(ji,jj,jk,jplgw,Krhs) = tr(ji,jj,jk,jplgw,Krhs) + (zgrarem - zgrarsig) * ldocz
-            zzligprod(ji,jj,jk) = (zgrarem - zgrarsig) * ldocz
          ENDIF
          !
          tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) - o2ut * zgrarsig
@@ -271,8 +264,9 @@ CONTAINS
           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(:,:,:))
+        IF( ln_ligand .AND. iom_use( "LPRODZ" ) ) THEN
+            zw3d(A2D(0),1:jpkm1) =  (zgrarem - zgrarsig) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
+           CALL iom_put( "LPRODZ", zw3d  * ldocz * 1e9 )
         ENDIF
         !
         DEALLOCATE( zw3d )
diff --git a/src/TOP/PISCES/P4Z/p4zopt.F90 b/src/TOP/PISCES/P4Z/p4zopt.F90
index 8daf2edf7..94ec459ad 100644
--- a/src/TOP/PISCES/P4Z/p4zopt.F90
+++ b/src/TOP/PISCES/P4Z/p4zopt.F90
@@ -88,7 +88,7 @@ CONTAINS
       !
       ! Computation of the light attenuation parameters based on a 
       ! look-up table
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          zchl =  ( tr(ji,jj,jk,jpnch,Kbb) + tr(ji,jj,jk,jpdch,Kbb) + rtrn ) * 1.e6
          IF( ln_p5z )   zchl = zchl + tr(ji,jj,jk,jppch,Kbb) * 1.e6
          zchl = MIN(  10. , MAX( 0.05, zchl )  )
@@ -240,7 +240,7 @@ CONTAINS
       heup   (:,:) = gdepw(:,:,2,Kmm)
       heup_01(:,:) = gdepw(:,:,2,Kmm)
 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr)
+      DO_3D( 0, 0, 0, 0, 2, nksr)
         IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >=  zqsr100(ji,jj) )  THEN
            neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
            !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
@@ -261,7 +261,7 @@ CONTAINS
       zetmp1 (:,:)   = 0.e0
       zetmp2 (:,:)   = 0.e0
 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
+      DO_3D( 0, 0, 0, 0, 1, nksr)
          IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
             zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Actual PAR for remineralisation
             zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Par averaged over 24h for production
@@ -272,7 +272,7 @@ CONTAINS
       emoy(:,:,:) = etot(:,:,:)       ! remineralisation
       zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle 
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
+      DO_3D( 0, 0, 0, 0, 1, nksr)
          IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
             z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
             emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep
@@ -286,7 +286,7 @@ CONTAINS
       zetmp3 (:,:)   = 0.e0
       zetmp4 (:,:)   = 0.e0
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
+      DO_3D( 0, 0, 0, 0, 1, nksr)
          IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
             zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Nanophytoplankton
             zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! Diatoms
@@ -296,7 +296,7 @@ CONTAINS
       enanom(:,:,:) = enano(:,:,:)
       ediatm(:,:,:) = ediat(:,:,:)
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
+      DO_3D( 0, 0, 0, 0, 1, nksr)
          IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
             z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
             enanom(ji,jj,jk) = zetmp3(ji,jj) * z1_dep
@@ -307,7 +307,7 @@ CONTAINS
       IF( ln_p5z ) THEN
          ! Picophytoplankton when using PISCES-QUOTA
          ALLOCATE( zetmp5(jpi,jpj) )  ;   zetmp5 (:,:) = 0.e0
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
+         DO_3D( 0, 0, 0, 0, 1, nksr)
             IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
                zetmp5(ji,jj)  = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
             ENDIF
@@ -315,7 +315,7 @@ CONTAINS
          !
          epicom(:,:,:) = epico(:,:,:)
          !
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr)
+         DO_3D( 0, 0, 0, 0, 1, nksr)
             IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
                z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
                epicom(ji,jj,jk) = zetmp5(ji,jj) * z1_dep
@@ -371,7 +371,7 @@ CONTAINS
          pe2(:,:,1) = zqsr(:,:)
          pe3(:,:,1) = zqsr(:,:)
          !
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr + 1)
+         DO_3D( 0, 0, 0, 0, 2, nksr + 1)
             pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -e3t(ji,jj,jk-1,Kmm) * xsi0r )
             pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -ekb  (ji,jj,jk-1 )        )
             pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -ekg  (ji,jj,jk-1 )        )
@@ -384,7 +384,7 @@ CONTAINS
         pe2(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekg(:,:,1) )
         pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) )
         !
-        DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr)
+        DO_3D( 0, 0, 0, 0, 2, nksr)
            pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) )
            pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) )
            pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( ekr(ji,jj,jk-1) + ekr(ji,jj,jk) ) )
diff --git a/src/TOP/PISCES/P4Z/p4zpoc.F90 b/src/TOP/PISCES/P4Z/p4zpoc.F90
index d729eff71..ab156817e 100644
--- a/src/TOP/PISCES/P4Z/p4zpoc.F90
+++ b/src/TOP/PISCES/P4Z/p4zpoc.F90
@@ -77,6 +77,7 @@ CONTAINS
       REAL(wp), DIMENSION(A2D(0)  )   :: totprod, totthick, totcons 
       REAL(wp), DIMENSION(A2D(0),jpk)   :: zremipoc, zremigoc, zorem3, ztremint, zfolimi
       REAL(wp), DIMENSION(A2D(0),jpk,jcpoc) :: alphag
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )  CALL timing_start('p4z_poc')
@@ -443,9 +444,25 @@ CONTAINS
      IF( lk_iomput ) THEN
         IF( knt == nrdttrc ) THEN
           zrfact2 = 1.e3 * rfact2r
-          CALL iom_put( "REMINP" , zremipoc(:,:,:) * tmask(:,:,:) )  ! Remineralisation rate of small particles
-          CALL iom_put( "REMING" , zremigoc(:,:,:) * tmask(:,:,:) )  ! Remineralisation rate of large particles
-          CALL iom_put( "REMINF" , zfolimi(:,:,:)  * tmask(:,:,:)  * 1.e+9 * zrfact2 )  ! Remineralisation of biogenic particulate iron
+          !
+          ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+          !
+          IF( iom_use ( "REMINP" ) ) THEN   ! Remineralisation rate of small particles
+            zw3d(A2D(0),1:jpkm1) = zremipoc(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+            CALL iom_put( "REMINP", zw3d )
+          ENDIF
+          !
+          IF( iom_use ( "REMING" ) ) THEN   ! Remineralisation rate of large particles
+            zw3d(A2D(0),1:jpkm1) = zremigoc(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+            CALL iom_put( "REMING", zw3d )
+          ENDIF
+          !
+          IF( iom_use ( "REMINF" ) ) THEN   ! Remineralisation of biogenic particulate iron
+            zw3d(A2D(0),1:jpkm1) = zfolimi(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1) * 1.e+9 * zrfact2
+            CALL iom_put( "REMINF", zw3d )
+          ENDIF
+          !
+          DEALLOCATE ( zw3d )  
         ENDIF
      ENDIF
 
diff --git a/src/TOP/PISCES/P4Z/p4zprod.F90 b/src/TOP/PISCES/P4Z/p4zprod.F90
index a48a6242a..f3189c38d 100644
--- a/src/TOP/PISCES/P4Z/p4zprod.F90
+++ b/src/TOP/PISCES/P4Z/p4zprod.F90
@@ -83,7 +83,7 @@ CONTAINS
       REAL(wp), DIMENSION(A2D(0),jpk) :: zprorcan, zprorcad, zprofed, zprofen
       REAL(wp), DIMENSION(A2D(0),jpk) :: zpronewn, zpronewd
       REAL(wp), DIMENSION(A2D(0),jpk) :: zmxl_fac, zmxl_chl
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zpligprod, zw3d
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_prod')
@@ -424,7 +424,7 @@ CONTAINS
        ENDIF
        !
        IF( ln_ligand .AND. ( iom_use( "LPRODP" ) .OR. iom_use( "LDETP" ) ) ) THEN
-           zw3d(A2D(0),1:jpkm1) = excretd * zprorcad(A2D(0),1:jpkm1) + excretn * zprorcan(A2D(0),1:jpkm1) &
+           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 )
            !
diff --git a/src/TOP/PISCES/P4Z/p4zsed.F90 b/src/TOP/PISCES/P4Z/p4zsed.F90
index ca6bdb281..040b05059 100644
--- a/src/TOP/PISCES/P4Z/p4zsed.F90
+++ b/src/TOP/PISCES/P4Z/p4zsed.F90
@@ -74,6 +74,7 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj    ) :: zsedcal, zsedsi, zsedc
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsoufer, zlight
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zpdep
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zw2d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )  CALL timing_start('p4z_sed')
@@ -92,7 +93,7 @@ CONTAINS
 
       ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments
       ! --------------------------------------------------------------------
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          ikt  = mbkt(ji,jj)
          zdep = e3t(ji,jj,ikt,Kmm) / xstep
          zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) )
@@ -103,7 +104,7 @@ CONTAINS
          ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used
          ! Computation of the fraction of organic matter that is permanently buried from Dunne's model
          ! -------------------------------------------------------
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
            IF( tmask(ji,jj,1) == 1 ) THEN
               ikt = mbkt(ji,jj)
               zflx = (  tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj)   &
@@ -129,7 +130,7 @@ CONTAINS
       ! ------------------------------------------------------
       IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac
 
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          ikt  = mbkt(ji,jj)
          zdep = xstep / e3t(ji,jj,ikt,Kmm) 
          zwsc = zwsbio4(ji,jj) * zdep
@@ -141,7 +142,7 @@ CONTAINS
       END_2D
       !
       IF( .NOT.lk_sed ) THEN
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
             ikt  = mbkt(ji,jj)
             zdep = xstep / e3t(ji,jj,ikt,Kmm) 
             zwsc = zwsbio4(ji,jj) * zdep
@@ -159,7 +160,7 @@ CONTAINS
          END_2D
       ENDIF
       !
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          ikt  = mbkt(ji,jj)
          zdep = xstep / e3t(ji,jj,ikt,Kmm) 
          zws4 = zwsbio4(ji,jj) * zdep
@@ -171,7 +172,7 @@ CONTAINS
       END_2D
       !
       IF( ln_p5z ) THEN
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
             ikt  = mbkt(ji,jj)
             zdep = xstep / e3t(ji,jj,ikt,Kmm) 
             zws4 = zwsbio4(ji,jj) * zdep
@@ -186,7 +187,7 @@ CONTAINS
       IF( .NOT.lk_sed ) THEN
          ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after
          ! denitrification in the sediments. Not very clever, but simpliest option.
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
             ikt  = mbkt(ji,jj)
             zdep = xstep / e3t(ji,jj,ikt,Kmm) 
             zws4 = zwsbio4(ji,jj) * zdep
@@ -223,7 +224,7 @@ CONTAINS
          zsoufer(:,:,jk) = zlight(:,:,jk) * 2E-11 / ( 2E-11 + biron(:,:,jk) )
       ENDDO
       IF( ln_p4z ) THEN
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
             !                      ! Potential nitrogen fixation dependant on temperature and iron
             ztemp = ts(ji,jj,jk,jp_tem,Kmm)
             zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) / rno3
@@ -239,7 +240,7 @@ CONTAINS
             nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk)
          END_3D
       ELSE       ! p5z
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
             !                      ! Potential nitrogen fixation dependant on temperature and iron
             ztemp = ts(ji,jj,jk,jp_tem,Kmm)
             zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
@@ -260,7 +261,7 @@ CONTAINS
       ! Nitrogen change due to nitrogen fixation
       ! ----------------------------------------
       IF( ln_p4z ) THEN
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
             zfact = nitrpot(ji,jj,jk) * nitrfix
             tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0
             tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0
@@ -278,7 +279,7 @@ CONTAINS
             &                     * 0.001 * tr(ji,jj,jk,jpdoc,Kbb) * xstep
          END_3D
       ELSE    ! p5z
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
             zfact = nitrpot(ji,jj,jk) * nitrfix
             tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0
             tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0
@@ -306,12 +307,34 @@ CONTAINS
       ENDIF
 
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-         zfact = 1.e+3 * rfact2r !  conversion from molC/l/kt  to molN/m3/s
-         CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation 
-         CALL iom_put( "SedCal", zsedcal(:,:) * zfact )
-         CALL iom_put( "SedSi" , zsedsi (:,:) * zfact )
-         CALL iom_put( "SedC"  , zsedc  (:,:) * zfact )
-         CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 )
+          zfact = 1.e+3 * rfact2r !  conversion from molC/l/kt  to molN/m3/s
+          !
+          ALLOCATE( zw2d(jpi,jpj) )  ;  zw2d(:,:) = 0._wp
+          !
+          IF( iom_use ( "Nfix" ) ) THEN   ! nitrogen fixation
+            CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )
+          ENDIF
+          !
+          IF( iom_use ( "Sdenit" ) ) THEN
+            CALL iom_put( "SedC", sdenit (:,:) * zfact * rno3 )
+          ENDIF
+          !
+          IF( iom_use ( "SedCal" ) ) THEN 
+            zw2d(:,:) = zsedcal(:,:) * zfact
+            CALL iom_put( "SedCal", zw2d )
+          ENDIF
+          !
+          IF( iom_use ( "SedSi" ) ) THEN
+            zw2d(:,:) = zsedsi(:,:) * zfact
+            CALL iom_put( "SedSi", zw2d )
+          ENDIF
+          !
+          IF( iom_use ( "SedC" ) ) THEN
+            zw2d(:,:) = zsedc(:,:) * zfact
+            CALL iom_put( "SedC", zw2d )
+          ENDIF
+          !
+          DEALLOCATE( zw2d )
       ENDIF
       !
       IF(sn_cfctl%l_prttrc) THEN  ! print mean trneds (USEd for debugging)
diff --git a/src/TOP/PISCES/P4Z/p4zsink.F90 b/src/TOP/PISCES/P4Z/p4zsink.F90
index e900c82ca..48bcf3f93 100644
--- a/src/TOP/PISCES/P4Z/p4zsink.F90
+++ b/src/TOP/PISCES/P4Z/p4zsink.F90
@@ -135,15 +135,41 @@ CONTAINS
      IF( lk_iomput .AND.  knt == nrdttrc ) THEN
        zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
        !
-       CALL iom_put( "EPC100"  , ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ) ! Export of carbon at 100m 
-       CALL iom_put( "EPFE100" , ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ) ! Export of iron at 100m 
-       CALL iom_put( "EPCAL100", sinkcal(:,:,ik100) * zfact * tmask(:,:,1) )      ! Export of calcite at 100m 
-       CALL iom_put( "EPSI100" , sinksil(:,:,ik100) * zfact * tmask(:,:,1) )          ! Export of bigenic silica at 100m 
-       CALL iom_put( "EXPC"    , ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ) ! Export of carbon in the water column 
-       CALL iom_put( "EXPFE"   , ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ) ! Export of iron  
-       CALL iom_put( "EXPCAL"  , sinkcal(:,:,:) * zfact * tmask(:,:,:) )      ! Export of calcite 
-       CALL iom_put( "EXPSI"   , sinksil(:,:,:) * zfact * tmask(:,:,:) )      ! Export of bigenic silica
-       CALL iom_put( "tcexp"   , t_oce_co2_exp * zfact )   ! molC/s
+       IF( iom_use ( "EPC100" ) ) THEN ! Export of carbon at 100m
+         CALL iom_put( "EPC100"  , ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ) 
+       ENDIF
+       !
+       IF( iom_use ( "EPFE100" ) ) THEN ! Export of iron at 100m
+         CALL iom_put( "EPFE100" , ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ) 
+       ENDIF
+       !
+       IF( iom_use ( "EPCAL100" ) ) THEN ! Export of calcite at 100m
+         CALL iom_put( "EPCAL100", sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ) 
+       ENDIF
+       !
+       IF( iom_use ( "EPSI100" ) ) THEN ! Export of bigenic silica at 100m
+         CALL iom_put( "EPSI100" , sinksil(:,:,ik100) * zfact * tmask(:,:,1) ) 
+       ENDIF
+       !
+       IF( iom_use ( "EXPC" ) ) THEN ! Export of carbon in the water column
+         CALL iom_put( "EXPC"    , ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ) 
+       ENDIF
+       !
+       IF( iom_use ( "EXPFE" ) ) THEN ! Export of iron
+         CALL iom_put( "EXPFE"   , ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ) 
+       ENDIF
+       !
+       IF( iom_use ( "EXPCAL" ) ) THEN ! Export of calcite
+         CALL iom_put( "EXPCAL"  , sinkcal(:,:,:) * zfact * tmask(:,:,:) ) 
+       ENDIF
+       !
+       IF( iom_use ( "EXPSI" ) ) THEN ! Export of bigenic silica
+         CALL iom_put( "EXPSI"   , sinksil(:,:,:) * zfact * tmask(:,:,:) ) 
+       ENDIF
+       !
+       IF( iom_use ( "EXPSI" ) ) THEN ! molC/s
+         CALL iom_put( "tcexp"   , t_oce_co2_exp * zfact ) 
+       ENDIF
        ! 
       ENDIF
       !
diff --git a/src/TOP/PISCES/P4Z/p4zsms.F90 b/src/TOP/PISCES/P4Z/p4zsms.F90
index 664d12064..4f1c7122b 100644
--- a/src/TOP/PISCES/P4Z/p4zsms.F90
+++ b/src/TOP/PISCES/P4Z/p4zsms.F90
@@ -140,7 +140,7 @@ CONTAINS
          ! ------------------------------------------------------------------
          xnegtr(:,:,:) = 1.e0
          DO jn = jp_pcs0, jp_pcs1
-            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk)
+            DO_3D( 0, 0, 0, 0, 1, jpk)
                IF( ( tr(ji,jj,jk,jn,Kbb) + tr(ji,jj,jk,jn,Krhs) ) < 0.e0 ) THEN
                   ztra             = ABS( tr(ji,jj,jk,jn,Kbb) ) / ( ABS( tr(ji,jj,jk,jn,Krhs) ) + rtrn )
                   xnegtr(ji,jj,jk) = MIN( xnegtr(ji,jj,jk),  ztra )
diff --git a/src/TOP/PISCES/P4Z/p5zlim.F90 b/src/TOP/PISCES/P4Z/p5zlim.F90
index 91e94183a..dffbe0384 100644
--- a/src/TOP/PISCES/P4Z/p5zlim.F90
+++ b/src/TOP/PISCES/P4Z/p5zlim.F90
@@ -136,7 +136,8 @@ CONTAINS
       REAL(wp) ::   zfvn, zfvp, zfvf, zsizen, zsizep, zsized, znanochl, zpicochl, zdiatchl
       REAL(wp) ::   zqfemn, zqfemp, zqfemd, zbactno3, zbactnh4, zbiron
       REAL(wp) ::   znutlimtot, zlimno3, zlimnh4, zlim1f, zsizetmp
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrassn, zrassp, zrassd
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zrassn, zrassp, zrassd
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p5z_lim')
@@ -144,7 +145,7 @@ CONTAINS
       zratchl = 6.0
       sizena(:,:,:) = 0.0  ;  sizepa(:,:,:) = 0.0  ;  sizeda(:,:,:) = 0.0
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          ! Computation of the Chl/C ratio of each phytoplankton group
          ! -------------------------------------------------------
          z1_trnphy   = 1. / ( tr(ji,jj,jk,jpphy,Kbb) + rtrn )
@@ -407,7 +408,7 @@ CONTAINS
       ! nutrient uptake pool and assembly machinery. DNA is assumed to represent 1% of the dry mass of 
       ! phytoplankton (see Daines et al., 2013). 
       ! --------------------------------------------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          ! Size estimation of nanophytoplankton based on total biomass
          ! Assumes that larger biomass implies addition of larger cells
          ! ------------------------------------------------------------
@@ -465,7 +466,7 @@ CONTAINS
       ! This is a purely adhoc formulation described in Aumont et al. (2015)
       ! This fraction depends on nutrient limitation, light, temperature
       ! --------------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          zlim1 =  tr(ji,jj,jk,jpnh4,Kbb) / ( tr(ji,jj,jk,jpnh4,Kbb) + concnnh4 ) + tr(ji,jj,jk,jpno3,Kbb)    &
          &        / ( tr(ji,jj,jk,jpno3,Kbb) + concnno3 ) * ( 1.0 - tr(ji,jj,jk,jpnh4,Kbb)   &
          &        / ( tr(ji,jj,jk,jpnh4,Kbb) + concnnh4 ) )
@@ -482,7 +483,7 @@ 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)
+      DO_3D( 0, 0, 0, 0, 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) )  )
@@ -490,19 +491,75 @@ CONTAINS
       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
-        CALL iom_put( "LPnut"  , xlimpic(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
-        CALL iom_put( "LDnut"  , xlimdia(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
-        CALL iom_put( "LNFe"   , xlimnfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "LPFe"   , xlimpfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "LDFe"   , xlimdfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "SIZEN"  , sizen  (:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "SIZEP"  , sizep  (:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "SIZED"  , sized  (:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "RASSN"  , zrassn (:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "RASSP"  , zrassp (:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "RASSD"  , zrassd (:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        !
+        ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+        !
+        IF( iom_use ( "xfracal" ) ) THEN   ! euphotic layer deptht
+          zw3d(A2D(0),1:jpkm1) = xfracal(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "xfracal",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "LNnut" ) ) THEN   ! Nutrient limitation term
+          zw3d(A2D(0),1:jpkm1) = xlimphy(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "LNnut",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "LPnut" ) ) THEN   ! Nutrient limitation term
+          zw3d(A2D(0),1:jpkm1) = xlimpic(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "LPnut",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "LDnut" ) ) THEN   ! Nutrient limitation term
+          zw3d(A2D(0),1:jpkm1) = xlimdia(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "LDnut",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "LNFe" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = xlimnfe(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "LNFe",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "LPFe" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = xlimpfe(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "LPFe",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "LDFe" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = xlimdfe(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "LDFe",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "SIZEN" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = sizen(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "SIZEN",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "SIZEP" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = sizep(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "SIZEP",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "SIZED" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = sized(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "SIZED",  zw3d)
+        ENDIF
+       !
+        IF( iom_use ( "RASSN" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = zrassn(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "RASSN",  zw3d)
+        ENDIF
+       !
+        IF( iom_use ( "RASSP" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = zrassp(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "RASSP",  zw3d)
+        ENDIF
+       !
+        IF( iom_use ( "RASSD" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = zrassd(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "RASSD",  zw3d)
+        ENDIF
+        !
+        DEALLOCATE( zw3d )
       ENDIF
       !
       IF( ln_timing )  CALL timing_stop('p5z_lim')
diff --git a/src/TOP/PISCES/P4Z/p5zmeso.F90 b/src/TOP/PISCES/P4Z/p5zmeso.F90
index 24ca0260d..7dc423073 100644
--- a/src/TOP/PISCES/P4Z/p5zmeso.F90
+++ b/src/TOP/PISCES/P4Z/p5zmeso.F90
@@ -103,14 +103,14 @@ CONTAINS
       REAL(wp) :: zmigreltime, zrum, zcodel, zargu, zval, zmigthick 
       CHARACTER (len=25) :: charout
       REAL(wp) :: zrfact2, zmetexcess, zsigma, zdiffdn
-      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(A2D(0),jpk) :: zgrazing, zfezoo2
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zgrarem, zgraref, zgrapoc, zgrapof
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zgrarep, zgraren, zgrapon, zgrapop
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zgradoc, zgradon, zgradop
       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(:,:,:) ::   zw3d
 
       !!---------------------------------------------------------------------
       !
@@ -136,7 +136,7 @@ CONTAINS
       zmetexcess = 0.0
       IF ( bmetexc2 ) zmetexcess = 1.0
 
-      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
 
@@ -427,7 +427,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, jpkm1)
+          DO_3D( 0, 0, 0, 0, 1, jpkm1)
              zmigreltime = (1. - strn(ji,jj))
              IF( gdept(ji,jj,jk,Kmm) <= heup(ji,jj) ) THEN
                 zmigthick  = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * ( 1. - zmigreltime ) 
@@ -460,7 +460,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)
@@ -490,7 +490,7 @@ CONTAINS
         !   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, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zgrarep(ji,jj,jk) 
          tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zgraren(ji,jj,jk)
          tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zgradoc(ji,jj,jk)
@@ -513,11 +513,30 @@ CONTAINS
       END_3D
       !
       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", zgradoc(:,:,:) * 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) = zgradoc(A2D(0),1:jpkm1) * ldocz * 1e9 * 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)
@@ -626,7 +645,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)
@@ -639,7 +658,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
@@ -648,7 +667,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,jppch,Kbb) + 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)
@@ -657,7 +676,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
@@ -669,7 +688,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/p5zmicro.F90 b/src/TOP/PISCES/P4Z/p5zmicro.F90
index 6edfde55e..e67b50197 100644
--- a/src/TOP/PISCES/P4Z/p5zmicro.F90
+++ b/src/TOP/PISCES/P4Z/p5zmicro.F90
@@ -88,8 +88,9 @@ CONTAINS
       REAL(wp) :: zgraznc, zgraznn, zgraznp, zgrazpoc, zgrazpon, zgrazpop, zgrazpof
       REAL(wp) :: zgrazdc, zgrazdn, zgrazdp, zgrazdf, zgraznf, zgrazz
       REAL(wp) :: zgrazpc, zgrazpn, zgrazpp, zgrazpf, zbeta, zrfact2, zmetexcess
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrazing, zfezoo, zzligprod
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zgrazing, zfezoo
       REAL(wp) :: zsigma, zdiffdn, zdiffpn, zdiffdp, zproport, zproport2
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zw3d
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
       !
@@ -99,7 +100,7 @@ CONTAINS
       zmetexcess = 0.0
       IF ( bmetexc ) zmetexcess = 1.0
       !
-      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
          ! Proportion of nano and diatoms that are within the size range
@@ -298,7 +299,6 @@ CONTAINS
          !
          IF( ln_ligand ) THEN 
             tr(ji,jj,jk,jplgw,Krhs) = tr(ji,jj,jk,jplgw,Krhs) + zgradoc * ldocz
-            zzligprod(ji,jj,jk) = zgradoc * ldocz
          ENDIF
          !
          tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zgradon
@@ -342,15 +342,25 @@ 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 ( "LPRODZ" ) ) THEN
+           zw3d(A2D(0),1:jpkm1) = zgradoc * ldocz  * 1e9 * 1.e+3 * rfact2r * tmask(A2D(0),1:jpkm1)
+           CALL iom_put( "LPRODZ", zw3d )
+        ENDIF
+        !
+        DEALLOCATE( zw3d )
       ENDIF
       !
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
diff --git a/src/TOP/PISCES/P4Z/p5zprod.F90 b/src/TOP/PISCES/P4Z/p5zprod.F90
index 848f85505..d8cf514b3 100644
--- a/src/TOP/PISCES/P4Z/p5zprod.F90
+++ b/src/TOP/PISCES/P4Z/p5zprod.F90
@@ -85,20 +85,20 @@ CONTAINS
       REAL(wp) ::   zqfpmax, zqfnmax, zqfdmax
       REAL(wp) ::   zfact, zrfact2, zmaxsi, zratiosi, zsizetmp, zlimfac, zsilim
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zmixnano, zmixpico, zmixdiat
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadp, zpislopeadd
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprnut, zprmaxp, zprmaxn, zprmaxd
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprbio, zprpic, zprdia, zysopt
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprchln, zprchlp, zprchld
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprorcan, zprorcap, zprorcad 
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprofed, zprofep, zprofen
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpronewn, zpronewp, zpronewd
-      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
+!      REAL(wp), DIMENSION(A2D(0)    ) :: zmixnano, zmixpico, zmixdiat
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zpislopeadn, zpislopeadp, zpislopeadd
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zprnut, zprmaxp, zprmaxn, zprmaxd
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zprbio, zprpic, zprdia, zysopt
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zprchln, zprchlp, zprchld
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zprorcan, zprorcap, zprorcad 
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zprofed, zprofep, zprofen
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zpronewn, zpronewp, zpronewd
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zproregn, zproregp, zproregd
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zpropo4n, zpropo4p, zpropo4d
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zprodopn, zprodopp, zprodopd
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zrespn, zrespp, zrespd
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zmxl_fac, zmxl_chl
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p5z_prod')
@@ -132,7 +132,7 @@ CONTAINS
 
       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 ))
@@ -144,7 +144,7 @@ CONTAINS
 
       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
@@ -166,7 +166,7 @@ CONTAINS
 
       ! Computation of the P-I slope for nanos, picos and diatoms
       ! The formulation proposed by Geider et al. (1997) has been used.
-      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
             ! Computation of the P-I slope for nanos and diatoms
             ztn         = MAX( 0., ts(ji,jj,jk,jp_tem,Kmm) - 15. )
@@ -209,7 +209,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)
           IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             ! Si/C of diatoms
             ! ------------------------
@@ -242,7 +242,7 @@ CONTAINS
       !  Sea-ice effect on production
       ! No production is assumed below sea ice
       ! -------------------------------------- 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 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) ) 
@@ -256,7 +256,7 @@ CONTAINS
       ! quota, uptake is downregulated according to a sigmoidal function 
       ! (power 2), as proposed by Flynn (2003)
       ! ---------------------------------------------------------------------------
-      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
             !  production terms for nanophyto.
             zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * tr(ji,jj,jk,jpphy,Kbb) * rfact2
@@ -307,7 +307,7 @@ CONTAINS
       ! quota, uptake is downregulated according to a sigmoidal function 
       ! (power 2), as proposed by Flynn (2003)
       ! ---------------------------------------------------------------------------
-      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
             !  production terms for picophyto.
             zprorcap(ji,jj,jk) = zprpic(ji,jj,jk)  * xlimpic(ji,jj,jk) * tr(ji,jj,jk,jppic,Kbb) * rfact2
@@ -357,7 +357,7 @@ CONTAINS
       ! quota, uptake is downregulated according to a sigmoidal function 
       ! (power 2), as proposed by Flynn (2003)
       ! ---------------------------------------------------------------------------
-      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
             !  production terms for diatomees
             zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * tr(ji,jj,jk,jpdia,Kbb) * rfact2
@@ -403,7 +403,7 @@ CONTAINS
       ! Production of Chlorophyll. The formulation proposed by Geider et al.
       ! is adopted here.
       ! --------------------------------------------------------------------
-      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
                !  production terms for nanophyto. ( chlorophyll )
             znanotot = enanom(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
@@ -428,7 +428,7 @@ CONTAINS
       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)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
         zpptot   = 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)
@@ -513,7 +513,7 @@ CONTAINS
      ! Shaked and Lis (2012)
      ! -------------------------------------------------------------------------
      IF( ln_ligand ) THEN
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
            zproddoc = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk) + excretp * zprorcap(ji,jj,jk)
            zprodfer = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) + texcretp * zprofep(ji,jj,jk)
            zprodlig = plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * lthet 
@@ -522,6 +522,7 @@ 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( 'p5zprod', ( zprorcan(:,:,:) + zprorcad(:,:,:) + zprorcap(:,:,:) ) * cvol(:,:,:) )
@@ -529,40 +530,134 @@ CONTAINS
     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( "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
-       CALL iom_put( "PBSi"    , zprmaxd (:,:,:) * zfact * tmask(:,:,:) * zysopt(:,:,:)  ) ! biogenic silica production
-       CALL iom_put( "PFeP"    , zprofep (:,:,:) * zfact * tmask(:,:,:)  ) ! biogenic iron production by picophyto
-       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 ( "PPPHYP" ) ) THEN   ! primary production by picophyto
+         zw3d(A2D(0),1:jpkm1) = zprorcap(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "PPPHYP", zw3d )
+       ENDIF
+       !
+       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 diatoms
+         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 ( "PPNEWP" ) ) THEN   ! new primary production by picophyto
+         zw3d(A2D(0),1:jpkm1) = zpronewp(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "PPNEWP", zw3d )
+       ENDIF
+       !
+       IF( iom_use ( "PPNEWN" ) ) THEN   ! new primary production by nanophyto
+         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 diatoms
+         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) = zprmaxd(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1) * zysopt(A2D(0),1:jpkm1)
+         CALL iom_put( "PBSi", zw3d )
+       ENDIF
+       !
+       IF( iom_use ( "PFeP" ) ) THEN   ! biogenic iron production by picophyto
+         zw3d(A2D(0),1:jpkm1) = zprofep(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1) 
+         CALL iom_put( "PFeP", 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 diatoms
+         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(:,:,:) + excretp * zprorcap(:,:,:)
-           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) &
+             &                   + excretp * zprorcap(A2D(0),1:jpkm1) ) * zfact * tmask(A2D(0),1:jpkm1)
+           CALL iom_put( "LPRODP"  , zw3d * ldocp * 1e9 )
            !
-           zpligprod(:,:,:) = ( texcretn * zprofen(:,:,:) + texcretd * zprofed(:,:,:) + texcretp * zprofep(:,:,:) ) & 
-             &                  * 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) &
+             &                  + texcretp * zprofep(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) = zprmaxn(A2D(0),1:jpkm1)  * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "Mumax", zw3d )
+       ENDIF
+       !
+       IF( iom_use ( "MuP" ) ) THEN    ! Realized growth rate for picophyto
+         zw3d(A2D(0),1:jpkm1) = zprpic(A2D(0),1:jpkm1) * xlimpic(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "MuP", zw3d )
+       ENDIF
+       !
+       !
+       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 ( "LPlight" ) ) THEN    ! light limitation term for pico
+         zw3d(A2D(0),1:jpkm1) = zprpic(A2D(0),1:jpkm1) / ( zprmaxp(A2D(0),1:jpkm1) + rtrn ) * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "LPlight", zw3d )
+       ENDIF
+       !
+       IF( iom_use ( "LNlight" ) ) THEN    ! light limitation term for nano
+         zw3d(A2D(0),1:jpkm1) = zprbio(A2D(0),1:jpkm1) / ( zprmaxn(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 diatoms
+         zw3d(A2D(0),1:jpkm1) = zprdia(A2D(0),1:jpkm1) / ( zprmaxd(A2D(0),1:jpkm1) + rtrn ) * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "LDlight", zw3d )
+       ENDIF
+       !
+       IF( iom_use ( "MunetP" ) ) THEN ! Realized growth rate for picophyto
+       CALL iom_put( "MunetP"  , ( tr(:,:,:,jppic,Krhs)/rfact2/(tr(:,:,:,jppic,Kbb)+ rtrn ) * tmask(:,:,:)) ) 
        ENDIF
-       CALL iom_put( "Mumax"   , zprmaxn(:,:,:) * tmask(:,:,:)  ) ! Maximum growth rate
-       CALL iom_put( "MuP"     , zprpic(:,:,:) * xlimpic(:,:,:) * tmask(:,:,:) ) ! Realized growth rate for picophyto
-       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( "LPlight" , zprpic(:,:,:) / (zprmaxp(:,:,:) + rtrn) * tmask(:,:,:)  )  ! light limitation term
-       CALL iom_put( "LNlight" , zprbio(:,:,:) / (zprmaxn(:,:,:) + rtrn) * tmask(:,:,:)  )  ! light limitation term
-       CALL iom_put( "LDlight" , zprdia(:,:,:) / (zprmaxd(:,:,:) + rtrn) * tmask(:,:,:)   )
-       CALL iom_put( "MunetP"  , ( tr(:,:,:,jppic,Krhs)/rfact2/(tr(:,:,:,jppic,Kbb)+ rtrn ) * tmask(:,:,:)) ) ! Realized growth rate for picophyto
-       CALL iom_put( "MunetN"  , ( tr(:,:,:,jpphy,Krhs)/rfact2/(tr(:,:,:,jpphy,Kbb)+ rtrn ) * tmask(:,:,:)) ) ! Realized growth rate for picophyto
-       CALL iom_put( "MunetD"  , ( tr(:,:,:,jpdia,Krhs)/rfact2/(tr(:,:,:,jpdia,Kbb)+ rtrn ) * tmask(:,:,:)) ) ! Realized growth rate for picophyto
-       CALL iom_put( "TPP"     , ( zprorcap(:,:,:) + zprorcan(:,:,:) + zprorcad(:,:,:) ) * zfact * tmask(:,:,:)  )  ! total primary production
-       CALL iom_put( "TPNEW"   , ( zpronewp(:,:,:) + zpronewn(:,:,:) + zpronewd(:,:,:) ) * zfact * tmask(:,:,:)  ) ! total new production
-       CALL iom_put( "TPBFE"   , ( zprofep (:,:,:) + zprofen (:,:,:) + zprofed (:,:,:) ) * zfact * tmask(:,:,:)  )  ! total biogenic iron production
-       CALL iom_put( "tintpp"  , tpp * zfact )  !  global total integrated primary production molC/s
+       !
+       IF( iom_use ( "MunetN" ) ) THEN ! Realized growth rate for nanophyto
+       CALL iom_put( "MunetN"  , ( tr(:,:,:,jpphy,Krhs)/rfact2/(tr(:,:,:,jpphy,Kbb)+ rtrn ) * tmask(:,:,:)) ) 
+       ENDIF
+       !
+       IF( iom_use ( "MunetD" ) ) THEN ! Realized growth rate for diatoms
+       CALL iom_put( "MunetD"  , ( tr(:,:,:,jpdia,Krhs)/rfact2/(tr(:,:,:,jpdia,Kbb)+ rtrn ) * tmask(:,:,:)) ) 
+       ENDIF
+       !
+       IF( iom_use ( "TPP" ) ) THEN   ! total primary production
+         zw3d(A2D(0),1:jpkm1) = ( zprorcap(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) = ( zpronewp(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) = ( zprofep(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)
diff --git a/src/TOP/PISCES/SED/sedchem.F90 b/src/TOP/PISCES/SED/sedchem.F90
index afdc80a1c..8c5d9e50f 100644
--- a/src/TOP/PISCES/SED/sedchem.F90
+++ b/src/TOP/PISCES/SED/sedchem.F90
@@ -138,7 +138,7 @@ CONTAINS
       IF (ln_sediment_offline) THEN
          CALL sed_chem_cst
       ELSE
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
             ikt = mbkt(ji,jj) 
             IF ( tmask(ji,jj,ikt) == 1 ) THEN
                zchem_data(ji,jj,1) = ak13  (ji,jj,ikt)
diff --git a/src/TOP/PISCES/SED/sedsfc.F90 b/src/TOP/PISCES/SED/sedsfc.F90
index 460d760d9..9e36c0727 100644
--- a/src/TOP/PISCES/SED/sedsfc.F90
+++ b/src/TOP/PISCES/SED/sedsfc.F90
@@ -49,7 +49,7 @@ CONTAINS
       CALL unpack_arr ( jpoce, trc_data(1:jpi,1:jpj,8), iarroce(1:jpoce), pwcp(1:jpoce,1,jwfe2) )
       CALL unpack_arr ( jpoce, trc_data(1:jpi,1:jpj,9), iarroce(1:jpoce), pwcp(1:jpoce,1,jwlgw) )
 
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          ikt = mbkt(ji,jj)
          IF ( tmask(ji,jj,ikt) == 1 ) THEN
             tr(ji,jj,ikt,jptal,Kbb) = trc_data(ji,jj,1)
diff --git a/src/TOP/PISCES/SED/trcdmp_sed.F90 b/src/TOP/PISCES/SED/trcdmp_sed.F90
index 190206e54..b17ee20af 100644
--- a/src/TOP/PISCES/SED/trcdmp_sed.F90
+++ b/src/TOP/PISCES/SED/trcdmp_sed.F90
@@ -92,7 +92,7 @@ CONTAINS
                jl = n_trc_index(jn) 
                CALL trc_dta( kt, jl, ztrcdta )   ! read tracer data at nit000
                !
-               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+               DO_2D( 0, 0, 0, 0 )
                   ikt = mbkt(ji,jj)
                   tr(ji,jj,ikt,jn,Kbb) = ztrcdta(ji,jj,ikt) + ( tr(ji,jj,ikt,jn,Kbb) -  ztrcdta(ji,jj,ikt) )     &
                   &                  * exp( -restosed(ji,jj,ikt) * dtsed )
-- 
GitLab