From aa622c73df4aca74b8b693283770b8680fbd8154 Mon Sep 17 00:00:00 2001
From: cetlod <Christian.Ethe@ipsl.fr>
Date: Mon, 8 Aug 2022 14:58:33 +0200
Subject: [PATCH] Finalizing the work on halo in TOP modules

---
 .../EXPREF/file_def_nemo-innerttrc.xml        |  15 +-
 cfgs/SHARED/field_def_nemo-innerttrc.xml      |   8 +-
 cfgs/SHARED/field_def_nemo-pisces.xml         | 122 +++++++--------
 src/TOP/C14/sms_c14.F90                       |   8 +-
 src/TOP/C14/trcatm_c14.F90                    |  19 ++-
 src/TOP/C14/trcsms_c14.F90                    |  12 +-
 src/TOP/C14/trcwri_c14.F90                    |  78 +++++-----
 src/TOP/CFC/trcini_cfc.F90                    |   2 +-
 src/TOP/CFC/trcsms_cfc.F90                    |  10 +-
 src/TOP/PISCES/P2Z/p2zbio.F90                 |  28 ++--
 src/TOP/PISCES/P2Z/p2zexp.F90                 |  31 ++--
 src/TOP/PISCES/P2Z/p2zopt.F90                 |  22 ++-
 src/TOP/PISCES/P2Z/p2zsed.F90                 |  27 ++--
 src/TOP/PISCES/P4Z/p4zbc.F90                  |  89 ++++++++---
 src/TOP/PISCES/P4Z/p4zche.F90                 |  26 ++--
 src/TOP/PISCES/P4Z/p4zfechem.F90              | 147 +++++++++++-------
 src/TOP/PISCES/P4Z/p4zflx.F90                 |  44 ++++--
 src/TOP/PISCES/P4Z/p4zint.F90                 |  12 +-
 src/TOP/PISCES/P4Z/p4zlim.F90                 |  20 +--
 src/TOP/PISCES/P4Z/p4zlys.F90                 |   6 +-
 src/TOP/PISCES/P4Z/p4zmeso.F90                |   4 +-
 src/TOP/PISCES/P4Z/p4zopt.F90                 | 147 ++++++++++--------
 src/TOP/PISCES/P4Z/p4zpoc.F90                 |   2 +-
 src/TOP/PISCES/P4Z/p4zprod.F90                |   6 +-
 src/TOP/PISCES/P4Z/p4zsed.F90                 |  69 ++++----
 src/TOP/PISCES/P4Z/p4zsink.F90                |  53 ++++---
 src/TOP/PISCES/P4Z/p4zsms.F90                 | 129 ++++++++++-----
 src/TOP/PISCES/P4Z/p5zlim.F90                 |  32 ++--
 src/TOP/PISCES/P4Z/p5zmeso.F90                |   2 +-
 src/TOP/PISCES/P4Z/p5zprod.F90                |  11 +-
 src/TOP/PISCES/sms_pisces.F90                 |  46 +++---
 src/TOP/TRP/trcdmp.F90                        |   2 +-
 src/TOP/TRP/trcsink.F90                       |  32 ++--
 src/TOP/trc.F90                               |   8 +-
 src/TOP/trcini.F90                            |   5 +-
 src/TOP/trcopt.F90                            | 125 ++++++++-------
 36 files changed, 814 insertions(+), 585 deletions(-)

diff --git a/cfgs/ORCA2_OFF_TRC/EXPREF/file_def_nemo-innerttrc.xml b/cfgs/ORCA2_OFF_TRC/EXPREF/file_def_nemo-innerttrc.xml
index 747e6cd31..d1d25f865 100644
--- a/cfgs/ORCA2_OFF_TRC/EXPREF/file_def_nemo-innerttrc.xml
+++ b/cfgs/ORCA2_OFF_TRC/EXPREF/file_def_nemo-innerttrc.xml
@@ -36,12 +36,23 @@
            <field field_ref="ssh"    name="zos"  />
          </file>
 
-	 <file id="file1" name_suffix="_trc" description="passive tracers variables" >
+	 <file id="file2" name_suffix="_trc" description="passive tracers variables" >
            <field field_ref="Age"     name="Age"     operation="average" freq_op="1y"  > @Age_e3t / @e3t   </field>		      
            <field field_ref="CFC11"   name="CFC11"   operation="average" freq_op="1y"  > @CFC11_e3t / @e3t </field>
            <field field_ref="CFC12"   name="CFC12"   operation="average" freq_op="1y"  > @CFC12_e3t / @e3t </field>
            <field field_ref="SF6"     name="SF6"     operation="average" freq_op="1y"  > @SF6_e3t / @e3t   </field>
-           <field field_ref="RC14"    name="RC14"    operation="average" freq_op="1y"  > @RC14_e3t / @e3t  </field>
+	   <field field_ref="RC14"    name="RC14"    operation="average" freq_op="1y"  > @RC14_e3t / @e3t  </field>
+           <field field_ref="qtr_CFC11"    />
+           <field field_ref="qint_CFC11"   />
+           <field field_ref="qtr_CFC12"    />
+           <field field_ref="qint_CFC12"  />
+           <field field_ref="qtr_SF6"     />
+           <field field_ref="qint_SF6"    />
+           <field field_ref="qtr_c14"    />
+           <field field_ref="qint_c14"  />
+           <field field_ref="DeltaC14"    />
+           <field field_ref="C14Age"      />
+           <field field_ref="RAge"       />
          </file>
 
       </file_group>
diff --git a/cfgs/SHARED/field_def_nemo-innerttrc.xml b/cfgs/SHARED/field_def_nemo-innerttrc.xml
index d54a667b2..ccf4b8b58 100644
--- a/cfgs/SHARED/field_def_nemo-innerttrc.xml
+++ b/cfgs/SHARED/field_def_nemo-innerttrc.xml
@@ -16,7 +16,7 @@
   -->
 
 
-  <field_group id="inerttrc" grid_ref="grid_T_2D">
+  <field_group id="inerttrc" grid_ref="grid_T_2D_inner">
 
     <!-- CFC11 : variables available with ln_cfc11 -->
     <field id="CFC11"        long_name="Chlorofluoro carbon11 Concentration"      unit="umol/m3"   grid_ref="grid_T_3D" />
@@ -39,8 +39,8 @@
     <!-- C14 : variables available with ln_c14 -->
     <field id="RC14"         long_name="Radiocarbon ratio"                        unit="-"         grid_ref="grid_T_3D"  />
     <field id="RC14_e3t"     long_name="RC14 * e3t"                               unit="m"         grid_ref="grid_T_3D"   > RC14 * e3t </field >
-    <field id="DeltaC14"     long_name="Delta C14"                                unit="permil"    grid_ref="grid_T_3D"  />
-    <field id="C14Age"       long_name="Radiocarbon age"                          unit="yr"        grid_ref="grid_T_3D"  />
+    <field id="DeltaC14"     long_name="Delta C14"                                unit="permil"    grid_ref="grid_T_3D_inner"  />
+    <field id="C14Age"       long_name="Radiocarbon age"                          unit="yr"        grid_ref="grid_T_3D_inner"  />
     <field id="RAge"         long_name="Reservoir Age"                            unit="yr"       />
     <field id="qtr_c14"      long_name="Air-sea flux of C14"                      unit="1/m2/s"   />
     <field id="qint_c14"     long_name="Cumulative air-sea flux of C14"           unit="1/m2"     />
@@ -52,7 +52,7 @@
 
     <!-- AGE : variables available with ln_age -->
     <field id="Age"          long_name="Sea water age since surface contact"      unit="yr"        grid_ref="grid_T_3D"  />
-    <field id="Age_e3t"      long_name="Age * e3t"                                unit="yr * m"   grid_ref="grid_T_3D"    > Age * e3t </field >
+    <field id="Age_e3t"      long_name="Age * e3t"                                unit="yr * m"    grid_ref="grid_T_3D"    > Age * e3t </field >
 
   </field_group>
 
diff --git a/cfgs/SHARED/field_def_nemo-pisces.xml b/cfgs/SHARED/field_def_nemo-pisces.xml
index a94dec2ab..54c42a51b 100644
--- a/cfgs/SHARED/field_def_nemo-pisces.xml
+++ b/cfgs/SHARED/field_def_nemo-pisces.xml
@@ -172,12 +172,12 @@
 
   <!-- PISCES additional diagnostics on T grid  -->
 
-  <field_group id="diad_T" grid_ref="grid_T_2D">
+  <field_group id="diad_T" grid_ref="grid_T_2D_inner" >
     <field id="PH"          long_name="PH"                                      unit="1"          grid_ref="grid_T_3D_inner" />
     <field id="CO3"         long_name="Bicarbonates"                            unit="mol/m3"     grid_ref="grid_T_3D_inner" />
     <field id="CO3sat"      long_name="CO3 saturation"                          unit="mol/m3"     grid_ref="grid_T_3D_inner" />
     <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="PAR"         long_name="Photosynthetically Available Radiation"  unit="W/m2"       grid_ref="grid_T_3D_inner" />
     <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_inner" />
     <field id="PPPHYD"      long_name="Primary production of diatoms"           unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
@@ -193,17 +193,17 @@
     <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="DENIT"       long_name="Anoxic remineralization of OM"           unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
     <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="Nfix"        long_name="Nitrogen fixation"                       unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
     <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_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" />
-    <field id="MunetD"      long_name="Net growth rate for diatomes"            unit="s-1"        grid_ref="grid_T_3D" />
+    <field id="MunetN"      long_name="Net growth rate for nanophyto"           unit="s-1"        grid_ref="grid_T_3D_inner" />
+    <field id="MunetP"      long_name="Net growth rate for picophyto"           unit="s-1"        grid_ref="grid_T_3D_inner" />
+    <field id="MunetD"      long_name="Net growth rate for diatomes"            unit="s-1"        grid_ref="grid_T_3D_inner" />
     <field id="LNnut"       long_name="Nutrient limitation term in Nanophyto"   unit=""           grid_ref="grid_T_3D_inner" />
     <field id="LPnut"       long_name="Nutrient limitation term in Picophyto"   unit="-"          grid_ref="grid_T_3D_inner" />
     <field id="LDnut"       long_name="Nutrient limitation term in Diatoms"     unit=""           grid_ref="grid_T_3D_inner" />
@@ -219,60 +219,60 @@
     <field id="RASSD"       long_name="Size of the protein machinery (Diat.)"   unit="-"          grid_ref="grid_T_3D_inner" />
     <field id="RASSN"       long_name="Size of the protein machinery (Nano.)"   unit="-"          grid_ref="grid_T_3D_inner" />
     <field id="RASSP"       long_name="Size of the protein machinery (Pico.)"   unit="-"          grid_ref="grid_T_3D_inner" />
-    <field id="Fe3"         long_name="Iron III concentration"                  unit="nmol/m3"    grid_ref="grid_T_3D" />
-    <field id="FeL1"        long_name="Complexed Iron concentration with L1"    unit="nmol/m3"    grid_ref="grid_T_3D" />
-    <field id="TL1"         long_name="Total L1 concentration"                  unit="nmol/m3"    grid_ref="grid_T_3D" />
-    <field id="pdust"       long_name="dust concentration"                      unit="g/m3"                            />
-    <field id="Totlig"      long_name="Total ligand concentation"               unit="nmol/m3"    grid_ref="grid_T_3D" />
-    <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"  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"                        />
-    <field id="EPSI100"     long_name="Export of Silicate at 100 m"             unit="mol/m2/s"                        />
-    <field id="EPCAL100"    long_name="Export of Calcite at 100 m"              unit="mol/m2/s"                        />
-    <field id="EXPC"        long_name="Export of carbon"                        unit="mol/m2/s"   grid_ref="grid_T_3D" />
-    <field id="EXPFE"       long_name="Export of biogenic iron"                 unit="mol/m2/s"   grid_ref="grid_T_3D" />
-    <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"   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"                               />
-    <field id="Irondep"     long_name="Iron deposition from dust"               unit="mol/m2/s"                        />
-    <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_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_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_inner" />
-    <field id="LIGPR"       long_name="Photochemical loss of ligands"           unit="nmol-L/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="Fe3"         long_name="Iron III concentration"                  unit="nmol/m3"    grid_ref="grid_T_3D_inner" />
+    <field id="FeL1"        long_name="Complexed Iron concentration with L1"    unit="nmol/m3"    grid_ref="grid_T_3D_inner" />
+    <field id="TL1"         long_name="Total L1 concentration"                  unit="nmol/m3"    grid_ref="grid_T_3D_inner" />
+    <field id="pdust"       long_name="dust concentration"                      unit="g/m3"         />
+    <field id="Totlig"      long_name="Total ligand concentation"               unit="nmol/m3"    grid_ref="grid_T_3D_inner" />
+    <field id="Biron"       long_name="Bioavailable iron"                       unit="nmol/m3"    grid_ref="grid_T_3D_inner" />
+    <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="HYDR"        long_name="Iron input from hydrothemal vents"       unit="mol/m2/s"   grid_ref="grid_T_3D_inner" />
+    <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"    />
+    <field id="EPSI100"     long_name="Export of Silicate at 100 m"             unit="mol/m2/s"    />
+    <field id="EPCAL100"    long_name="Export of Calcite at 100 m"              unit="mol/m2/s"    />
+    <field id="EXPC"        long_name="Export of carbon"                        unit="mol/m2/s"   grid_ref="grid_T_3D_inner" />
+    <field id="EXPFE"       long_name="Export of biogenic iron"                 unit="mol/m2/s"   grid_ref="grid_T_3D_inner" />
+    <field id="EXPSI"       long_name="Export of Silicate"                      unit="mol/m2/s"   grid_ref="grid_T_3D_inner" />
+    <field id="EXPCAL"      long_name="Export of Calcite"                       unit="mol/m2/s"   grid_ref="grid_T_3D_inner" />
+    <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="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"          />
+    <field id="Irondep"     long_name="Iron deposition from dust"               unit="mol/m2/s"    />
+    <field id="Ironsed"     long_name="Iron deposition from sediment"           unit="mol/m2/s"            grid_ref="grid_T_3D_inner" />
+    <field id="FESCAV"      long_name="Scavenging of Iron"                      unit="mmol-Fe/m3/s"        grid_ref="grid_T_3D_inner" />
+    <field id="FECOLL"      long_name="Colloidal Pumping of FeL"                unit="mmol-FeL/m3/s"       grid_ref="grid_T_3D_inner" />
+    <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_inner" />
+    <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_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_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" />
+    <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_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" />
 
     <!-- PISCES tracers trends -->
-    <field id="INTdtAlk"    long_name="Vertically int. of change of alkalinity"             unit="mol/m2/s"                       />
-    <field id="INTdtDIC"    long_name="Vertically int. of change of dissic    "             unit="mol/m2/s"                       />
-    <field id="INTdtFer"    long_name="Vertically int. of change of iron      "             unit="mol/m2/s"                       />
-    <field id="INTdtDIN"    long_name="Vertically int. of change of nitrogen  "             unit="mol/m2/s"                       />
-    <field id="INTdtDIP"    long_name="Vertically int. of change of phophate  "             unit="mol/m2/s"                       />
-    <field id="INTdtSil"    long_name="Vertically int. of change of silicon   "             unit="mol/m2/s"                       />
+    <field id="INTdtAlk"    long_name="Vertically int. of change of alkalinity"             unit="mol/m2/s" />
+    <field id="INTdtDIC"    long_name="Vertically int. of change of dissic    "             unit="mol/m2/s" />
+    <field id="INTdtFer"    long_name="Vertically int. of change of iron      "             unit="mol/m2/s" />
+    <field id="INTdtDIN"    long_name="Vertically int. of change of nitrogen  "             unit="mol/m2/s" />
+    <field id="INTdtDIP"    long_name="Vertically int. of change of phophate  "             unit="mol/m2/s" />
+    <field id="INTdtSil"    long_name="Vertically int. of change of silicon   "             unit="mol/m2/s" />
 
 
     <!-- dbio_T on T grid : variables available with diaar5 -->
@@ -292,9 +292,9 @@
     <field id="INTPBSI"     long_name="Vertically integrated of biogenic Si production"     unit="mol/m2/s"  grid_ref="grid_T_vsum"  detect_missing_value="true" > PBSi * e3t </field >
 
     <!-- PISCES light : variables available with key_pisces_reduced -->
-    <field id="FNO3PHY"     long_name="FNO3PHY"                                 unit=""          grid_ref="grid_T_3D" />
-    <field id="FNH4PHY"     long_name="FNH4PHY"                                 unit=""          grid_ref="grid_T_3D" />
-    <field id="FNH4NO3"     long_name="FNH4NO3"                                 unit=""          grid_ref="grid_T_3D" />
+    <field id="FNO3PHY"     long_name="FNO3PHY"                                 unit=""          grid_ref="grid_T_3D_inner" />
+    <field id="FNH4PHY"     long_name="FNH4PHY"                                 unit=""          grid_ref="grid_T_3D_inner" />
+    <field id="FNH4NO3"     long_name="FNH4NO3"                                 unit=""          grid_ref="grid_T_3D_inner" />
     <field id="TNO3PHY"     long_name="TNO3PHY"                                 unit=""  />
     <field id="TNH4PHY"     long_name="TNH4PHY"                                 unit=""  />
     <field id="TPHYDOM"     long_name="TPHYDOM"                                 unit=""  />
diff --git a/src/TOP/C14/sms_c14.F90 b/src/TOP/C14/sms_c14.F90
index bdea77765..a9971f6a2 100644
--- a/src/TOP/C14/sms_c14.F90
+++ b/src/TOP/C14/sms_c14.F90
@@ -51,6 +51,8 @@ MODULE sms_c14
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   spco2      ! Atmospheric CO2
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   tyrco2     ! Time (yr) atmospheric CO2 data
 
+         !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
    !! $Id: sms_c14.F90 10071 2018-08-28 14:49:04Z nicolasmartin $ 
@@ -64,9 +66,9 @@ CONTAINS
       !!                  ***  ROUTINE trc_sms_c14_alloc  ***
       !!----------------------------------------------------------------------
       sms_c14_alloc = 0
-      ALLOCATE( exch_c14(jpi,jpj)        ,  exch_co2(jpi,jpj)        ,   &
-         &      qtr_c14(jpi,jpj)         ,  qint_c14(jpi,jpj)        ,   &
-         &      c14sbc(jpi,jpj)          ,  STAT = sms_c14_alloc )
+      ALLOCATE( exch_c14(A2D(0))        ,  exch_co2(A2D(0))        ,   &
+         &      qtr_c14(A2D(0))         ,  qint_c14(A2D(0))        ,   &
+         &      c14sbc(A2D(0))          ,  STAT = sms_c14_alloc )
          !
       !
    END FUNCTION sms_c14_alloc
diff --git a/src/TOP/C14/trcatm_c14.F90 b/src/TOP/C14/trcatm_c14.F90
index 54a768067..d09dcf5df 100644
--- a/src/TOP/C14/trcatm_c14.F90
+++ b/src/TOP/C14/trcatm_c14.F90
@@ -59,6 +59,13 @@ CONTAINS
       !
       tyrc14_now = 0._wp   ! initialize
       !
+      IF( kc14typ == 0) THEN
+         co2sbc=pco2at
+         DO_2D( 0, 0, 0, 0 )
+            c14sbc(ji,jj) = rc14at
+         END_2D
+      ENDIF
+      !
       IF(kc14typ >= 1) THEN  ! Transient atmospheric forcing: CO2
       !
          clfile = TRIM( cfileco2 )
@@ -116,10 +123,10 @@ CONTAINS
        ! Linear  interpolation of the C-14 source fonction
        ! in linear latitude bands  (20N,40N) and (20S,40S)
        !------------------------------------------------------
-            ALLOCATE( fareaz  (jpi,jpj ,nc14zon) , STAT=ierr3 )
+            ALLOCATE( fareaz(A2D(0) ,nc14zon) , STAT=ierr3 )
             IF( ierr3 /= 0 )   CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate fareaz' )
       !
-            DO_2D( 1, 1, 1, 1 )                 ! from C14b package
+            DO_2D( 0, 0, 0, 0 )                 ! from C14b package
               IF( gphit(ji,jj) >= yn40 ) THEN
                  fareaz(ji,jj,1) = 0.
                  fareaz(ji,jj,2) = 0.
@@ -205,9 +212,9 @@ CONTAINS
       !! ** Action  :   atmospheric values interpolated at time-step kt
       !!----------------------------------------------------------------------
       INTEGER                 , INTENT(in   )   ::   kt       ! ocean time-step
-      REAL(wp), DIMENSION(:,:), INTENT(  out)   ::   c14sbc   ! atm c14 ratio
+      REAL(wp), DIMENSION(A2D(0)), INTENT(  out)   ::   c14sbc   ! atm c14 ratio
       REAL(wp),                 INTENT(  out)   ::   co2sbc   ! atm co2 p
-      INTEGER                                   ::   jz       ! dummy loop indice
+      INTEGER                                   ::   ji, jj, jz       ! dummy loop indice
       REAL(wp)                              ::   zdint,zint   ! work
       REAL(wp), DIMENSION(nc14zon)              ::   zonbc14  ! work
       !
@@ -215,10 +222,6 @@ CONTAINS
       !
       IF( ln_timing )   CALL timing_start('trc_atm_c14')
       !
-      IF( kc14typ == 0) THEN
-         co2sbc=pco2at
-         c14sbc(:,:)=rc14at
-      ENDIF
       !
       IF(kc14typ >= 1) THEN  ! Transient C14 & CO2
       !
diff --git a/src/TOP/C14/trcsms_c14.F90 b/src/TOP/C14/trcsms_c14.F90
index 9ce0a584a..0135dc8c9 100644
--- a/src/TOP/C14/trcsms_c14.F90
+++ b/src/TOP/C14/trcsms_c14.F90
@@ -80,7 +80,7 @@ CONTAINS
       !  CO2 solubility (Weiss, 1974; Wanninkhof, 2014) 
       ! -------------------------------------------------------------------
 
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          IF( tmask(ji,jj,1) >  0. ) THEN
             !
             zt   = MIN( 40. , ts(ji,jj,1,jp_tem,Kmm) )
@@ -121,21 +121,21 @@ CONTAINS
       !
       ! Flux of C-14 from air-to-sea; units: (C14/C ratio) x m/s
       !                               already masked
-      qtr_c14(:,:) = exch_c14(:,:) * ( c14sbc(:,:) - tr(:,:,1,jp_c14,Kbb) )
+      DO_2D( 0, 0, 0, 0 )
+         qtr_c14(ji,jj) = exch_c14(ji,jj) * ( c14sbc(ji,jj) - tr(ji,jj,1,jp_c14,Kbb) )
+      END_2D
             
       ! cumulation of air-to-sea flux at each time step
       qint_c14(:,:) = qint_c14(:,:) + qtr_c14(:,:) * rn_Dt
       !
       ! Add the surface flux to the trend of jp_c14
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          tr(ji,jj,1,jp_c14,Krhs) = tr(ji,jj,1,jp_c14,Krhs) + qtr_c14(ji,jj) / e3t(ji,jj,1,Kmm) 
       END_2D
       !
       ! Computation of decay effects on jp_c14
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
-         !
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          tr(ji,jj,jk,jp_c14,Krhs) = tr(ji,jj,jk,jp_c14,Krhs) - rlam14 * tr(ji,jj,jk,jp_c14,Kbb) * tmask(ji,jj,jk) 
-         !
       END_3D
       !
       IF( lrst_trc ) THEN
diff --git a/src/TOP/C14/trcwri_c14.F90 b/src/TOP/C14/trcwri_c14.F90
index 7f3de9f0c..7fcda51a3 100644
--- a/src/TOP/C14/trcwri_c14.F90
+++ b/src/TOP/C14/trcwri_c14.F90
@@ -38,8 +38,8 @@ CONTAINS
       CHARACTER (len=20)   :: cltra         ! short title for tracer
       INTEGER              :: ji,jj,jk,jn   ! dummy loop indexes
       REAL(wp)             :: zage,zarea,ztemp   ! temporary
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zres, z2d ! temporary storage 2D
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z3d , zz3d ! temporary storage 3D
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: z2d ! temporary storage 2D
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z3d ! temporary storage 3D
       !!---------------------------------------------------------------------
  
       ! write the tracer concentrations in the file
@@ -49,41 +49,35 @@ CONTAINS
 
       ! compute and write the tracer diagnostic in the file
       ! ---------------------------------------
+      IF( iom_use("qtr_c14") ) CALL iom_put( "qtr_c14" , rsiyea * qtr_c14(:,:)  )   !  Radiocarbon surf flux [./m2/yr]
+                               CALL iom_put( "qint_c14", qint_c14(:,:)  )         ! cumulative flux [./m2]
       
       IF( iom_use("DeltaC14") .OR. iom_use("C14Age") .OR. iom_use("RAge")   ) THEN
          !
-         ALLOCATE( z2d(jpi,jpj), zres(jpi,jpj) )
-         ALLOCATE( z3d(jpi,jpj,jpk), zz3d(jpi,jpj,jpk) )
+         ALLOCATE( z2d(A2D(0)), z3d(A2D(0),jpk) )
          !
          zage = -1._wp / rlam14 / rsiyea  ! factor for radioages in year
          z3d(:,:,:)  = 1._wp
-         zz3d(:,:,:) = 0._wp
          !
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             IF( tmask(ji,jj,jk) > 0._wp) THEN
-               z3d (ji,jj,jk) = tr(ji,jj,jk,jp_c14,Kmm)
-               zz3d(ji,jj,jk) = LOG( z3d(ji,jj,jk) )
+               z3d(ji,jj,jk) = tr(ji,jj,jk,jp_c14,Kmm)
             ENDIF
          END_3D
-         zres(:,:) = z3d(:,:,1)
+         CALL iom_put( "C14Age", zage * LOG( z3d(:,:,:) ) )            !  Radiocarbon age [yr]
 
          ! Reservoir age [yr]
-         z2d(:,:) =0._wp
-         jk = 1
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-            ztemp = zres(ji,jj) / c14sbc(ji,jj)
-            IF( ztemp > 0._wp .AND. tmask(ji,jj,jk) > 0._wp ) z2d(ji,jj) = LOG( ztemp )
+         z2d(:,:) = 0._wp
+         DO_2D( 0, 0, 0, 0 )
+            ztemp = z3d(ji,jj,1) / c14sbc(ji,jj)
+            IF( ztemp > 0._wp .AND. tmask(ji,jj,1) > 0._wp ) z2d(ji,jj) = LOG( ztemp )
          END_2D
+         CALL iom_put( "RAge" , zage * z2d(:,:) )                     ! Reservoir age [yr]
          !
          z3d(:,:,:) = 1.d03 * ( z3d(:,:,:) - 1._wp )
          CALL iom_put( "DeltaC14" , z3d(:,:,:)  )  ! Delta C14 [permil]
-         CALL iom_put( "C14Age"   , zage * zz3d(:,:,:) )            !  Radiocarbon age [yr]
-
-         CALL iom_put( "qtr_c14", rsiyea * qtr_c14(:,:)  )            !  Radiocarbon surf flux [./m2/yr]
-         CALL iom_put( "qint_c14" , qint_c14  )                       ! cumulative flux [./m2]
-         CALL iom_put( "RAge" , zage * z2d(:,:) )                     ! Reservoir age [yr]
          !
-         DEALLOCATE( z2d, zres, z3d, zz3d )
+         DEALLOCATE( z2d, z3d )
          !
       ENDIF
       !
@@ -91,23 +85,35 @@ CONTAINS
       !
       CALL iom_put( "AtmCO2", co2sbc )  !     global atmospheric CO2 [ppm]
     
-      IF( iom_use("AtmC14") ) THEN
-         zarea = glob_sum( 'trcwri_c14', e1e2t(:,:) )           ! global ocean surface
-         ztemp = glob_sum( 'trcwri_c14', c14sbc(:,:) * e1e2t(:,:) )
-         ztemp = ( ztemp / zarea - 1._wp ) * 1000._wp
-         CALL iom_put( "AtmC14" , ztemp )   ! Global atmospheric DeltaC14 [permil]
-      ENDIF
-      IF( iom_use("K_C14") ) THEN
-         ztemp = glob_sum ( 'trcwri_c14', exch_c14(:,:) * e1e2t(:,:) )
-         ztemp = rsiyea * ztemp / zarea
-         CALL iom_put( "K_C14" , ztemp )   ! global mean exchange velocity for C14/C ratio [m/yr]
-      ENDIF
-      IF( iom_use("K_CO2") ) THEN
+      IF( iom_use("AtmC14") .OR. iom_use("K_C14") .OR. iom_use("K_CO2") ) THEN
          zarea = glob_sum( 'trcwri_c14', e1e2t(:,:) )           ! global ocean surface
-         ztemp = glob_sum ( 'trcwri_c14', exch_co2(:,:) * e1e2t(:,:) )
-         ztemp = 360000._wp * ztemp / zarea       ! cm/h units: directly comparable with literature
-         CALL iom_put( "K_CO2", ztemp )  !  global mean CO2 piston velocity [cm/hr]
-      ENDIF
+         ALLOCATE( z2d(A2D(0)) )
+         IF( iom_use("AtmC14") ) THEN
+            DO_2D( 0, 0, 0, 0 )
+               z2d(ji,jj) = c14sbc(ji,jj) * e1e2t(ji,jj)
+            END_2D
+            ztemp = glob_sum( 'trcwri_c14', z2d(:,:) )
+            ztemp = ( ztemp / zarea - 1._wp ) * 1000._wp
+            CALL iom_put( "AtmC14" , ztemp )   ! Global atmospheric DeltaC14 [permil]
+         ENDIF
+         IF( iom_use("K_C14") ) THEN
+             DO_2D( 0, 0, 0, 0 )
+                z2d(ji,jj) = exch_c14(ji,jj) * e1e2t(ji,jj)
+             END_2D
+             ztemp = glob_sum( 'trcwri_c14', z2d(:,:) )
+             ztemp = rsiyea * ztemp / zarea
+             CALL iom_put( "K_C14" , ztemp )   ! global mean exchange velocity for C14/C ratio [m/yr]
+         ENDIF
+         IF( iom_use("K_CO2") ) THEN
+            DO_2D( 0, 0, 0, 0 )
+               z2d(ji,jj) = exch_co2(ji,jj) * e1e2t(ji,jj)
+            END_2D
+            ztemp = glob_sum( 'trcwri_c14', z2d(:,:) )
+            ztemp = 360000._wp * ztemp / zarea       ! cm/h units: directly comparable with literature
+            CALL iom_put( "K_CO2", ztemp )  !  global mean CO2 piston velocity [cm/hr]
+         ENDIF
+         DEALLOCATE( z2d )
+      END IF
       IF( iom_use("C14Inv") ) THEN
          ztemp = glob_sum( 'trcwri_c14', tr(:,:,:,jp_c14,Kmm) * cvol(:,:,:) )
          ztemp = atomc14 * xdicsur * ztemp
diff --git a/src/TOP/CFC/trcini_cfc.F90 b/src/TOP/CFC/trcini_cfc.F90
index cacfee081..5d759a136 100644
--- a/src/TOP/CFC/trcini_cfc.F90
+++ b/src/TOP/CFC/trcini_cfc.F90
@@ -131,7 +131,7 @@ CONTAINS
       ! Linear interpolation between 2 hemispheric function of latitud between ylats and ylatn
       !---------------------------------------------------------------------------------------
       zyd = ylatn - ylats      
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          IF(     gphit(ji,jj) >= ylatn ) THEN   ;   xphem(ji,jj) = 1.e0
          ELSEIF( gphit(ji,jj) <= ylats ) THEN   ;   xphem(ji,jj) = 0.e0
          ELSE                                   ;   xphem(ji,jj) = ( gphit(ji,jj) - ylats) / zyd
diff --git a/src/TOP/CFC/trcsms_cfc.F90 b/src/TOP/CFC/trcsms_cfc.F90
index 32e9b63eb..e473c088c 100644
--- a/src/TOP/CFC/trcsms_cfc.F90
+++ b/src/TOP/CFC/trcsms_cfc.F90
@@ -124,9 +124,9 @@ CONTAINS
                &           +  atm_cfc(iyear_end, jm, jl) * REAL(im2, wp) ) / 12.
          END DO
          
-         !                                                         !------------!
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )                   !  i-j loop  !
-            !                                                      !------------!
+         !                                     !------------!
+         DO_2D( 0, 0, 0, 0 )                   !  i-j loop  !
+            !                                  !------------!
             ! space interpolation
             zpp_cfc  =       xphem(ji,jj)   * zpatm(1,jl)   &
                &     + ( 1.- xphem(ji,jj) ) * zpatm(2,jl)
@@ -309,8 +309,8 @@ CONTAINS
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE trc_sms_cfc_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( xphem   (jpi,jpj)        , atm_cfc(jpyear,jphem,jp_cfc)  ,    &
-         &      qtr_cfc (jpi,jpj,jp_cfc) , qint_cfc(jpi,jpj,jp_cfc)      ,    &
+      ALLOCATE( xphem   (A2D(0))        , atm_cfc(jpyear,jphem,jp_cfc)  ,    &
+         &      qtr_cfc (A2D(0),jp_cfc) , qint_cfc(A2D(0),jp_cfc)      ,    &
          &      soa(4,jp_cfc)    ,  sob(3,jp_cfc)   ,  sca(5,jp_cfc)     ,    &
          &      STAT=trc_sms_cfc_alloc )
          !
diff --git a/src/TOP/PISCES/P2Z/p2zbio.F90 b/src/TOP/PISCES/P2Z/p2zbio.F90
index ecbfc9f47..babf325b7 100644
--- a/src/TOP/PISCES/P2Z/p2zbio.F90
+++ b/src/TOP/PISCES/P2Z/p2zbio.F90
@@ -104,7 +104,6 @@ CONTAINS
       !
       IF( ln_timing )   CALL timing_start('p2z_bio')
       !
-      IF( lk_iomput )   ALLOCATE( zw2d(jpi,jpj,17), zw3d(jpi,jpj,jpk,3) )
 
       IF( kt == nittrc000 ) THEN
          IF(lwp) WRITE(numout,*)
@@ -112,18 +111,18 @@ CONTAINS
          IF(lwp) WRITE(numout,*) ' ~~~~~~~'
       ENDIF
 
-      xksi(:,:) = 0.e0        ! zooplakton closure ( fbod)
       IF( lk_iomput ) THEN
-         zw2d  (:,:,:) = 0._wp
-         zw3d(:,:,:,:) = 0._wp
+         ALLOCATE( zw3d(A2D(0),jpk,3) )   ;   zw3d(:,:,jpk,:) = 0._wp
+         ALLOCATE( zw2d(A2D(0),17) )
       ENDIF
+      !
+      xksi(:,:) = 0.e0        ! zooplakton closure ( fbod)
 
       !                                      ! -------------------------- !
-      DO jk = 1, jpkbm1                      !  Upper ocean (bio-layers)  !
+      DO_3D( 0, 0, 0, 0, 1, jpkbm1 )         !  Upper ocean (bio-layers)  !
          !                                   ! -------------------------- !
-         DO_2D( 0, 0, 0, 0 )
-            ! trophic variables( det, zoo, phy, no3, nh4, dom)
-            ! ------------------------------------------------
+         ! trophic variables( det, zoo, phy, no3, nh4, dom)
+         ! ------------------------------------------------
 
             ! negative trophic variables DO not contribute to the fluxes
             zdet = MAX( 0.e0, tr(ji,jj,jk,jpdet,Kmm) )
@@ -235,13 +234,11 @@ CONTAINS
                zw3d(ji,jj,jk,3) = znh4no3 * 86400   
                 ! 
              ENDIF
-         END_2D
-      END DO
+      END_3D
 
       !                                      ! -------------------------- !
-      DO jk = jpkb, jpkm1                    !  Upper ocean (bio-layers)  !
+      DO_3D( 0, 0, 0, 0, jpkb, jpkm1 )       !  Upper ocean (bio-layers)  !
          !                                   ! -------------------------- !
-         DO_2D( 0, 0, 0, 0 )
             ! remineralisation of all quantities towards nitrate 
 
             !    trophic variables( det, zoo, phy, no3, nh4, dom)
@@ -334,12 +331,9 @@ CONTAINS
                zw3d(ji,jj,jk,3) = znh4no3 * 86400._wp
                !
             ENDIF
-         END_2D
-      END DO
+      END_3D
       !
       IF( lk_iomput ) THEN
-         CALL lbc_lnk( 'p2zbio', zw2d(:,:,:),'T', 1.0_wp )
-         CALL lbc_lnk( 'p2zbio', zw3d(:,:,:,1),'T', 1.0_wp, zw3d(:,:,:,2),'T', 1.0_wp, zw3d(:,:,:,3),'T', 1.0_wp )
          ! Save diagnostics
          CALL iom_put( "TNO3PHY", zw2d(:,:,1) )
          CALL iom_put( "TNH4PHY", zw2d(:,:,2) )
@@ -362,6 +356,8 @@ CONTAINS
          CALL iom_put( "FNH4PHY", zw3d(:,:,:,2) )
          CALL iom_put( "FNH4NO3", zw3d(:,:,:,3) )
          !
+         DEALLOCATE( zw2d, zw3d )
+         !
       ENDIF
 
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
diff --git a/src/TOP/PISCES/P2Z/p2zexp.F90 b/src/TOP/PISCES/P2Z/p2zexp.F90
index d5fe6fe11..71ed7709e 100644
--- a/src/TOP/PISCES/P2Z/p2zexp.F90
+++ b/src/TOP/PISCES/P2Z/p2zexp.F90
@@ -65,7 +65,7 @@ CONTAINS
       !!
       INTEGER  ::   ji, jj, jk, jl, ikt
       REAL(wp) ::   zgeolpoc, zfact, zwork, ze3t, zsedpocd, zmaskt
-      REAL(wp), DIMENSION(jpi,jpj)   ::  zsedpoca
+      REAL(wp), DIMENSION(A2D(0))   ::  zsedpoca
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
       !
@@ -106,7 +106,7 @@ CONTAINS
          tr(ji,jj,1,jpno3,Krhs) = tr(ji,jj,1,jpno3,Krhs) + zgeolpoc * cmask(ji,jj) / areacot / e3t(ji,jj,1,Kmm)
       END_2D
 
-      CALL lbc_lnk( 'p2zexp', sedpocn, 'T', 1.0_wp )
+!      CALL lbc_lnk( 'p2zexp', sedpocn, 'T', 1.0_wp )
  
       ! Oa & Ek: diagnostics depending on jpdia2d !          left as example
       IF( lk_iomput )  CALL iom_put( "SEDPOC" , sedpocn )
@@ -120,7 +120,7 @@ CONTAINS
         !                                              
       ELSE
         !
-        DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+        DO_2D( 0, 0, 0, 0 )
            zsedpocd = zsedpoca(ji,jj) - 2. * sedpocn(ji,jj) + sedpocb(ji,jj)      ! time laplacian on tracers
            sedpocb(ji,jj) = sedpocn(ji,jj) + rn_atfp * zsedpocd                     ! sedpocb <-- filtered sedpocn
            sedpocn(ji,jj) = zsedpoca(ji,jj)                                       ! sedpocn <-- sedpoca
@@ -156,8 +156,8 @@ CONTAINS
       INTEGER, INTENT(in)  ::  Kmm      ! time level index
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zmaskt, zfluo, zfluu
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zrro
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdm0
+      REAL(wp), DIMENSION(A2D(0)    ) :: zrro, zarea
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zdm0
       !!---------------------------------------------------------------------
       !
       IF(lwp) THEN
@@ -171,9 +171,9 @@ CONTAINS
       ! Calculate vertical distribution of newly formed biogenic poc
       ! in the water column in the case of max. possible bottom depth
       ! ------------------------------------------------------------
-      zdm0 = 0._wp
-      zrro = 1._wp
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, jpkb, jpkm1 )
+      zdm0(:,:,:) = 0._wp
+      zrro(:,:) = 1._wp
+      DO_3D( 0, 0, 0, 0, jpkb, jpkm1 )
          zfluo = ( gdepw(ji,jj,jk  ,Kmm) / gdepw(ji,jj,jpkb,Kmm) )**xhr
          zfluu = ( gdepw(ji,jj,jk+1,Kmm) / gdepw(ji,jj,jpkb,Kmm) )**xhr
          IF( zfluo.GT.1. )   zfluo = 1._wp
@@ -190,14 +190,14 @@ CONTAINS
       ! ----------------------------------------------------------------------
       dminl(:,:)   = 0._wp
       dmin3(:,:,:) = zdm0
-      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) == 0._wp ) THEN
             dminl(ji,jj) = dminl(ji,jj) + dmin3(ji,jj,jk)
             dmin3(ji,jj,jk) = 0._wp
          ENDIF
       END_3D
 
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          IF( tmask(ji,jj,1) == 0 )   dmin3(ji,jj,1) = 0._wp
       END_2D
 
@@ -209,8 +209,11 @@ CONTAINS
             IF( zmaskt == 0. )   cmask(ji,jj) = 1._wp
          END IF
       END_2D
-      CALL lbc_lnk( 'p2zexp', cmask , 'T', 1.0_wp )      ! lateral boundary conditions on cmask   (sign unchanged)
-      areacot = glob_sum( 'p2zexp', e1e2t(:,:) * cmask(:,:) )
+!      CALL lbc_lnk( 'p2zexp', cmask , 'T', 1.0_wp )      ! lateral boundary conditions on cmask   (sign unchanged)
+      DO_2D( 0, 0, 0, 0 )
+         zarea(ji,jj) = e1e2t(ji,jj) * cmask(ji,jj)  
+      END_2D
+      areacot = glob_sum( 'p2zexp', zarea(:,:) )
       !
       IF( ln_rsttr ) THEN
          CALL iom_get( numrtr, jpdom_auto, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) )
@@ -226,8 +229,8 @@ CONTAINS
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE p2z_exp_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( cmask(jpi,jpj) , dminl(jpi,jpj) , dmin3(jpi,jpj,jpk), &
-         &      sedpocb(jpi,jpj) , sedpocn(jpi,jpj),   STAT=p2z_exp_alloc )
+      ALLOCATE( cmask(A2D(0)) , dminl(A2D(0)) , dmin3(A2D(0),jpk), &
+         &      sedpocb(A2D(0)) , sedpocn(A2D(0)),   STAT=p2z_exp_alloc )
       IF( p2z_exp_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p2z_exp_alloc : failed to allocate arrays.' )
       !
    END FUNCTION p2z_exp_alloc
diff --git a/src/TOP/PISCES/P2Z/p2zopt.F90 b/src/TOP/PISCES/P2Z/p2zopt.F90
index 85f8b3e2a..4c8268dea 100644
--- a/src/TOP/PISCES/P2Z/p2zopt.F90
+++ b/src/TOP/PISCES/P2Z/p2zopt.F90
@@ -70,8 +70,8 @@ CONTAINS
       REAL(wp) ::   zpig                ! log of the total pigment
       REAL(wp) ::   zkr, zkg            ! total absorption coefficient in red and green
       REAL(wp) ::   zcoef               ! temporary scalar
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zpar100, zpar0m
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zparr, zparg
+      REAL(wp), DIMENSION(A2D(0)    ) :: zpar100, zpar0m
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zparr, zparg
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p2z_opt')
@@ -85,8 +85,14 @@ CONTAINS
 
       !                                          ! surface irradiance
       !                                          ! ------------------
-      IF( ln_dm2dc ) THEN   ;   zpar0m(:,:) = qsr_mean(:,:) * 0.43
-      ELSE                  ;   zpar0m(:,:) = qsr     (:,:) * 0.43
+      IF( ln_dm2dc ) THEN  
+         DO_2D( 0, 0, 0, 0 )
+            zpar0m(ji,jj) = qsr_mean(ji,jj) * 0.43
+         END_2D
+      ELSE  
+         DO_2D( 0, 0, 0, 0 )
+            zpar0m(ji,jj) = qsr(ji,jj) * 0.43
+         END_2D
       ENDIF
       zpar100(:,:)   = zpar0m(:,:) * 0.01
       zparr  (:,:,1) = zpar0m(:,:) * 0.5
@@ -94,14 +100,14 @@ CONTAINS
 
       !                                          ! Photosynthetically Available Radiation (PAR)
       zcoef = 12 * redf / rcchl / rpig           ! --------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpk )
+      DO_3D( 0, 0, 0, 0, 2, jpk )
          zpig = LOG(  MAX( TINY(0.), tr(ji,jj,jk-1,jpphy,Kmm) ) * zcoef  )
          zkr  = xkr0 + xkrp * EXP( xlr * zpig )
          zkg  = xkg0 + xkgp * EXP( xlg * zpig )
          zparr(ji,jj,jk) = zparr(ji,jj,jk-1) * EXP( -zkr * e3t(ji,jj,jk-1,Kmm) )
          zparg(ji,jj,jk) = zparg(ji,jj,jk-1) * EXP( -zkg * e3t(ji,jj,jk-1,Kmm) )
       END_3D
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) ! mean par at t-levels
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! mean par at t-levels
          zpig = LOG(  MAX( TINY(0.), tr(ji,jj,jk,jpphy,Kmm) ) * zcoef  )
          zkr  = xkr0 + xkrp * EXP( xlr * zpig )
          zkg  = xkg0 + xkgp * EXP( xlg * zpig )
@@ -113,11 +119,11 @@ CONTAINS
       !                                          ! Euphotic layer
       !                                          ! --------------
       neln(:,:) = 1                                   ! euphotic layer level
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )  ! (i.e. 1rst T-level strictly below EL bottom)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )  ! (i.e. 1rst T-level strictly below EL bottom)
         IF( etot(ji,jj,jk) >= zpar100(ji,jj) )   neln(ji,jj) = jk + 1 
       END_3D
       !                                               ! Euphotic layer depth
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
+      DO_2D( 0, 0, 0, 0 ) 
          heup(ji,jj) = gdepw(ji,jj,neln(ji,jj),Kmm)
       END_2D
 
diff --git a/src/TOP/PISCES/P2Z/p2zsed.F90 b/src/TOP/PISCES/P2Z/p2zsed.F90
index 66f24308c..50d7e7696 100644
--- a/src/TOP/PISCES/P2Z/p2zsed.F90
+++ b/src/TOP/PISCES/P2Z/p2zsed.F90
@@ -61,10 +61,11 @@ CONTAINS
       INTEGER, INTENT( in ) ::   kt         ! ocean time-step index      
       INTEGER, INTENT( in ) ::   Kmm, Krhs  ! time level indices
       !
-      INTEGER  ::   ji, jj, jk, jl, ierr
+      INTEGER  ::   ji, jj, jk
       CHARACTER (len=25) :: charout
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zw2d
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwork, ztra
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zwork
+      REAL(wp) :: ztra
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p2z_sed')
@@ -83,22 +84,26 @@ CONTAINS
       zwork(:,:,jpk) = 0.e0      ! bottom value  set to zero
 
       ! tracer flux at w-point: we use -vsed (downward flux)  with simplification : no e1*e2
-      DO jk = 2, jpkm1
-         zwork(:,:,jk) = -vsed * tr(:,:,jk-1,jpdet,Kmm)
-      END DO
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
+         zwork(ji,jj,jk) = -vsed * tr(ji,jj,jk-1,jpdet,Kmm)
+      END_3D
 
       ! tracer flux divergence at t-point added to the general trend
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
-         ztra(ji,jj,jk)  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
-         tr(ji,jj,jk,jpdet,Krhs) = tr(ji,jj,jk,jpdet,Krhs) + ztra(ji,jj,jk) 
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
+         ztra  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
+         tr(ji,jj,jk,jpdet,Krhs) = tr(ji,jj,jk,jpdet,Krhs) + ztra
       END_3D
 
       IF( lk_iomput )  THEN
          IF( iom_use( "TDETSED" ) ) THEN
-            ALLOCATE( zw2d(jpi,jpj) )
-            zw2d(:,:) =  ztra(:,:,1) * e3t(:,:,1,Kmm) * 86400._wp
+            ALLOCATE( zw2d(A2D(0)) )
+            DO_2D( 0, 0, 0, 0 ) 
+               zw2d(ji,jj)  = - ( zwork(ji,jj,1) - zwork(ji,jj,2) ) * 86400._wp
+            END_2D
             DO jk = 2, jpkm1
-               zw2d(:,:) = zw2d(:,:) + ztra(:,:,jk) * e3t(:,:,jk,Kmm) * 86400._wp
+               DO_2D( 0, 0, 0, 0 ) 
+                  zw2d(ji,jj) = zw2d(ji,jj) - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) * 86400._wp
+               END_2D
             END DO
             CALL iom_put( "TDETSED", zw2d )
             DEALLOCATE( zw2d )
diff --git a/src/TOP/PISCES/P4Z/p4zbc.F90 b/src/TOP/PISCES/P4Z/p4zbc.F90
index f0bd7da61..e4b20adbd 100644
--- a/src/TOP/PISCES/P4Z/p4zbc.F90
+++ b/src/TOP/PISCES/P4Z/p4zbc.F90
@@ -13,6 +13,7 @@ MODULE p4zbc
    USE sms_pisces      !  PISCES Source Minus Sink variables
    USE iom             !  I/O manager
    USE fldread         !  time interpolation
+   USE prtctl          !  print control for debugging
    USE trcbc
 
    IMPLICIT NONE
@@ -36,7 +37,6 @@ MODULE p4zbc
    LOGICAL , PUBLIC ::   ll_river     !: boolean for river input of nutrients
    LOGICAL , PUBLIC ::   ll_ndepo     !: boolean for atmospheric deposition of N
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_dust      ! structure of input dust
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ironsed   ! structure of input iron from sediment
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_hydrofe   ! structure of input iron from sediment
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dust    !: dust fields
@@ -74,6 +74,8 @@ CONTAINS
       REAL(wp) ::  zcoef, zwdust, zrivdin, zdustdep, zndep
       !
       CHARACTER (len=25) :: charout
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zw2d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_bc')
@@ -84,7 +86,9 @@ CONTAINS
       IF( ll_dust )  THEN
          !
          CALL fld_read( kt, 1, sf_dust )
-         dust(:,:) = MAX( rtrn, sf_dust(1)%fnow(:,:,1) )
+         DO_2D( 0, 0, 0, 0 )
+            dust(ji,jj) = MAX( rtrn, sf_dust(1)%fnow(ji,jj,1) )
+         END_2D
          !
          ! Iron solubilization of particles in the water column
          ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/d
@@ -107,9 +111,11 @@ CONTAINS
 
             IF( lk_iomput ) THEN
                 ! surface downward dust depo of iron
+                ALLOCATE( zw2d(A2D(0)) ) 
                 jl = n_trc_indsbc(jpfer)
-                CALL iom_put( "Irondep", ( rf_trsfac(jl) * sf_trcsbc(jl)%fnow(:,:,1) / rn_sbc_time ) * 1.e+3 * tmask(:,:,1) )
-
+                zw2d(A2D(0)) = rf_trsfac(jl) * ( sf_trcsbc(jl)%fnow(A2D(0),1) / rn_sbc_time ) * 1.e+3 * tmask(A2D(0),1)
+                CALL iom_put( "Irondep", zw2d )
+                DEALLOCATE( zw2d )
             ENDIF
 
          ENDIF
@@ -135,7 +141,10 @@ CONTAINS
          !
          IF( lk_iomput ) THEN
              ! dust concentration at surface
-             CALL iom_put( "pdust"  , dust(:,:) / ( wdust / rday ) * tmask(:,:,1) ) ! dust concentration at surface
+             ALLOCATE( zw2d(A2D(0)) )
+             zw2d(A2D(0)) = dust(A2D(0)) / ( wdust / rday ) * tmask(A2D(0),1)
+             CALL iom_put( "pdust", zw2d )
+             DEALLOCATE( zw2d )
          ENDIF
       ENDIF
 
@@ -190,34 +199,64 @@ CONTAINS
             tr(ji,jj,1,jpfer,Krhs) = tr(ji,jj,1,jpfer,Krhs) + zironice
          END_2D
          !
-         ! iron flux from ice
-         IF( lk_iomput ) &
-         & CALL iom_put( "Ironice", MAX( -0.99 * tr(:,:,1,jpfer,Kbb), (-1.*fmmflx(:,:)/1000._wp )*icefeinput*1.e+3*tmask(:,:,1)) )
+         IF( lk_iomput ) THEN
+            ! iron flux from ice
+            ALLOCATE( zw2d(A2D(0)) )
+            zw2d(A2D(0)) = MAX( -0.99 * tr(A2D(0),1,jpfer,Kbb), (-1.*fmmflx(A2D(0))/1000._wp )*icefeinput*1.e+3*tmask(A2D(0),1))
+            CALL iom_put( "Ironice", zw2d )
+            DEALLOCATE( zw2d )
+         ENDIF
          !
       ENDIF
 
       ! Add the external input of iron from sediment mobilization
       ! ------------------------------------------------------
       IF( ln_ironsed .AND. .NOT.lk_sed ) THEN
-          tr(:,:,:,jpfer,Krhs) = tr(:,:,:,jpfer,Krhs) + ironsed(:,:,:) * rfact
-          !
-          IF( lk_iomput )  CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) 
+          DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+             tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + ironsed(ji,jj,jk) * rfact
+         END_3D
+         !
+         IF( lk_iomput ) THEN
+            ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+            zw3d(A2D(0),1:jpkm1) = ironsed(A2D(0),1:jpkm1) * 1.e+3 * tmask(A2D(0),1:jpkm1)
+            CALL iom_put( "Ironsed", zw3d )
+            DEALLOCATE( zw3d )
+         ENDIF
       ENDIF
 
       ! Add the external input of iron from hydrothermal vents
       ! ------------------------------------------------------
       IF( ln_hydrofe ) THEN
          CALL fld_read( kt, 1, sf_hydrofe )
-         DO jk = 1, jpk
-            hydrofe(:,:,jk) = ( MAX( rtrn, sf_hydrofe(1)%fnow(:,:,jk) ) * hratio ) &
-              &              / ( e1e2t(:,:) * e3t(:,:,jk,Kmm) * ryyss + rtrn ) / 1000._wp &
-              &              * tmask(:,:,jk)
-         ENDDO
-                         tr(:,:,:,jpfer,Krhs) = tr(:,:,:,jpfer,Krhs) + hydrofe(:,:,:) * rfact
-         IF( ln_ligand ) tr(:,:,:,jplgw,Krhs) = tr(:,:,:,jplgw,Krhs) + ( hydrofe(:,:,:) * lgw_rath ) * rfact
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            hydrofe(ji,jj,jk) = ( MAX( rtrn, sf_hydrofe(1)%fnow(ji,jj,jk) ) * hratio ) &
+              &              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * ryyss + rtrn ) / 1000._wp &
+              &              * tmask(ji,jj,jk)
+         END_3D
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+              tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + hydrofe(ji,jj,jk) * rfact
+         END_3D
+         IF( ln_ligand ) THEN
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+               tr(ji,jj,jk,jplgw,Krhs) = tr(ji,jj,jk,jplgw,Krhs) + ( hydrofe(ji,jj,jk) * lgw_rath ) * rfact
+            END_3D
+         ENDIF
          !
-         IF( lk_iomput ) CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input
+         IF( lk_iomput ) THEN
+            ! hydrothermal iron input
+            ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+            zw3d(A2D(0),1:jpkm1) = hydrofe(A2D(0),1:jpkm1) * 1.e+3 * tmask(A2D(0),1:jpkm1)
+            CALL iom_put( "HYDR", zw3d )
+            DEALLOCATE( zw3d )
+         ENDIF
       ENDIF
+      !
+      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
+         WRITE(charout, FMT="('bc')")
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
+      ENDIF
+      !
       IF( ln_timing )  CALL timing_stop('p4z_bc')
       !
    END SUBROUTINE p4z_bc
@@ -303,7 +342,7 @@ CONTAINS
          IF(lwp) WRITE(numout,*) '    initialize dust input from atmosphere '
          IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
          !
-         ALLOCATE( dust(jpi,jpj) ) 
+         ALLOCATE( dust(A2D(0)) ) 
          !
          ALLOCATE( sf_dust(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
          IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_bc_init: unable to allocate sf_dust structure' )
@@ -321,7 +360,7 @@ CONTAINS
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '   ==>>>   ln_ironsed=T , computation of an island mask to enhance coastal supply of iron'
          !
-         ALLOCATE( ironsed(jpi,jpj,jpk) )    ! allocation
+         ALLOCATE( ironsed(A2D(0),jpk) )    ! allocation
          !
          CALL iom_open ( TRIM( sn_ironsed%clname ), numiron )
          ALLOCATE( zcmask(jpi,jpj,jpk) )
@@ -358,9 +397,9 @@ CONTAINS
          ! Coastal supply of iron
          ! -------------------------
          ironsed(:,:,jpk) = 0._wp
-         DO jk = 1, jpkm1
-            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_0(:,:,jk) * rday )
-         END DO
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            ironsed(ji,jj,jk) = sedfeinput * zcmask(ji,jj,jk) / ( e3t_0(ji,jj,jk) * rday )
+         END_3D
          DEALLOCATE( zcmask)
       ENDIF
       !
@@ -371,7 +410,7 @@ CONTAINS
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '   ==>>>   ln_hydrofe=T , Input of iron from hydrothermal vents'
          !
-         ALLOCATE( hydrofe(jpi,jpj,jpk) )    ! allocation
+         ALLOCATE( hydrofe(A2D(0),jpk) )    ! allocation
          !
          ALLOCATE( sf_hydrofe(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
          IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_bc_init: unable to allocate sf_hydro structure' )
diff --git a/src/TOP/PISCES/P4Z/p4zche.F90 b/src/TOP/PISCES/P4Z/p4zche.F90
index cacce4c2a..8d9f68fe0 100644
--- a/src/TOP/PISCES/P4Z/p4zche.F90
+++ b/src/TOP/PISCES/P4Z/p4zche.F90
@@ -168,9 +168,13 @@ CONTAINS
       ! practical salinity
       ! -------------------------------------------------------------
       IF (neos == -1) THEN
-         salinprac(:,:,:) = ts(:,:,:,jp_sal,Kmm) * 35.0 / 35.16504
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            salinprac(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * 35.0 / 35.16504
+         END_3D
       ELSE
-         salinprac(:,:,:) = ts(:,:,:,jp_sal,Kmm)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            salinprac(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm)
+         END_3D
       ENDIF
 
       !
@@ -797,17 +801,17 @@ CONTAINS
 
       ierr(:) = 0
 
-      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), STAT=ierr(1) )
+      ALLOCATE( sio3eq(A2D(0),jpk), fekeq(A2D(0),jpk), chemc(A2D(0),3), chemo2(A2D(0),jpk), STAT=ierr(1) )
 
-      ALLOCATE( akb3(jpi,jpj,jpk)     , tempis(jpi,jpj,jpk),       &
-         &      akw3(jpi,jpj,jpk)     , borat (jpi,jpj,jpk)  ,       &
-         &      aks3(jpi,jpj,jpk)     , akf3(jpi,jpj,jpk)    ,       &
-         &      ak1p3(jpi,jpj,jpk)    , ak2p3(jpi,jpj,jpk)   ,       &
-         &      ak3p3(jpi,jpj,jpk)    , aksi3(jpi,jpj,jpk)   ,       &
-         &      fluorid(jpi,jpj,jpk)  , sulfat(jpi,jpj,jpk)  ,       &
-         &      salinprac(jpi,jpj,jpk),                 STAT=ierr(2) )
+      ALLOCATE( akb3(A2D(0),jpk)     , tempis(A2D(0),jpk),       &
+         &      akw3(A2D(0),jpk)     , borat (A2D(0),jpk)  ,       &
+         &      aks3(A2D(0),jpk)     , akf3(A2D(0),jpk)    ,       &
+         &      ak1p3(A2D(0),jpk)    , ak2p3(A2D(0),jpk)   ,       &
+         &      ak3p3(A2D(0),jpk)    , aksi3(A2D(0),jpk)   ,       &
+         &      fluorid(A2D(0),jpk)  , sulfat(A2D(0),jpk)  ,       &
+         &      salinprac(A2D(0),jpk),                 STAT=ierr(2) )
 
-      ALLOCATE( fesol(jpi,jpj,jpk,5), STAT=ierr(3) )
+      ALLOCATE( fesol(A2D(0),jpk,5), STAT=ierr(3) )
 
       !* Variable for chemistry of the CO2 cycle
       p4z_che_alloc = MAXVAL( ierr )
diff --git a/src/TOP/PISCES/P4Z/p4zfechem.F90 b/src/TOP/PISCES/P4Z/p4zfechem.F90
index f98370b0a..8be5e7864 100644
--- a/src/TOP/PISCES/P4Z/p4zfechem.F90
+++ b/src/TOP/PISCES/P4Z/p4zfechem.F90
@@ -61,33 +61,42 @@ CONTAINS
       REAL(wp) ::   zprecip, zprecipno3,  zconsfe, za1
       REAL(wp) ::   zrfact2
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zTL1, zFe3, ztotlig, zfeprecip, zFeL1, zfecoll
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zcoll3d, zscav3d, zlcoll3d
+      REAL(wp), DIMENSION(A2D(0),jpk) ::   zTL1, zFe3, ztotlig, zfeprecip, zFeL1, zfecoll
+      REAL(wp), DIMENSION(A2D(0),jpk) ::   zcoll3d, zscav3d, zlcoll3d
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_fechem')
       !
-      zFe3     (:,:,jpk) = 0.
-      zFeL1    (:,:,jpk) = 0.
-      zTL1     (:,:,jpk) = 0.
-      zfeprecip(:,:,jpk) = 0.
-      zcoll3d  (:,:,jpk) = 0.
-      zscav3d  (:,:,jpk) = 0.
-      zlcoll3d (:,:,jpk) = 0.
-      zfecoll  (:,:,jpk) = 0.
-      xfecolagg(:,:,jpk) = 0.
-      xcoagfe  (:,:,jpk) = 0.
+!      zFe3     (:,:,jpk) = 0.
+!      zFeL1    (:,:,jpk) = 0.
+!      zTL1     (:,:,jpk) = 0.
+!      zfeprecip(:,:,jpk) = 0.
+!      zcoll3d  (:,:,jpk) = 0.
+!      zscav3d  (:,:,jpk) = 0.
+!      zlcoll3d (:,:,jpk) = 0.
+!      zfecoll  (:,:,jpk) = 0.
+!      xfecolagg(:,:,jpk) = 0.
+!      xcoagfe  (:,:,jpk) = 0.
       !
       ! Total ligand concentration : Ligands can be chosen to be constant or variable
       ! Parameterization from Pham and Ito (2018)
       ! -------------------------------------------------
-      xfecolagg(:,:,:) = ligand * 1E9 + MAX(0., chemo2(:,:,:) - tr(:,:,:,jpoxy,Kbb) ) / 400.E-6
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+         xfecolagg(ji,jj,jk) = ligand * 1E9 + MAX(0., chemo2(ji,jj,jk) - tr(ji,jj,jk,jpoxy,Kbb) ) / 400.E-6
+      END_3D
       IF( ln_ligvar ) THEN
-         ztotlig(:,:,:) =  0.09 * 0.667 * tr(:,:,:,jpdoc,Kbb) * 1E6 + xfecolagg(:,:,:)
-         ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
+            ztotlig(ji,jj,jk) =  0.09 * 0.667 * tr(ji,jj,jk,jpdoc,Kbb) * 1E6 + xfecolagg(ji,jj,jk)
+            ztotlig(ji,jj,jk) =  MIN( ztotlig(ji,jj,jk), 10. )
+         END_3D
       ELSE
-        IF( ln_ligand ) THEN  ;   ztotlig(:,:,:) = tr(:,:,:,jplgw,Kbb) * 1E9
-        ELSE                  ;   ztotlig(:,:,:) = ligand * 1E9 
+        IF( ln_ligand ) THEN  
+           DO_3D( 0, 0, 0, 0, 1, jpkm1)
+              ztotlig(ji,jj,jk) = tr(ji,jj,jk,jplgw,Kbb) * 1E9
+           END_3D
+        ELSE
+             ztotlig(:,:,:) = ligand * 1E9 
         ENDIF
       ENDIF
 
@@ -109,7 +118,9 @@ CONTAINS
           zFeL1(ji,jj,jk) = MAX( 0., tr(ji,jj,jk,jpfer,Kbb) - zFe3(ji,jj,jk) )
       END_3D
       !
-      plig(:,:,:) =  MAX( 0., ( zFeL1(:,:,:) / ( tr(:,:,:,jpfer,Kbb) + rtrn ) ) )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+         plig(ji,jj,jk) =  MAX( 0., ( zFeL1(ji,jj,jk) / ( tr(ji,jj,jk,jpfer,Kbb) + rtrn ) ) )
+      END_3D
       !
       zdust = 0.         ! if no dust available
 
@@ -125,12 +136,21 @@ CONTAINS
       ! to coagulate
       ! ----------------------------------------------------------------------
       IF (ln_ligand) THEN
-         zfecoll(:,:,:) = 0.5 * zFeL1(:,:,:) * MAX(0., tr(:,:,:,jplgw,Kbb) - xfecolagg(:,:,:) * 1.0E-9 ) / ( tr(:,:,:,jplgw,Kbb) + rtrn ) 
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
+            zfecoll(ji,jj,jk) = 0.5 * zFeL1(ji,jj,jk) * MAX(0., tr(ji,jj,jk,jplgw,Kbb) - xfecolagg(ji,jj,jk) * 1.0E-9 ) &
+                &                / ( tr(ji,jj,jk,jplgw,Kbb) + rtrn ) 
+         END_3D
       ELSE
          IF (ln_ligvar) THEN
-            zfecoll(:,:,:) = 0.5 * zFeL1(:,:,:) * MAX(0., tr(:,:,:,jplgw,Kbb) - xfecolagg(:,:,:) * 1.0E-9 ) / ( tr(:,:,:,jplgw,Kbb) + rtrn )   
+            DO_3D( 0, 0, 0, 0, 1, jpkm1)
+            !Ch Buguuuugggggggg : tr(ji,jj,jk,jplgw ?????????
+               zfecoll(ji,jj,jk) = 0.5 * zFeL1(ji,jj,jk) * MAX(0., tr(ji,jj,jk,jplgw,Kbb) - xfecolagg(ji,jj,jk) * 1.0E-9 ) &
+                &                   / ( tr(ji,jj,jk,jplgw,Kbb) + rtrn )   
+            END_3D
          ELSE
-            zfecoll(:,:,:) = 0.5 * zFeL1(:,:,:)
+            DO_3D( 0, 0, 0, 0, 1, jpkm1)
+               zfecoll(ji,jj,jk) = 0.5 * zFeL1(ji,jj,jk)
+            END_3D
          ENDIF
       ENDIF
 
@@ -199,44 +219,59 @@ CONTAINS
       !
       !  Define the bioavailable fraction of iron
       !  ----------------------------------------
-      biron(:,:,:) = tr(:,:,:,jpfer,Kbb) 
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+         biron(ji,jj,jk) = tr(ji,jj,jk,jpfer,Kbb) 
+      END_3D
       !
       !  Output of some diagnostics variables
       !     ---------------------------------
-      IF( lk_iomput .AND. knt == nrdttrc ) THEN
-         zrfact2 = 1.e3 * rfact2r  ! conversion from mol/L/timestep into mol/m3/s
-         !
-         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
+      IF( lk_iomput .AND.  knt == nrdttrc ) THEN
+        !
+        ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+        zrfact2 = 1.e3 * rfact2r  ! conversion from mol/L/timestep into mol/m3/s
+        !
+        IF( iom_use ( "Fe3" ) ) THEN   ! Fe3+
+           zw3d(A2D(0),1:jpkm1) = zFe3(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "Fe3", zw3d )
+        ENDIF
+        !
+        IF( iom_use ( "FeL1" ) ) THEN   ! FeL1
+           zw3d(A2D(0),1:jpkm1) = zFeL1(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "FeL1", zw3d )
+        ENDIF
+        !
+        IF( iom_use ( "TL1" ) ) THEN   ! TL1
+           zw3d(A2D(0),1:jpkm1) = zTL1(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "TL1", zw3d )
+        ENDIF
+        !
+        IF( iom_use ( "Totlig" ) ) THEN   ! Totlig
+           zw3d(A2D(0),1:jpkm1) = ztotlig(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "Totlig", zw3d )
+        ENDIF
+        !
+        IF( iom_use ( "Biron" ) ) THEN   ! biron
+           zw3d(A2D(0),1:jpkm1) = biron(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "Biron", zw3d )
+        ENDIF
+        !
+        IF( iom_use ( "FESCAV" ) ) THEN   ! FESCAV
+           zw3d(A2D(0),1:jpkm1) = zscav3d(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1) * zrfact2
+          CALL iom_put( "FESCAV", zw3d )
+        ENDIF
+        !
+        IF( iom_use ( "FECOLL" ) ) THEN   ! FECOLL
+           zw3d(A2D(0),1:jpkm1) = zcoll3d(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1) * zrfact2
+          CALL iom_put( "FECOLL", zw3d )
+        ENDIF
+        !
+        IF( iom_use ( "FEPREC" ) ) THEN   ! FEPREC
+           zw3d(A2D(0),1:jpkm1) = zfeprecip(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1) * zrfact2
+          CALL iom_put( "FEPREC", zw3d )
+        ENDIF
+        !
+        DEALLOCATE( zw3d )
+        !
       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 bd632bab8..7902a63b0 100644
--- a/src/TOP/PISCES/P4Z/p4zflx.F90
+++ b/src/TOP/PISCES/P4Z/p4zflx.F90
@@ -109,7 +109,11 @@ CONTAINS
          satmco2(:,:) = atcco2 
       ENDIF
 
-      IF( l_co2cpl )   satmco2(:,:) = atm_co2(:,:)
+      IF( l_co2cpl ) THEN
+         DO_2D( 0, 0, 0, 0 )
+            satmco2(ji,jj) = atm_co2(ji,jj)
+         END_2D
+      ENDIF
 
       DO_2D( 0, 0, 0, 0 )
          ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
@@ -172,11 +176,15 @@ CONTAINS
       END_2D
 
       IF( iom_use("tcflx") .OR. iom_use("tcflxcum") .OR. kt == nitrst   &
-         &                 .OR. (ln_check_mass .AND. kt == nitend) )    &
-         t_oce_co2_flx  = glob_sum( 'p4zflx', oce_co2(:,:) * e1e2t(:,:) * 1000. )                    !  Total Flux of Carbon
-      t_oce_co2_flx_cum = t_oce_co2_flx_cum + t_oce_co2_flx       !  Cumulative Total Flux of Carbon
-!      t_atm_co2_flx     = glob_sum( 'p4zflx', satmco2(:,:) * e1e2t(:,:) )       ! Total atmospheric pCO2
-      t_atm_co2_flx     =  atcco2      ! Total atmospheric pCO2
+         &                 .OR. (ln_check_mass .AND. kt == nitend) )  THEN
+         ALLOCATE( zw2d(A2D(0)) )
+         zw2d(A2D(0)) = oce_co2(A2D(0)) * e1e2t(A2D(0)) * 1000._wp
+         t_oce_co2_flx  = glob_sum( 'p4zflx',  zw2d(:,:) )                    !  Total Flux of Carbon
+         t_oce_co2_flx_cum = t_oce_co2_flx_cum + t_oce_co2_flx       !  Cumulative Total Flux of Carbon
+!        t_atm_co2_flx     = glob_sum( 'p4zflx', satmco2(:,:) * e1e2t(:,:) )       ! Total atmospheric pCO2
+         t_atm_co2_flx     =  atcco2      ! Total atmospheric pCO2
+         DEALLOCATE( zw2d )
+      ENDIF
  
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('flx ')")
@@ -186,18 +194,21 @@ CONTAINS
 
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
         !
-        ALLOCATE( zw2d(A2D(0)) )  ;  zw2d(A2D(0)) = 0._wp
+        ALLOCATE( zw2d(A2D(0)) ) 
         !
         IF( iom_use ( "AtmCo2" ) )  THEN  ! Atmospheric CO2 concentration
-          CALL iom_put( "AtmCo2"  , satmco2(:,:) * tmask(:,:,1) ) 
+          zw2d(A2D(0)) = satmco2(A2D(0)) * tmask(A2D(0),1)
+          CALL iom_put( "AtmCo2", zw2d )
         ENDIF
         !
         IF( iom_use ( "Cflx" ) ) THEN 
-          CALL iom_put( "Cflx"    , oce_co2(:,:) * 1000. ) 
+          zw2d(A2D(0)) = oce_co2(A2D(0)) * 1000._wp
+          CALL iom_put( "Cflx", zw2d )
         ENDIF
         !
         IF( iom_use ( "Dpo2" ) ) THEN
-          CALL iom_put( "Dpo2"    , ( atcox * patm(:,:) - atcox * tr(:,:,1,jpoxy,Kbb) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) )
+          zw2d(A2D(0)) =  ( atcox * patm(A2D(0)) - atcox * tr(A2D(0),1,jpoxy,Kbb) / ( chemo2(A2D(0),1) + rtrn ) ) * tmask(A2D(0),1) 
+          CALL iom_put( "Dpo2", zw2d )
         ENDIF
         !
         IF( iom_use ( "tcflx" ) ) THEN  ! molC/s
@@ -303,7 +314,7 @@ CONTAINS
          IF(lwp) WRITE(numout,*) '    Spatialized Atmospheric pCO2 from an external file'
       ENDIF
       !
-      oce_co2(:,:)  = 0._wp                ! Initialization of Flux of Carbon
+!      oce_co2(:,:)  = 0._wp                ! Initialization of Flux of Carbon
       t_oce_co2_flx = 0._wp
       t_atm_co2_flx = 0._wp
       !
@@ -324,6 +335,7 @@ CONTAINS
       CHARACTER(len=100) ::   cn_dir      ! Root directory for location of ssr files
       TYPE(FLD_N)        ::   sn_patm     ! informations about the fields to be read
       TYPE(FLD_N)        ::   sn_atmco2   ! informations about the fields to be read
+      INTEGER  :: ji, jj
       !!
       NAMELIST/nampisatm/ ln_presatm, ln_presatmco2, sn_patm, sn_atmco2, cn_dir
       !!----------------------------------------------------------------------
@@ -373,12 +385,16 @@ CONTAINS
       !
       IF( ln_presatm ) THEN
          CALL fld_read( kt, 1, sf_patm )               !* input Patm provided at kt + 1/2
-         patm(:,:) = sf_patm(1)%fnow(:,:,1)/101325.0     ! atmospheric pressure
+         DO_2D( 0, 0, 0, 0 )
+            patm(ji,jj) = sf_patm(1)%fnow(ji,jj,1)/101325.0     ! atmospheric pressure
+         END_2D
       ENDIF
       !
       IF( ln_presatmco2 ) THEN
          CALL fld_read( kt, 1, sf_atmco2 )               !* input atmco2 provided at kt + 1/2
-         satmco2(:,:) = sf_atmco2(1)%fnow(:,:,1)                        ! atmospheric pressure
+         DO_2D( 0, 0, 0, 0 )
+            satmco2(ji,jj) = sf_atmco2(1)%fnow(ji,jj,1)                        ! atmospheric pressure
+         END_2D
       ELSE
          satmco2(:,:) = atcco2    ! Initialize atmco2 if no reading from a file
       ENDIF
@@ -390,7 +406,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_flx_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( satmco2(jpi,jpj), patm(jpi,jpj), STAT=p4z_flx_alloc )
+      ALLOCATE( satmco2(A2D(0)), patm(A2D(0)), STAT=p4z_flx_alloc )
       !
       IF( p4z_flx_alloc /= 0 )   CALL ctl_stop( 'STOP', 'p4z_flx_alloc : failed to allocate arrays' )
       !
diff --git a/src/TOP/PISCES/P4Z/p4zint.F90 b/src/TOP/PISCES/P4Z/p4zint.F90
index 3b19ae86d..c0243fc4c 100644
--- a/src/TOP/PISCES/P4Z/p4zint.F90
+++ b/src/TOP/PISCES/P4Z/p4zint.F90
@@ -44,10 +44,12 @@ CONTAINS
       !
       ! Computation of phyto and zoo metabolic rate
       ! -------------------------------------------
-      ! Generic temperature dependence (Eppley, 1972)
-      tgfunc (:,:,:) = EXP( 0.0631 * ts(:,:,:,jp_tem,Kmm) )
-      ! Temperature dependence of mesozooplankton (Buitenhuis et al. (2005))
-      tgfunc2(:,:,:) = EXP( 0.0761 * ts(:,:,:,jp_tem,Kmm) )
+      DO_3D( 0, 0, 0, 0, 1, jpk )
+         ! Generic temperature dependence (Eppley, 1972)
+         tgfunc (ji,jj,jk) = EXP( 0.0631 * ts(ji,jj,jk,jp_tem,Kmm) )
+         ! Temperature dependence of mesozooplankton (Buitenhuis et al. (2005))
+         tgfunc2(ji,jj,jk) = EXP( 0.0761 * ts(ji,jj,jk,jp_tem,Kmm) )
+      END_3D
 
 
       ! Computation of the silicon dependant half saturation  constant for silica uptake
@@ -73,7 +75,7 @@ CONTAINS
       zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
 
       ! day length in hours
-      strn(:,:) = 0.
+!      strn(:,:) = 0.
       DO_2D( 0, 0, 0, 0 )
          zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
          zargu = MAX( -1., MIN(  1., zargu ) )
diff --git a/src/TOP/PISCES/P4Z/p4zlim.F90 b/src/TOP/PISCES/P4Z/p4zlim.F90
index 93a988c71..b8aaeffeb 100644
--- a/src/TOP/PISCES/P4Z/p4zlim.F90
+++ b/src/TOP/PISCES/P4Z/p4zlim.F90
@@ -388,16 +388,16 @@ CONTAINS
       !!----------------------------------------------------------------------
 
       !*  Biological arrays for phytoplankton growth
-      ALLOCATE( xnanono3(jpi,jpj,jpk), xdiatno3(jpi,jpj,jpk),       &
-         &      xnanonh4(jpi,jpj,jpk), xdiatnh4(jpi,jpj,jpk),       &
-         &      xnanopo4(jpi,jpj,jpk), xdiatpo4(jpi,jpj,jpk),       &
-         &      xnanofer(jpi,jpj,jpk), xdiatfer(jpi,jpj,jpk),       &
-         &      xlimphy (jpi,jpj,jpk), xlimdia (jpi,jpj,jpk),       &
-         &      xlimnfe (jpi,jpj,jpk), xlimdfe (jpi,jpj,jpk),       &
-         &      xlimbac (jpi,jpj,jpk), xlimbacl(jpi,jpj,jpk),       &
-         &      concnfe (jpi,jpj,jpk), concdfe (jpi,jpj,jpk),       &
-         &      xqfuncfecn(jpi,jpj,jpk), xqfuncfecd(jpi,jpj,jpk),   &
-         &      xlimsi  (jpi,jpj,jpk), STAT=p4z_lim_alloc )
+      ALLOCATE( xnanono3(A2D(0),jpk), xdiatno3(A2D(0),jpk),       &
+         &      xnanonh4(A2D(0),jpk), xdiatnh4(A2D(0),jpk),       &
+         &      xnanopo4(A2D(0),jpk), xdiatpo4(A2D(0),jpk),       &
+         &      xnanofer(A2D(0),jpk), xdiatfer(A2D(0),jpk),       &
+         &      xlimphy (A2D(0),jpk), xlimdia (A2D(0),jpk),       &
+         &      xlimnfe (A2D(0),jpk), xlimdfe (A2D(0),jpk),       &
+         &      xlimbac (A2D(0),jpk), xlimbacl(A2D(0),jpk),       &
+         &      concnfe (A2D(0),jpk), concdfe (A2D(0),jpk),       &
+         &      xqfuncfecn(A2D(0),jpk), xqfuncfecd(A2D(0),jpk),   &
+         &      xlimsi  (A2D(0),jpk), STAT=p4z_lim_alloc )
       !
       IF( p4z_lim_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_lim_alloc : failed to allocate arrays.' )
       !
diff --git a/src/TOP/PISCES/P4Z/p4zlys.F90 b/src/TOP/PISCES/P4Z/p4zlys.F90
index 9450e39f3..d7559a945 100644
--- a/src/TOP/PISCES/P4Z/p4zlys.F90
+++ b/src/TOP/PISCES/P4Z/p4zlys.F90
@@ -69,13 +69,15 @@ CONTAINS
 !! CE later      REAL(wp), DIMENSION(A2D(0),jpk) ::   zco3, zcaldiss, zhinit, zhi, zco3sat
 !!  because of the routine solve_at_general in p4zche.F90
       REAL(wp), DIMENSION(A2D(0),jpk) :: zco3, zcaldiss, zco3sat
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhinit, zhi
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zhinit, zhi
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)  :: zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )  CALL timing_start('p4z_lys')
       !
-      zhinit  (:,:,:) = hi(:,:,:) / ( rhd(:,:,:) + 1._wp )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+         zhinit(ji,jj,jk) = hi(ji,jj,jk) / ( rhd(ji,jj,jk) + 1._wp )
+      END_3D
       !
       !     -------------------------------------------
       !     COMPUTE [CO3--] and [H+] CONCENTRATIONS
diff --git a/src/TOP/PISCES/P4Z/p4zmeso.F90 b/src/TOP/PISCES/P4Z/p4zmeso.F90
index 8dc46de78..cdebf9a28 100644
--- a/src/TOP/PISCES/P4Z/p4zmeso.F90
+++ b/src/TOP/PISCES/P4Z/p4zmeso.F90
@@ -522,7 +522,7 @@ CONTAINS
       INTEGER  :: ji, jj, jk
       !
       REAL(wp) :: ztotchl, z1dep
-      REAL(wp), DIMENSION(jpi,jpj) :: oxymoy, tempmoy, zdepmoy
+      REAL(wp), DIMENSION(A2D(0)) :: oxymoy, tempmoy, zdepmoy
 
       !!---------------------------------------------------------------------
       !
@@ -600,7 +600,7 @@ CONTAINS
       !!                     ***  ROUTINE p4z_meso_alloc  ***
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( depmig(jpi,jpj), kmig(jpi,jpj), STAT= p4z_meso_alloc  )
+      ALLOCATE( depmig(A2D(0)), kmig(A2D(0)), STAT= p4z_meso_alloc  )
       !
       IF( p4z_meso_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_meso_alloc : failed to allocate arrays.' )
       !
diff --git a/src/TOP/PISCES/P4Z/p4zopt.F90 b/src/TOP/PISCES/P4Z/p4zopt.F90
index 94ec459ad..5ec1ddae1 100644
--- a/src/TOP/PISCES/P4Z/p4zopt.F90
+++ b/src/TOP/PISCES/P4Z/p4zopt.F90
@@ -63,10 +63,12 @@ CONTAINS
       INTEGER  ::   irgb
       REAL(wp) ::   zchl
       REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zetmp5
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zqsr100, zqsr_corr
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpar, ze0, ze1, ze2, ze3
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zetmp5
+      REAL(wp), DIMENSION(A2D(0)    ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4
+      REAL(wp), DIMENSION(A2D(0)    ) :: zqsr100, zqsr_corr
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zpar, ze0, ze1, ze2, ze3
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zw2d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_opt')
@@ -116,36 +118,40 @@ CONTAINS
             ! not fully correct with LIM3 and SI3 but no information is 
             ! currently available to do a better job. SHould be improved in the 
             ! (near) future.
-            zqsr_corr(:,:) = qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
+            DO_2D( 0, 0, 0, 0 )
+               zqsr_corr(ji,jj) = qsr_mean(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
+            END_2D
             !
             CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )
             !
             ! Used PAR is computed for each phytoplankton species
             ! etot_ndcy is PAR at level jk averaged over 24h.
             ! Due to their size, they have different light absorption characteristics
-            DO jk = 1, nksr
-               etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
-            END DO
+            DO_3D( 0, 0, 0, 0, 1, nksr )
+               etot_ndcy(ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
+            END_3D
             !
             ! SW over the ice free zone of the grid cell. This assumes that
             ! SW is zero below sea ice which is a very crude assumption that is 
             ! not fully correct with LIM3 and SI3 but no information is 
             ! currently available to do a better job. SHould be improved in the 
             ! (near) future.
-            zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
+            DO_2D( 0, 0, 0, 0 )
+               zqsr_corr(ji,jj) = qsr(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
+            END_2D
             !
             CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3 )
             !
             ! Total PAR computation at level jk that includes the diurnal cycle
-            DO jk = 1, nksr
-               etot (:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
-               enano(:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
-               ediat(:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)
-            END DO
+            DO_3D( 0, 0, 0, 0, 1, nksr )
+               etot (ji,jj,jk) =         ze1(ji,jj,jk) +        ze2(ji,jj,jk) +        ze3(ji,jj,jk)
+               enano(ji,jj,jk) =  1.85 * ze1(ji,jj,jk) + 0.69 * ze2(ji,jj,jk) + 0.46 * ze3(ji,jj,jk)
+               ediat(ji,jj,jk) =  1.62 * ze1(ji,jj,jk) + 0.74 * ze2(ji,jj,jk) + 0.63 * ze3(ji,jj,jk)
+            END_3D
             IF( ln_p5z ) THEN
-               DO jk = 1, nksr
-                  epico  (:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
-               END DO
+               DO_3D( 0, 0, 0, 0, 1, nksr )
+                  epico(ji,jj,jk) =  1.94 * ze1(ji,jj,jk) + 0.66 * ze2(ji,jj,jk) + 0.4 * ze3(ji,jj,jk)
+               END_3D
             ENDIF
 
          ELSE ! No diurnal cycle in PISCES
@@ -157,22 +163,24 @@ CONTAINS
             ! not fully correct with LIM3 and SI3 but no information is 
             ! currently available to do a better job. SHould be improved in the 
             ! (near) future.
-            zqsr_corr(:,:) = qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
+            DO_2D( 0, 0, 0, 0 )
+               zqsr_corr(ji,jj) = qsr_mean(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
+            END_2D
             !
             CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 ) 
             !
             ! Used PAR is computed for each phytoplankton species
             ! etot_ndcy is PAR at level jk averaged over 24h.
             ! Due to their size, they have different light absorption characteristics
-            DO jk = 1, nksr      
-               etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
-               enano    (:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
-               ediat    (:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)
-            END DO
+            DO_3D( 0, 0, 0, 0, 1, nksr )
+               etot_ndcy(ji,jj,jk) =         ze1(ji,jj,jk) +        ze2(ji,jj,jk) +        ze3(ji,jj,jk)
+               enano    (ji,jj,jk) =  1.85 * ze1(ji,jj,jk) + 0.69 * ze2(ji,jj,jk) + 0.46 * ze3(ji,jj,jk)
+               ediat    (ji,jj,jk) =  1.62 * ze1(ji,jj,jk) + 0.74 * ze2(ji,jj,jk) + 0.63 * ze3(ji,jj,jk)
+            END_3D
             IF( ln_p5z ) THEN
-               DO jk = 1, nksr      
-                  epico  (:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
-               END DO
+               DO_3D( 0, 0, 0, 0, 1, nksr )
+                  epico(ji,jj,jk) =  1.94 * ze1(ji,jj,jk) + 0.66 * ze2(ji,jj,jk) + 0.4 * ze3(ji,jj,jk)
+               END_3D
             ENDIF
             !
             ! SW over the ice free zone of the grid cell. This assumes that
@@ -180,14 +188,16 @@ CONTAINS
             ! not fully correct with LIM3 and SI3 but no information is 
             ! currently available to do a better job. SHould be improved in the 
             ! (near) future.
-            zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
+            DO_2D( 0, 0, 0, 0 )
+               zqsr_corr(ji,jj) = qsr(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
+            END_2D
             !
             CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3 ) 
             !
             ! Total PAR computation at level jk that includes the diurnal cycle
-            DO jk = 1, nksr      
-               etot(:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
-            END DO
+            DO_3D( 0, 0, 0, 0, 1, nksr )
+               etot(ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
+            END_3D
          ENDIF
          !
       ELSE   ! no diurnal cycle
@@ -198,22 +208,24 @@ CONTAINS
          ! not fully correct with LIM3 and SI3 but no information is 
          ! currently available to do a better job. SHould be improved in the 
          ! (near) future.
-         zqsr_corr(:,:) = qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
+         DO_2D( 0, 0, 0, 0 )
+            zqsr_corr(ji,jj) = qsr(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
+         END_2D
          !
          CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100  ) 
          !
 
          ! Used PAR is computed for each phytoplankton species
          ! Due to their size, they have different light absorption characteristics
-         DO jk = 1, nksr      
-            etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)    ! Total PAR
-            enano(:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)  ! Nanophytoplankton
-            ediat(:,:,jk) =  1.62 * ze1(:,:,jk) + 0.74 * ze2(:,:,jk) + 0.63 * ze3(:,:,jk)  ! Diatoms
-         END DO
+         DO_3D( 0, 0, 0, 0, 1, nksr )
+            etot (ji,jj,jk) =         ze1(ji,jj,jk) +        ze2(ji,jj,jk) +        ze3(ji,jj,jk)
+            enano(ji,jj,jk) =  1.85 * ze1(ji,jj,jk) + 0.69 * ze2(ji,jj,jk) + 0.46 * ze3(ji,jj,jk)
+            ediat(ji,jj,jk) =  1.62 * ze1(ji,jj,jk) + 0.74 * ze2(ji,jj,jk) + 0.63 * ze3(ji,jj,jk)
+         END_3D
          IF( ln_p5z ) THEN
-            DO jk = 1, nksr      
-              epico(:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)  ! Picophytoplankton (PISCES-QUOTA)
-            END DO
+            DO_3D( 0, 0, 0, 0, 1, nksr )
+               epico(ji,jj,jk) =  1.94 * ze1(ji,jj,jk) + 0.66 * ze2(ji,jj,jk) + 0.4 * ze3(ji,jj,jk) ! Picophytoplankton (PISCES-QUOTA)
+            END_3D
          ENDIF
          etot_ndcy(:,:,:) =  etot(:,:,:) 
       ENDIF
@@ -224,10 +236,12 @@ CONTAINS
          !                                     !  ------------------------
          CALL p4z_opt_par( kt, Kmm, qsr, ze1, ze2, ze3, pe0=ze0 )
          !
-         etot3(:,:,1) =  qsr(:,:) * tmask(:,:,1)
-         DO jk = 2, nksr + 1
-            etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
-         END DO
+         DO_2D( 0, 0, 0, 0 )
+            etot3(ji,jj,1) =  qsr(ji,jj) * tmask(ji,jj,1)
+         END_2D
+         DO_3D( 0, 0, 0, 0, 2, nksr+1 )
+            etot3(ji,jj,jk) =  ( ze0(ji,jj,jk) + ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk) ) * tmask(ji,jj,jk)
+         END_3D
          !                                     !  ------------------------
       ENDIF
       
@@ -236,9 +250,11 @@ CONTAINS
       ! (1) The classical definition based on the relative threshold value
       ! (2) An alternative definition based on a absolute threshold value.
       ! -------------------------------------------------------------------
-      neln(:,:) = 1
-      heup   (:,:) = gdepw(:,:,2,Kmm)
-      heup_01(:,:) = gdepw(:,:,2,Kmm)
+      DO_2D( 0, 0, 0, 0 )
+         neln   (ji,jj) = 1
+         heup   (ji,jj) = gdepw(ji,jj,2,Kmm)
+         heup_01(ji,jj) = gdepw(ji,jj,2,Kmm)
+      END_2D
 
       DO_3D( 0, 0, 0, 0, 2, nksr)
         IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >=  zqsr100(ji,jj) )  THEN
@@ -306,7 +322,7 @@ CONTAINS
       !
       IF( ln_p5z ) THEN
          ! Picophytoplankton when using PISCES-QUOTA
-         ALLOCATE( zetmp5(jpi,jpj) )  ;   zetmp5 (:,:) = 0.e0
+         ALLOCATE( zetmp5(A2D(0)) )  ;   zetmp5 (:,:) = 0.e0
          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)
@@ -325,10 +341,18 @@ CONTAINS
       ENDIF
       !
       IF( lk_iomput .AND.  knt == nrdttrc ) THEN
-         CALL iom_put( "Heup" , heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
+         IF( iom_use( "Heup" ) ) THEN
+           ALLOCATE( zw2d(A2D(0)) ) 
+           zw2d(A2D(0)) = heup(A2D(0)) * tmask(A2D(0),1)
+           CALL iom_put( "Heup", zw2d )  ! Euphotic layer depth
+           DEALLOCATE( zw2d ) 
+        ENDIF
          IF( iom_use( "PAR" ) ) THEN
-            zpar(:,:,1) = zpar(:,:,1) * ( 1._wp - fr_i(:,:) )
-            CALL iom_put( "PAR", zpar(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
+           ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+            zw3d(A2D(0),1) = zpar(A2D(0),1) * ( 1._wp - fr_i(A2D(0)) ) * tmask(A2D(0),1)
+            zw3d(A2D(0),2:jpkm1) = zpar(A2D(0),2:jpkm1) * tmask(A2D(0),2:jpkm1)
+            CALL iom_put( "PAR", zw3d )  ! Photosynthetically Available Radiation
+            DEALLOCATE( zw3d ) 
          ENDIF
       ENDIF
       !
@@ -345,15 +369,15 @@ CONTAINS
       !!                for a given shortwave radiation
       !!
       !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in)              ::   kt                ! ocean time-step
-      INTEGER                         , INTENT(in)              ::   Kmm               ! ocean time-index
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   )           ::   pqsr              ! shortwave
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::   pe1 , pe2 , pe3   ! PAR ( R-G-B)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::   pe0               !
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out), OPTIONAL ::   pqsr100           !
+      INTEGER                        , INTENT(in)              ::   kt                ! ocean time-step
+      INTEGER                        , INTENT(in)              ::   Kmm               ! ocean time-index
+      REAL(wp), DIMENSION(A2D(0))    , INTENT(in   )           ::   pqsr              ! shortwave
+      REAL(wp), DIMENSION(A2D(0),jpk), INTENT(inout)           ::   pe1 , pe2 , pe3   ! PAR ( R-G-B)
+      REAL(wp), DIMENSION(A2D(0),jpk), INTENT(inout), OPTIONAL ::   pe0               !
+      REAL(wp), DIMENSION(A2D(0))    , INTENT(  out), OPTIONAL ::   pqsr100           !
       !
       INTEGER    ::   ji, jj, jk     ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj) ::  zqsr   ! shortwave
+      REAL(wp), DIMENSION(A2D(0)) ::  zqsr   ! shortwave
       !!----------------------------------------------------------------------
 
       !  Real shortwave
@@ -419,7 +443,9 @@ CONTAINS
       IF( ln_varpar ) THEN
          IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
             CALL fld_read( kt, 1, sf_par )
-            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0
+            DO_2D( 0, 0, 0, 0 )
+               par_varsw(ji,jj) = ( sf_par(1)%fnow(ji,jj,1) ) / 3.0
+            END_2D
          ENDIF
       ENDIF
       !
@@ -479,7 +505,7 @@ CONTAINS
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '   ==>>>   initialize variable par fraction (ln_varpar=T)'
          !
-         ALLOCATE( par_varsw(jpi,jpj) )
+         ALLOCATE( par_varsw(A2D(0)) )
          !
          ALLOCATE( sf_par(1), STAT=ierr )           !* allocate and fill sf_sst (forcing structure) with sn_sst
          IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_opt_init: unable to allocate sf_par structure' )
@@ -510,8 +536,7 @@ CONTAINS
       !!                     ***  ROUTINE p4z_opt_alloc  ***
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( ekb(jpi,jpj,jpk), ekr(jpi,jpj,jpk),  &
-                ekg(jpi,jpj,jpk), STAT= p4z_opt_alloc  ) 
+      ALLOCATE( ekb(A2D(0),jpk), ekr(A2D(0),jpk), ekg(A2D(0),jpk), STAT= p4z_opt_alloc  ) 
       !
       IF( p4z_opt_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_opt_alloc : failed to allocate arrays.' )
       !
diff --git a/src/TOP/PISCES/P4Z/p4zpoc.F90 b/src/TOP/PISCES/P4Z/p4zpoc.F90
index ab156817e..b96011e11 100644
--- a/src/TOP/PISCES/P4Z/p4zpoc.F90
+++ b/src/TOP/PISCES/P4Z/p4zpoc.F90
@@ -525,7 +525,7 @@ CONTAINS
       ! Discretization along the lability space
       ! ---------------------------------------
       !
-      ALLOCATE( alphan(jcpoc) , reminp(jcpoc) , alphap(jpi,jpj,jpk,jcpoc) )
+      ALLOCATE( alphan(jcpoc) , reminp(jcpoc) , alphap(A2D(0),jpk,jcpoc) )
       !
       IF (jcpoc > 1) THEN  ! Case when more than one lability class is used
          !
diff --git a/src/TOP/PISCES/P4Z/p4zprod.F90 b/src/TOP/PISCES/P4Z/p4zprod.F90
index f3189c38d..d6da5eee9 100644
--- a/src/TOP/PISCES/P4Z/p4zprod.F90
+++ b/src/TOP/PISCES/P4Z/p4zprod.F90
@@ -543,8 +543,8 @@ CONTAINS
       texcretd  = 1._wp - excretd
       tpp       = 0._wp
       ! CE not really needed ; tempory, shoub be removed when quotan( A2D(0),jpk ) 
-      quotan(:,:,:) = 0._wp 
-      quotad(:,:,:) = 0._wp 
+!      quotan(:,:,:) = 0._wp 
+!      quotad(:,:,:) = 0._wp 
       !
    END SUBROUTINE p4z_prod_init
 
@@ -553,7 +553,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_prod_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( quotan(jpi,jpj,jpk), quotad(jpi,jpj,jpk), STAT = p4z_prod_alloc )
+      ALLOCATE( quotan(A2D(0),jpk), quotad(A2D(0),jpk), STAT = p4z_prod_alloc )
       !
       IF( p4z_prod_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_prod_alloc : failed to allocate arrays.' )
       !
diff --git a/src/TOP/PISCES/P4Z/p4zsed.F90 b/src/TOP/PISCES/P4Z/p4zsed.F90
index 040b05059..6b4949b08 100644
--- a/src/TOP/PISCES/P4Z/p4zsed.F90
+++ b/src/TOP/PISCES/P4Z/p4zsed.F90
@@ -69,27 +69,26 @@ CONTAINS
       REAL(wp) ::  xdiano3, xdianh4
       !
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zdenit2d, zbureff, zwork
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zwsbio3, zwsbio4
-      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
+      REAL(wp), DIMENSION(A2D(0)) :: zdenit2d, zbureff, zwork
+      REAL(wp), DIMENSION(A2D(0)) :: zwsbio3, zwsbio4
+      REAL(wp), DIMENSION(A2D(0)) :: zsedcal, zsedsi, zsedc
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zsoufer, zlight
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )  CALL timing_start('p4z_sed')
       !
 
       ! Allocate temporary workspace
-      ALLOCATE( ztrpo4(jpi,jpj,jpk) )
-      IF( ln_p5z )    ALLOCATE( ztrdop(jpi,jpj,jpk) )
+      ALLOCATE( ztrpo4(A2D(0),jpk) )
+      IF( ln_p5z )    ALLOCATE( ztrdop(A2D(0),jpk) )
 
       zdenit2d(:,:) = 0.e0
       zbureff (:,:) = 0.e0
-      zwork   (:,:) = 0.e0
-      zsedsi  (:,:) = 0.e0
-      zsedcal (:,:) = 0.e0
-      zsedc   (:,:) = 0.e0
+ !     zwork   (:,:) = 0.e0
+ !     zsedsi  (:,:) = 0.e0
+ !     zsedcal (:,:) = 0.e0
+ !     zsedc   (:,:) = 0.e0
 
       ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments
       ! --------------------------------------------------------------------
@@ -219,10 +218,11 @@ CONTAINS
       ! Nitrogen fixation process
       ! Small source iron from particulate inorganic iron
       !-----------------------------------
-      DO jk = 1, jpkm1
-         zlight (:,:,jk) =  ( 1.- EXP( -etot_ndcy(:,:,jk) / diazolight ) ) * ( 1. - fr_i(:,:) ) 
-         zsoufer(:,:,jk) = zlight(:,:,jk) * 2E-11 / ( 2E-11 + biron(:,:,jk) )
-      ENDDO
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+         zlight (ji,jj,jk) =  ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) * ( 1. - fr_i(ji,jj) ) 
+         zsoufer(ji,jj,jk) = zlight(ji,jj,jk) * 2E-11 / ( 2E-11 + biron(ji,jj,jk) )
+      END_3D
+      !
       IF( ln_p4z ) THEN
          DO_3D( 0, 0, 0, 0, 1, jpkm1)
             !                      ! Potential nitrogen fixation dependant on temperature and iron
@@ -309,32 +309,18 @@ CONTAINS
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
           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(:,:,:) )
+            ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+            zw3d(A2D(0),1:jpkm1) = nitrpot(A2D(0),1:jpkm1) * nitrfix * rno3 * zfact * tmask(A2D(0),1:jpkm1)
+            CALL iom_put( "Nfix", zw3d )
+            DEALLOCATE( zw3d )
           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
+          IF( iom_use( "Sdenit" ) )  CALL iom_put( "Sdenit", sdenit(:,:)  * rno3 * zfact )
+          IF( iom_use( "SedCal" ) )  CALL iom_put( "SedCal", zsedcal(:,:)        * zfact )
+          IF( iom_use( "SedSi"  ) )  CALL iom_put( "SedSi" , zsedsi(:,:)         * zfact )
+          IF( iom_use( "SedC"   ) )  CALL iom_put( "SedC"  , zsedc(:,:)          * zfact )
           !
-          DEALLOCATE( zw2d )
       ENDIF
       !
       IF(sn_cfctl%l_prttrc) THEN  ! print mean trneds (USEd for debugging)
@@ -343,7 +329,8 @@ CONTAINS
          CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
-      IF( ln_p5z )    DEALLOCATE( ztrpo4, ztrdop )
+                    DEALLOCATE( ztrpo4 )
+      IF( ln_p5z )  DEALLOCATE( ztrdop )
       !
       IF( ln_timing )  CALL timing_stop('p4z_sed')
       !
@@ -390,7 +377,7 @@ CONTAINS
       !
       lk_sed = ln_sediment .AND. ln_sed_2way 
       !
-      nitrpot(:,:,jpk) = 0._wp   ! define last level for iom_put
+!      nitrpot(:,:,jpk) = 0._wp   ! define last level for iom_put
       !
    END SUBROUTINE p4z_sed_init
 
@@ -398,7 +385,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_sed_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( nitrpot(jpi,jpj,jpk), sdenit(jpi,jpj), STAT=p4z_sed_alloc )
+      ALLOCATE( nitrpot(A2D(0),jpk), sdenit(A2D(0)), STAT=p4z_sed_alloc )
       !
       IF( p4z_sed_alloc /= 0 )   CALL ctl_stop( 'STOP', 'p4z_sed_alloc: failed to allocate arrays' )
       !
diff --git a/src/TOP/PISCES/P4Z/p4zsink.F90 b/src/TOP/PISCES/P4Z/p4zsink.F90
index 48bcf3f93..a810685f2 100644
--- a/src/TOP/PISCES/P4Z/p4zsink.F90
+++ b/src/TOP/PISCES/P4Z/p4zsink.F90
@@ -68,6 +68,8 @@ CONTAINS
       INTEGER  ::   ji, jj, jk
       CHARACTER (len=25) :: charout
       REAL(wp) :: zmax, zfact
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zw2d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_sink')
@@ -129,54 +131,69 @@ CONTAINS
       ENDIF
 
      ! Total carbon export per year
-     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
-        &   t_oce_co2_exp = glob_sum( 'p4zsink', ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
+     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  THEN
+         ALLOCATE( zw2d(A2D(0)) )
+         zw2d(A2D(0)) = ( sinking(A2D(0),ik100) + sinking2(A2D(0),ik100) ) * e1e2t(A2D(0)) * tmask(A2D(0),1)
+         t_oce_co2_exp = glob_sum( 'p4zsink', zw2d(:,:) )
+         DEALLOCATE( zw2d )
+     ENDIF
      !
      IF( lk_iomput .AND.  knt == nrdttrc ) THEN
        zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
        !
+       ALLOCATE( zw3d(A2D(0),jpk), zw2d(A2D(0)) )  ;  zw3d(A2D(0),jpk) = 0._wp
+       !
        IF( iom_use ( "EPC100" ) ) THEN ! Export of carbon at 100m
-         CALL iom_put( "EPC100"  , ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ) 
+         zw2d(A2D(0)) = ( sinking(A2D(0),ik100) + sinking2(A2D(0),ik100) ) * zfact * tmask(A2D(0),1)
+         CALL iom_put( "EPC100",  zw2d ) 
        ENDIF
        !
        IF( iom_use ( "EPFE100" ) ) THEN ! Export of iron at 100m
-         CALL iom_put( "EPFE100" , ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ) 
+         zw2d(A2D(0)) = ( sinkfer(A2D(0),ik100) + sinkfer2(A2D(0),ik100) ) * zfact * tmask(A2D(0),1)
+         CALL iom_put( "EPFE100",  zw2d ) 
        ENDIF
        !
        IF( iom_use ( "EPCAL100" ) ) THEN ! Export of calcite at 100m
-         CALL iom_put( "EPCAL100", sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ) 
+         zw2d(A2D(0)) = sinkcal(A2D(0),ik100)  * zfact * tmask(A2D(0),1)
+         CALL iom_put( "EPCAL100",  zw2d ) 
        ENDIF
        !
        IF( iom_use ( "EPSI100" ) ) THEN ! Export of bigenic silica at 100m
-         CALL iom_put( "EPSI100" , sinksil(:,:,ik100) * zfact * tmask(:,:,1) ) 
+         zw2d(A2D(0)) = sinksil(A2D(0),ik100) * zfact * tmask(A2D(0),1)
+         CALL iom_put( "EPSI100",  zw2d ) 
        ENDIF
        !
        IF( iom_use ( "EXPC" ) ) THEN ! Export of carbon in the water column
-         CALL iom_put( "EXPC"    , ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ) 
+         zw3d(A2D(0),1:jpkm1) = ( sinking(A2D(0),1:jpkm1) + sinking2(A2D(0),1:jpkm1) ) * zfact * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "EXPC",  zw3d ) 
        ENDIF
        !
        IF( iom_use ( "EXPFE" ) ) THEN ! Export of iron
-         CALL iom_put( "EXPFE"   , ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ) 
+         zw3d(A2D(0),1:jpkm1) = ( sinkfer(A2D(0),1:jpkm1) + sinkfer2(A2D(0),1:jpkm1) ) * zfact * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "EXPFE",  zw3d ) 
        ENDIF
        !
        IF( iom_use ( "EXPCAL" ) ) THEN ! Export of calcite
-         CALL iom_put( "EXPCAL"  , sinkcal(:,:,:) * zfact * tmask(:,:,:) ) 
+         zw3d(A2D(0),1:jpkm1) = sinkcal(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "EXPCAL",  zw3d ) 
        ENDIF
        !
        IF( iom_use ( "EXPSI" ) ) THEN ! Export of bigenic silica
-         CALL iom_put( "EXPSI"   , sinksil(:,:,:) * zfact * tmask(:,:,:) ) 
+         zw3d(A2D(0),1:jpkm1) = sinksil(A2D(0),1:jpkm1) * zfact * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "EXPSI",  zw3d ) 
        ENDIF
        !
-       IF( iom_use ( "EXPSI" ) ) THEN ! molC/s
+       IF( iom_use ( "tcexp" ) ) THEN ! molC/s
          CALL iom_put( "tcexp"   , t_oce_co2_exp * zfact ) 
        ENDIF
        ! 
+       DEALLOCATE( zw3d, zw2d )
       ENDIF
       !
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('sink')")
          CALL prt_ctl_info( charout, cdcomp = 'top' )
-         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Kbb), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
       IF( ln_timing )   CALL timing_stop('p4z_sink')
@@ -218,13 +235,13 @@ CONTAINS
       !
       ierr(:) = 0
       !
-      ALLOCATE( sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk)                    ,     &                
-         &      sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk)                    ,     &                
-         &      sinkfer2(jpi,jpj,jpk)                                           ,     &                
-         &      sinkfer(jpi,jpj,jpk)                                            , STAT=ierr(1) )                
+      ALLOCATE( sinking(A2D(0),jpk) , sinking2(A2D(0),jpk)                    ,     &                
+         &      sinkcal(A2D(0),jpk) , sinksil (A2D(0),jpk)                    ,     &                
+         &      sinkfer2(A2D(0),jpk)                                           ,     &                
+         &      sinkfer(A2D(0),jpk)                                            , STAT=ierr(1) )                
          !
-      IF( ln_p5z    ) ALLOCATE( sinkingn(jpi,jpj,jpk), sinking2n(jpi,jpj,jpk)   ,     &
-         &                      sinkingp(jpi,jpj,jpk), sinking2p(jpi,jpj,jpk)   , STAT=ierr(2) )
+      IF( ln_p5z    ) ALLOCATE( sinkingn(A2D(0),jpk), sinking2n(A2D(0),jpk)   ,     &
+         &                      sinkingp(A2D(0),jpk), sinking2p(A2D(0),jpk)   , STAT=ierr(2) )
       !
       p4z_sink_alloc = MAXVAL( ierr )
       IF( p4z_sink_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_sink_alloc : failed to allocate arrays.' )
diff --git a/src/TOP/PISCES/P4Z/p4zsms.F90 b/src/TOP/PISCES/P4Z/p4zsms.F90
index 4f1c7122b..dc2b32ecc 100644
--- a/src/TOP/PISCES/P4Z/p4zsms.F90
+++ b/src/TOP/PISCES/P4Z/p4zsms.F90
@@ -157,45 +157,56 @@ CONTAINS
         IF(  iom_use( 'INTdtAlk' ) .OR. iom_use( 'INTdtDIC' ) .OR. iom_use( 'INTdtFer' ) .OR.  &
           &  iom_use( 'INTdtDIN' ) .OR. iom_use( 'INTdtDIP' ) .OR. iom_use( 'INTdtSil' ) )  THEN
           !
-          ALLOCATE( zw3d(jpi,jpj,jpk), zw2d(jpi,jpj) )
-          zw3d(:,:,jpk) = 0.
-          DO jk = 1, jpkm1
-              zw3d(:,:,jk) = xnegtr(:,:,jk) * xfact * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
-          ENDDO
+          ALLOCATE( zw3d(A2D(0),jpk), zw2d(A2D(0)) )
+          DO_3D( 0, 0, 0, 0, 1, jpkm1)
+             zw3d(ji,jj,jk) = xnegtr(ji,jj,jk) * xfact * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+          END_3D
           !
           zw2d(:,:) = 0.
           DO jk = 1, jpkm1
-             zw2d(:,:) = zw2d(:,:) + zw3d(:,:,jk) * tr(:,:,jk,jptal,Krhs)
+             DO_2D( 0, 0, 0, 0 )
+                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) *  tr(ji,jj,jk,jptal,Krhs) 
+             END_2D
           ENDDO
           CALL iom_put( 'INTdtAlk', zw2d )
           !
           zw2d(:,:) = 0.
           DO jk = 1, jpkm1
-             zw2d(:,:) = zw2d(:,:) + zw3d(:,:,jk) * tr(:,:,jk,jpdic,Krhs)
+             DO_2D( 0, 0, 0, 0 )
+                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) *  tr(ji,jj,jk,jpdic,Krhs) 
+             END_2D
           ENDDO
           CALL iom_put( 'INTdtDIC', zw2d )
           !
           zw2d(:,:) = 0.
           DO jk = 1, jpkm1
-             zw2d(:,:) = zw2d(:,:) + zw3d(:,:,jk) * rno3 * ( tr(:,:,jk,jpno3,Krhs) + tr(:,:,jk,jpnh4,Krhs) )
+             DO_2D( 0, 0, 0, 0 )
+                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) * rno3 * ( tr(ji,jj,jk,jpno3,Krhs) + tr(ji,jj,jk,jpnh4,Krhs) ) 
+             END_2D
           ENDDO
           CALL iom_put( 'INTdtDIN', zw2d )
           !
           zw2d(:,:) = 0.
           DO jk = 1, jpkm1
-             zw2d(:,:) = zw2d(:,:) + zw3d(:,:,jk) * po4r * tr(:,:,jk,jppo4,Krhs)
+             DO_2D( 0, 0, 0, 0 )
+                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) * po4r * tr(ji,jj,jk,jppo4,Krhs) 
+             END_2D
           ENDDO
           CALL iom_put( 'INTdtDIP', zw2d )
           !
           zw2d(:,:) = 0.
           DO jk = 1, jpkm1
-             zw2d(:,:) = zw2d(:,:) + zw3d(:,:,jk) * tr(:,:,jk,jpfer,Krhs)
+             DO_2D( 0, 0, 0, 0 )
+                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) *  tr(ji,jj,jk,jpfer,Krhs) 
+             END_2D
           ENDDO
           CALL iom_put( 'INTdtFer', zw2d )
           !
           zw2d(:,:) = 0.
           DO jk = 1, jpkm1
-             zw2d(:,:) = zw2d(:,:) + zw3d(:,:,jk) * tr(:,:,jk,jpsil,Krhs)
+             DO_2D( 0, 0, 0, 0 )
+                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) *  tr(ji,jj,jk,jpsil,Krhs) 
+             END_2D
           ENDDO
           CALL iom_put( 'INTdtSil', zw2d )
           !
@@ -522,8 +533,9 @@ CONTAINS
       INTEGER, INTENT( in ) ::   Kmm     ! time level indices
       REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot
       CHARACTER(LEN=100)   ::   cltxt
-      INTEGER :: jk
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwork
+      INTEGER :: ji, jj, jk
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zw2d
       !!----------------------------------------------------------------------
       !
       IF( kt == nittrc000 ) THEN 
@@ -542,82 +554,113 @@ CONTAINS
 
       ! Compute the budget of NO3
       IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
+         ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
          IF( ln_p4z ) THEN
-            zwork(:,:,:) =    tr(:,:,:,jpno3,Kmm) + tr(:,:,:,jpnh4,Kmm)                      &
-               &          +   tr(:,:,:,jpphy,Kmm) + tr(:,:,:,jpdia,Kmm)                      &
-               &          +   tr(:,:,:,jppoc,Kmm) + tr(:,:,:,jpgoc,Kmm)  + tr(:,:,:,jpdoc,Kmm)  &        
-               &          +   tr(:,:,:,jpzoo,Kmm) + tr(:,:,:,jpmes,Kmm) 
+            DO_3D( 0, 0, 0, 0, 1, jpk)
+            zw3d(ji,jj,jk)  =  ( tr(ji,jj,jk,jpno3,Kmm) + tr(ji,jj,jk,jpnh4,Kmm)                      &
+               &             +   tr(ji,jj,jk,jpphy,Kmm) + tr(ji,jj,jk,jpdia,Kmm)                      &
+               &             +   tr(ji,jj,jk,jppoc,Kmm) + tr(ji,jj,jk,jpgoc,Kmm)  + tr(ji,jj,jk,jpdoc,Kmm)  &        
+               &             +   tr(ji,jj,jk,jpzoo,Kmm) + tr(ji,jj,jk,jpmes,Kmm)  ) * cvol(ji,jj,jk)
+            END_3D
         ELSE
-            zwork(:,:,:) =    tr(:,:,:,jpno3,Kmm) + tr(:,:,:,jpnh4,Kmm) + tr(:,:,:,jpnph,Kmm)   &
-               &          +   tr(:,:,:,jpndi,Kmm) + tr(:,:,:,jpnpi,Kmm)                      & 
-               &          +   tr(:,:,:,jppon,Kmm) + tr(:,:,:,jpgon,Kmm) + tr(:,:,:,jpdon,Kmm)   &
-               &          + ( tr(:,:,:,jpzoo,Kmm) + tr(:,:,:,jpmes,Kmm) ) * no3rat3 
+            DO_3D( 0, 0, 0, 0, 1, jpk)
+            zw3d(ji,jj,jk) =    ( tr(ji,jj,jk,jpno3,Kmm) + tr(ji,jj,jk,jpnh4,Kmm) + tr(ji,jj,jk,jpnph,Kmm)   &
+               &             +   tr(ji,jj,jk,jpndi,Kmm) + tr(ji,jj,jk,jpnpi,Kmm)                      & 
+               &             +   tr(ji,jj,jk,jppon,Kmm) + tr(ji,jj,jk,jpgon,Kmm) + tr(ji,jj,jk,jpdon,Kmm)   &
+               &             + ( tr(ji,jj,jk,jpzoo,Kmm) + tr(ji,jj,jk,jpmes,Kmm) ) * no3rat3 ) * cvol(ji,jj,jk)
+            END_3D
         ENDIF
         !
-        no3budget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )  
+        no3budget = glob_sum( 'p4zsms', zw3d(:,:,:)  )  
         no3budget = no3budget / areatot
         CALL iom_put( "pno3tot", no3budget )
+        DEALLOCATE( zw3d )
       ENDIF
       !
       ! Compute the budget of PO4
       IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
+         ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
          IF( ln_p4z ) THEN
-            zwork(:,:,:) =    tr(:,:,:,jppo4,Kmm)                                         &
-               &          +   tr(:,:,:,jpphy,Kmm) + tr(:,:,:,jpdia,Kmm)                      &
-               &          +   tr(:,:,:,jppoc,Kmm) + tr(:,:,:,jpgoc,Kmm)  + tr(:,:,:,jpdoc,Kmm)  &        
-               &          +   tr(:,:,:,jpzoo,Kmm) + tr(:,:,:,jpmes,Kmm) 
-        ELSE
-            zwork(:,:,:) =    tr(:,:,:,jppo4,Kmm) + tr(:,:,:,jppph,Kmm)                      &
-               &          +   tr(:,:,:,jppdi,Kmm) + tr(:,:,:,jpppi,Kmm)                      & 
-               &          +   tr(:,:,:,jppop,Kmm) + tr(:,:,:,jpgop,Kmm) + tr(:,:,:,jpdop,Kmm)   &
-               &          + ( tr(:,:,:,jpzoo,Kmm) + tr(:,:,:,jpmes,Kmm) ) * po4rat3 
+            DO_3D( 0, 0, 0, 0, 1, jpk)
+               zw3d(ji,jj,jk) =    ( tr(ji,jj,jk,jppo4,Kmm)                                         &
+                  &             +   tr(ji,jj,jk,jpphy,Kmm) + tr(ji,jj,jk,jpdia,Kmm)                      &
+                  &             +   tr(ji,jj,jk,jppoc,Kmm) + tr(ji,jj,jk,jpgoc,Kmm)  + tr(ji,jj,jk,jpdoc,Kmm)  &        
+                  &             +   tr(ji,jj,jk,jpzoo,Kmm) + tr(ji,jj,jk,jpmes,Kmm) ) * cvol(ji,jj,jk)
+            END_3D
+         ELSE
+            DO_3D( 0, 0, 0, 0, 1, jpk)
+               zw3d(ji,jj,jk) =    ( tr(ji,jj,jk,jppo4,Kmm) + tr(ji,jj,jk,jppph,Kmm)                      &
+                  &             +   tr(ji,jj,jk,jppdi,Kmm) + tr(ji,jj,jk,jpppi,Kmm)                      & 
+                  &             +   tr(ji,jj,jk,jppop,Kmm) + tr(ji,jj,jk,jpgop,Kmm) + tr(ji,jj,jk,jpdop,Kmm)   &
+                  &             + ( tr(ji,jj,jk,jpzoo,Kmm) + tr(ji,jj,jk,jpmes,Kmm) ) * po4rat3 ) * cvol(ji,jj,jk)
+            END_3D
         ENDIF
         !
-        po4budget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )  
+        po4budget = glob_sum( 'p4zsms', zw3d(:,:,:)  )  
         po4budget = po4budget / areatot
         CALL iom_put( "ppo4tot", po4budget )
+        DEALLOCATE( zw3d )
       ENDIF
       !
       ! Compute the budget of SiO3
       IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
-         zwork(:,:,:) =  tr(:,:,:,jpsil,Kmm) + tr(:,:,:,jpgsi,Kmm) + tr(:,:,:,jpdsi,Kmm) 
+         ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+         DO_3D( 0, 0, 0, 0, 1, jpk)
+            zw3d(ji,jj,jk) =  ( tr(ji,jj,jk,jpsil,Kmm) + tr(ji,jj,jk,jpgsi,Kmm) + tr(ji,jj,jk,jpdsi,Kmm) ) * cvol(ji,jj,jk)
+         END_3D
          !
-         silbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )  
+         silbudget = glob_sum( 'p4zsms', zw3d(:,:,:)  )  
          silbudget = silbudget / areatot
          CALL iom_put( "psiltot", silbudget )
+         DEALLOCATE( zw3d )
       ENDIF
       !
       IF( iom_use( "palktot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
-         zwork(:,:,:) =  tr(:,:,:,jpno3,Kmm) * rno3 + tr(:,:,:,jptal,Kmm) + tr(:,:,:,jpcal,Kmm) * 2.              
+         ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+         DO_3D( 0, 0, 0, 0, 1, jpk)
+            zw3d(ji,jj,jk) =  ( tr(ji,jj,jk,jpno3,Kmm) * rno3 + tr(ji,jj,jk,jptal,Kmm) + tr(ji,jj,jk,jpcal,Kmm) * 2. ) * cvol(ji,jj,jk)             
+         END_3D
          !
-         alkbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )         !
+         alkbudget = glob_sum( 'p4zsms', zw3d(:,:,:)  )         !
          alkbudget = alkbudget / areatot
          CALL iom_put( "palktot", alkbudget )
+         DEALLOCATE( zw3d )
       ENDIF
       !
       ! Compute the budget of Iron
       IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
-         zwork(:,:,:) =   tr(:,:,:,jpfer,Kmm) + tr(:,:,:,jpnfe,Kmm) + tr(:,:,:,jpdfe,Kmm)   &
-            &         +   tr(:,:,:,jpbfe,Kmm) + tr(:,:,:,jpsfe,Kmm)                      &
-            &         + ( tr(:,:,:,jpzoo,Kmm) * feratz + tr(:,:,:,jpmes,Kmm) ) * feratm   
+         ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+         DO_3D( 0, 0, 0, 0, 1, jpk)
+            zw3d(ji,jj,jk) =   ( tr(ji,jj,jk,jpfer,Kmm) + tr(ji,jj,jk,jpnfe,Kmm) + tr(ji,jj,jk,jpdfe,Kmm)   &
+               &            +    tr(ji,jj,jk,jpbfe,Kmm) + tr(ji,jj,jk,jpsfe,Kmm)                      &
+               &            +    tr(ji,jj,jk,jpzoo,Kmm) * feratz + tr(ji,jj,jk,jpmes,Kmm) * feratm ) * cvol(ji,jj,jk)  
+         END_3D
          !
-         ferbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )  
+         ferbudget = glob_sum( 'p4zsms', zw3d(:,:,:) )  
          ferbudget = ferbudget / areatot
          CALL iom_put( "pfertot", ferbudget )
+         DEALLOCATE( zw3d )
       ENDIF
       !
       ! Global budget of N SMS : denitrification in the water column and in the sediment
       !                          nitrogen fixation by the diazotrophs
       ! --------------------------------------------------------------------------------
       IF( iom_use( "tnfix" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
-         znitrpottot  = glob_sum ( 'p4zsms', nitrpot(:,:,:) * nitrfix * cvol(:,:,:) )
+         ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+         zw3d(A2D(0),1:jpkm1) =  nitrpot(A2D(0),1:jpkm1) * nitrfix * cvol(A2D(0),1:jpkm1)
+         znitrpottot  = glob_sum ( 'p4zsms',  zw3d)
          CALL iom_put( "tnfix"  , znitrpottot * xfact3 )  ! Global  nitrogen fixation molC/l  to molN/m3 
+         DEALLOCATE( zw3d )
       ENDIF
       !
       IF( iom_use( "tdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
-         zrdenittot = glob_sum ( 'p4zsms', denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
-         zsdenittot = glob_sum ( 'p4zsms', sdenit(:,:) * e1e2t(:,:) * tmask(:,:,1) )
+         ALLOCATE( zw3d(A2D(0),jpk), zw2d(A2D(0)) )  ;  zw3d(A2D(0),jpk) = 0._wp
+         zw3d(A2D(0),1:jpkm1) =  denitr(A2D(0),1:jpkm1) *  rdenit * xnegtr(A2D(0),1:jpkm1) * cvol(A2D(0),1:jpkm1)
+         zw2d(A2D(0))         =  sdenit(A2D(0)) *  e1e2t(A2D(0)) * tmask(A2D(0),1)
+         zrdenittot = glob_sum ( 'p4zsms', zw3d )
+         zsdenittot = glob_sum ( 'p4zsms', Zw2d )
          CALL iom_put( "tdenit" , ( zrdenittot + zsdenittot ) * xfact3 )  ! Total denitrification molC/l to molN/m3 
+         DEALLOCATE( zw3d, zw2d )
       ENDIF
       !
       IF( ln_check_mass .AND. kt == nitend ) THEN   ! Compute the budget of NO3, ALK, Si, Fer
diff --git a/src/TOP/PISCES/P4Z/p5zlim.F90 b/src/TOP/PISCES/P4Z/p5zlim.F90
index dffbe0384..ff8040915 100644
--- a/src/TOP/PISCES/P4Z/p5zlim.F90
+++ b/src/TOP/PISCES/P4Z/p5zlim.F90
@@ -692,24 +692,24 @@ CONTAINS
       ierr(:) = 0
       !
       !*  Biological arrays for phytoplankton growth
-      ALLOCATE( xpicono3(jpi,jpj,jpk), xpiconh4(jpi,jpj,jpk),       &
-         &      xpicopo4(jpi,jpj,jpk), xpicodop(jpi,jpj,jpk),       &
-         &      xnanodop(jpi,jpj,jpk), xdiatdop(jpi,jpj,jpk),       &
-         &      xpicofer(jpi,jpj,jpk), xlimpfe (jpi,jpj,jpk),       &
-         &      fvnuptk (jpi,jpj,jpk), fvduptk (jpi,jpj,jpk),       &
-         &      xlimphys(jpi,jpj,jpk), xlimdias(jpi,jpj,jpk),       &
-         &      xlimnpp (jpi,jpj,jpk), xlimnpn (jpi,jpj,jpk),       &
-         &      xlimnpd (jpi,jpj,jpk),                              &
-         &      xlimpics(jpi,jpj,jpk), xqfuncfecp(jpi,jpj,jpk),     &
-         &      fvpuptk (jpi,jpj,jpk), xlimpic (jpi,jpj,jpk),    STAT=ierr(1) )
+      ALLOCATE( xpicono3(A2D(0),jpk), xpiconh4(A2D(0),jpk),       &
+         &      xpicopo4(A2D(0),jpk), xpicodop(A2D(0),jpk),       &
+         &      xnanodop(A2D(0),jpk), xdiatdop(A2D(0),jpk),       &
+         &      xpicofer(A2D(0),jpk), xlimpfe (A2D(0),jpk),       &
+         &      fvnuptk (A2D(0),jpk), fvduptk (A2D(0),jpk),       &
+         &      xlimphys(A2D(0),jpk), xlimdias(A2D(0),jpk),       &
+         &      xlimnpp (A2D(0),jpk), xlimnpn (A2D(0),jpk),       &
+         &      xlimnpd (A2D(0),jpk),                              &
+         &      xlimpics(A2D(0),jpk), xqfuncfecp(A2D(0),jpk),     &
+         &      fvpuptk (A2D(0),jpk), xlimpic (A2D(0),jpk),    STAT=ierr(1) )
          !
       !*  Minimum/maximum quotas of phytoplankton
-      ALLOCATE( xqnnmin (jpi,jpj,jpk), xqnnmax(jpi,jpj,jpk),       &
-         &      xqpnmin (jpi,jpj,jpk), xqpnmax(jpi,jpj,jpk),       &
-         &      xqnpmin (jpi,jpj,jpk), xqnpmax(jpi,jpj,jpk),       &
-         &      xqppmin (jpi,jpj,jpk), xqppmax(jpi,jpj,jpk),       &
-         &      xqndmin (jpi,jpj,jpk), xqndmax(jpi,jpj,jpk),       &
-         &      xqpdmin (jpi,jpj,jpk), xqpdmax(jpi,jpj,jpk),     STAT=ierr(2) )
+      ALLOCATE( xqnnmin (A2D(0),jpk), xqnnmax(A2D(0),jpk),       &
+         &      xqpnmin (A2D(0),jpk), xqpnmax(A2D(0),jpk),       &
+         &      xqnpmin (A2D(0),jpk), xqnpmax(A2D(0),jpk),       &
+         &      xqppmin (A2D(0),jpk), xqppmax(A2D(0),jpk),       &
+         &      xqndmin (A2D(0),jpk), xqndmax(A2D(0),jpk),       &
+         &      xqpdmin (A2D(0),jpk), xqpdmax(A2D(0),jpk),     STAT=ierr(2) )
          !
       p5z_lim_alloc = MAXVAL( ierr )
       !
diff --git a/src/TOP/PISCES/P4Z/p5zmeso.F90 b/src/TOP/PISCES/P4Z/p5zmeso.F90
index 7dc423073..ce7eefc81 100644
--- a/src/TOP/PISCES/P4Z/p5zmeso.F90
+++ b/src/TOP/PISCES/P4Z/p5zmeso.F90
@@ -708,7 +708,7 @@ CONTAINS
       !!                     ***  ROUTINE p5z_meso_alloc  ***
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( depmig(jpi,jpj), kmig(jpi,jpj), STAT= p5z_meso_alloc  )
+      ALLOCATE( depmig(A2D(0)), kmig(A2D(0)), STAT= p5z_meso_alloc  )
       !
       IF( p5z_meso_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p5z_meso_alloc : failed to allocate arrays.' )
       !
diff --git a/src/TOP/PISCES/P4Z/p5zprod.F90 b/src/TOP/PISCES/P4Z/p5zprod.F90
index bc54b6e30..4430242e5 100644
--- a/src/TOP/PISCES/P4Z/p5zprod.F90
+++ b/src/TOP/PISCES/P4Z/p5zprod.F90
@@ -628,15 +628,18 @@ CONTAINS
        ENDIF
        !
        IF( iom_use ( "MunetP" ) ) THEN ! Realized growth rate for picophyto
-       CALL iom_put( "MunetP"  , ( tr(:,:,:,jppic,Krhs)/rfact2/(tr(:,:,:,jppic,Kbb)+ rtrn ) * tmask(:,:,:)) ) 
+         zw3d(A2D(0),1:jpkm1) = tr(A2D(0),1:jpkm1,jppic,Krhs)/rfact2/(tr(A2D(0),1:jpkm1,jppic,Kbb) + rtrn ) * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "MunetP", zw3d )
        ENDIF
        !
        IF( iom_use ( "MunetN" ) ) THEN ! Realized growth rate for nanophyto
-       CALL iom_put( "MunetN"  , ( tr(:,:,:,jpphy,Krhs)/rfact2/(tr(:,:,:,jpphy,Kbb)+ rtrn ) * tmask(:,:,:)) ) 
+         zw3d(A2D(0),1:jpkm1) = tr(A2D(0),1:jpkm1,jpphy,Krhs)/rfact2/(tr(A2D(0),1:jpkm1,jpphy,Kbb) + rtrn ) * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "MunetN", zw3d )
        ENDIF
        !
        IF( iom_use ( "MunetD" ) ) THEN ! Realized growth rate for diatoms
-       CALL iom_put( "MunetD"  , ( tr(:,:,:,jpdia,Krhs)/rfact2/(tr(:,:,:,jpdia,Kbb)+ rtrn ) * tmask(:,:,:)) ) 
+         zw3d(A2D(0),1:jpkm1) = tr(A2D(0),1:jpkm1,jpdia,Krhs)/rfact2/(tr(A2D(0),1:jpkm1,jpdia,Kbb) + rtrn ) * tmask(A2D(0),1:jpkm1)
+         CALL iom_put( "MunetD", zw3d )
        ENDIF
        !
        IF( iom_use ( "TPP" ) ) THEN   ! total primary production
@@ -726,7 +729,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE p5z_prod_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( zdaylen(jpi,jpj), STAT = p5z_prod_alloc )
+      ALLOCATE( zdaylen(A2D(0)), STAT = p5z_prod_alloc )
       !
       IF( p5z_prod_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p5z_prod_alloc : failed to allocate arrays.' )
       !
diff --git a/src/TOP/PISCES/sms_pisces.F90 b/src/TOP/PISCES/sms_pisces.F90
index baf78a5d7..e4862e949 100644
--- a/src/TOP/PISCES/sms_pisces.F90
+++ b/src/TOP/PISCES/sms_pisces.F90
@@ -142,52 +142,52 @@ CONTAINS
       !!----------------------------------------------------------------------
       ierr(:) = 0
       !*  Biological fluxes for light : shared variables for pisces & lobster
-      ALLOCATE( xksi(jpi,jpj), strn(jpi,jpj),  STAT=ierr(1) )
+      ALLOCATE( xksi(A2D(0)), strn(A2D(0)),  STAT=ierr(1) )
 
       IF( ln_p4z .OR. ln_p5z ) THEN
 
          !* Optics
-         ALLOCATE(  enano(jpi,jpj,jpk) , ediat(jpi,jpj,jpk) ,   &
-           &        enanom(jpi,jpj,jpk), ediatm(jpi,jpj,jpk),   &
-           &        emoy(jpi,jpj,jpk)  , etotm(jpi,jpj,jpk),   STAT=ierr(2) )
+         ALLOCATE(  enano(A2D(0),jpk) , ediat(A2D(0),jpk) ,   &
+           &        enanom(A2D(0),jpk), ediatm(A2D(0),jpk),   &
+           &        emoy(A2D(0),jpk)  , etotm(A2D(0),jpk),   STAT=ierr(2) )
 
          !* Biological SMS
-         ALLOCATE( xksimax(jpi,jpj)  , biron(jpi,jpj,jpk)      ,  STAT=ierr(3) )
+         ALLOCATE( xksimax(A2D(0))  , biron(A2D(0),jpk)      ,  STAT=ierr(3) )
 
          ! Biological SMS
-         ALLOCATE( xfracal  (jpi,jpj,jpk), orem    (jpi,jpj,jpk),  &
-            &      nitrfac  (jpi,jpj,jpk), nitrfac2(jpi,jpj,jpk),  &
-            &      prodcal  (jpi,jpj,jpk), xdiss   (jpi,jpj,jpk),  &
-            &      prodpoc  (jpi,jpj,jpk), conspoc (jpi,jpj,jpk),  &
-            &      prodgoc  (jpi,jpj,jpk), consgoc (jpi,jpj,jpk),  &
-            &      blim     (jpi,jpj,jpk), consfe3 (jpi,jpj,jpk),  &
-            &      xfecolagg(jpi,jpj,jpk), xcoagfe (jpi,jpj,jpk), STAT=ierr(4) )
+         ALLOCATE( xfracal  (A2D(0),jpk), orem    (A2D(0),jpk),  &
+            &      nitrfac  (A2D(0),jpk), nitrfac2(A2D(0),jpk),  &
+            &      prodcal  (A2D(0),jpk), xdiss   (A2D(0),jpk),  &
+            &      prodpoc  (A2D(0),jpk), conspoc (A2D(0),jpk),  &
+            &      prodgoc  (A2D(0),jpk), consgoc (A2D(0),jpk),  &
+            &      blim     (A2D(0),jpk), consfe3 (A2D(0),jpk),  &
+            &      xfecolagg(A2D(0),jpk), xcoagfe (A2D(0),jpk), STAT=ierr(4) )
 
          !* Carbonate chemistry
-         ALLOCATE( ak13  (jpi,jpj,jpk)  ,                          &
-            &      ak23(jpi,jpj,jpk)    , aksp  (jpi,jpj,jpk) ,    &
-            &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,    &
-            &      aphscale(jpi,jpj,jpk),                         STAT=ierr(5) )
+         ALLOCATE( ak13(A2D(0),jpk),                          &
+            &      ak23(A2D(0),jpk), aksp  (A2D(0),jpk) ,    &
+            &      hi  (A2D(0),jpk), excess(A2D(0),jpk) ,    &
+            &      aphscale(A2D(0),jpk),                         STAT=ierr(5) )
          !
          !* Temperature dependency of SMS terms
-         ALLOCATE( tgfunc (jpi,jpj,jpk) , tgfunc2(jpi,jpj,jpk),   STAT=ierr(6) )
+         ALLOCATE( tgfunc (A2D(0),jpk) , tgfunc2(A2D(0),jpk),   STAT=ierr(6) )
          !
          !* Sinking speed
-         ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4 (jpi,jpj,jpk),   STAT=ierr(7) )
+         ALLOCATE( wsbio3 (A2D(0),jpk) , wsbio4 (A2D(0),jpk),   STAT=ierr(7) )
 
          !*  Size of phytoplankton cells
-         ALLOCATE( sizen (jpi,jpj,jpk), sized (jpi,jpj,jpk),        &
-           &       sizena(jpi,jpj,jpk), sizeda(jpi,jpj,jpk),      STAT=ierr(8) )
+         ALLOCATE( sizen (A2D(0),jpk), sized (A2D(0),jpk),        &
+           &       sizena(A2D(0),jpk), sizeda(A2D(0),jpk),      STAT=ierr(8) )
          ! 
-         ALLOCATE( plig(jpi,jpj,jpk)  ,                           STAT=ierr(9) )
+         ALLOCATE( plig(A2D(0),jpk)  ,                           STAT=ierr(9) )
       ENDIF
       !
       IF( ln_p5z ) THEN
          ! PISCES-QUOTA specific part      
-         ALLOCATE( epico(jpi,jpj,jpk)   , epicom(jpi,jpj,jpk) ,   STAT=ierr(10) ) 
+         ALLOCATE( epico(A2D(0),jpk)   , epicom(A2D(0),jpk) ,   STAT=ierr(10) ) 
 
          !*  Size of phytoplankton cells
-         ALLOCATE( sizep(jpi,jpj,jpk), sizepa(jpi,jpj,jpk),       STAT=ierr(11) )
+         ALLOCATE( sizep(A2D(0),jpk), sizepa(A2D(0),jpk),       STAT=ierr(11) )
       ENDIF
       !
       sms_pisces_alloc = MAXVAL( ierr )
diff --git a/src/TOP/TRP/trcdmp.F90 b/src/TOP/TRP/trcdmp.F90
index 2b11fd629..03535a63c 100644
--- a/src/TOP/TRP/trcdmp.F90
+++ b/src/TOP/TRP/trcdmp.F90
@@ -57,7 +57,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE trc_dmp_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( restotr(jpi,jpj,jpk) , STAT=trc_dmp_alloc )
+      ALLOCATE( restotr(A2D(0),jpk) , STAT=trc_dmp_alloc )
       !
       IF( trc_dmp_alloc /= 0 )   CALL ctl_warn('trc_dmp_alloc: failed to allocate array')
       !
diff --git a/src/TOP/TRP/trcsink.F90 b/src/TOP/TRP/trcsink.F90
index 5cee9b0ba..572f501b7 100644
--- a/src/TOP/TRP/trcsink.F90
+++ b/src/TOP/TRP/trcsink.F90
@@ -50,12 +50,12 @@ CONTAINS
       INTEGER , INTENT(in)  :: Kbb, Kmm
       INTEGER , INTENT(in)  :: jp_tra    ! tracer index index      
       REAL(wp), INTENT(in)  :: rsfact    ! time step duration
-      REAL(wp), INTENT(in)   , DIMENSION(jpi,jpj,jpk) :: pwsink
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: psinkflx
+      REAL(wp), INTENT(in)   , DIMENSION(A2D(0),jpk) :: pwsink
+      REAL(wp), INTENT(inout), DIMENSION(A2D(0),jpk) :: psinkflx
       INTEGER  ::   ji, jj, jk
-      INTEGER, DIMENSION(jpi, jpj) ::   iiter
+      INTEGER, DIMENSION(A2D(0)) ::   iiter
       REAL(wp) ::   zfact, zwsmax, zmax
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwsink
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zwsink
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('trc_sink')
@@ -73,7 +73,7 @@ CONTAINS
       IF( nitermax == 1 ) THEN
          iiter(:,:) = 1
       ELSE
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
             iiter(ji,jj) = 1
             DO jk = 1, jpkm1
                IF( tmask(ji,jj,jk) == 1.0 ) THEN
@@ -85,7 +85,7 @@ CONTAINS
          iiter(:,:) = MIN( iiter(:,:), nitermax )
       ENDIF
 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          IF( tmask(ji,jj,jk) == 1.0 ) THEN
            zwsmax = 0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact
            zwsink(ji,jj,jk) = MIN( pwsink(ji,jj,jk), zwsmax * REAL( iiter(ji,jj), wp ) )
@@ -121,23 +121,25 @@ CONTAINS
       INTEGER,  INTENT(in   )                         ::   Kbb, Kmm  ! time level indices
       INTEGER,  INTENT(in   )                         ::   jp_tra    ! tracer index index      
       REAL(wp), INTENT(in   )                         ::   rsfact    ! duration of time step
-      INTEGER,  INTENT(in   ), DIMENSION(jpi,jpj)     ::   kiter     ! number of iterations for time-splitting 
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
+      INTEGER,  INTENT(in   ), DIMENSION(A2D(0))     ::   kiter     ! number of iterations for time-splitting 
+      REAL(wp), INTENT(in   ), DIMENSION(A2D(0),jpk) ::   pwsink    ! sinking speed
+      REAL(wp), INTENT(inout), DIMENSION(A2D(0),jpk) ::   psinkflx  ! sinking fluxe
       !
       INTEGER  ::   ji, jj, jk, jn, jt
       REAL(wp) ::   zigma,zew,zign, zflx, zstep
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb, psinking 
+      REAL(wp), DIMENSION(A2D(0),jpk) :: ztraz, zakz, zwsink2, ztrb, psinking 
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('trc_sink2')
       !
-      DO jk = 1, jpkm1
-         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) 
-      END DO
-      zwsink2(:,:,1) = 0.e0
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         zwsink2(ji,jj,jk+1) = -pwsink(ji,jj,jk) / rday * tmask(ji,jj,jk+1) 
+      END_3D
+      DO_2D( 0, 0, 0, 0 )
+         zwsink2(ji,jj,1) = 0.e0
+      END_2D
 
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          ! Vertical advective flux
          zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2.
          DO jt = 1, kiter(ji,jj)
diff --git a/src/TOP/trc.F90 b/src/TOP/trc.F90
index 6958508cc..4836efb7e 100644
--- a/src/TOP/trc.F90
+++ b/src/TOP/trc.F90
@@ -161,11 +161,11 @@ CONTAINS
       ALLOCATE( tr(jpi,jpj,jpk,jptra,jpt)                                         ,       &  
          &      gtru (jpi,jpj,jptra) , gtrv (jpi,jpj,jptra)                       ,       &
          &      gtrui(jpi,jpj,jptra) , gtrvi(jpi,jpj,jptra)                       ,       &
-         &      trc_i(jpi,jpj,jptra)  , trc_o(jpi,jpj,jptra)                        ,       &
+         &      trc_i(A2D(0),jptra)  , trc_o(A2D(0),jptra)                        ,       &
          &      trc_ice_ratio(jptra) , trc_ice_prescr(jptra) , cn_trc_o(jptra)    ,       &
-         &      neln(jpi,jpj)         , heup(jpi,jpj)         , heup_01(jpi,jpj)     ,       &
-         &      etot(jpi,jpj,jpk)     , etot_ndcy(jpi,jpj,jpk)                      ,       &
-         &      sbc_trc_b(jpi,jpj,jptra), sbc_trc(jpi,jpj,jptra)                    ,       &  
+         &      neln(A2D(0))         , heup(A2D(0))         , heup_01(A2D(0))     ,       &
+         &      etot(A2D(0),jpk)     , etot_ndcy(A2D(0),jpk)                      ,       &
+         &      sbc_trc_b(A2D(0),jptra), sbc_trc(A2D(0),jptra)                    ,       &  
          &      cvol(jpi,jpj,jpk)    , trai(jptra)                                ,       &
          &      ctrcnm(jptra)        , ctrcln(jptra)         , ctrcun(jptra)      ,       &
          &      ln_trc_ini(jptra)    ,                                                    &
diff --git a/src/TOP/trcini.F90 b/src/TOP/trcini.F90
index 394f1cd16..f5275ead3 100644
--- a/src/TOP/trcini.F90
+++ b/src/TOP/trcini.F90
@@ -32,6 +32,8 @@ MODULE trcini
    
    PUBLIC   trc_init   ! called by opa
 
+      !! * Substitutions
+#  include "do_loop_substitute.h90"
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -93,9 +95,8 @@ CONTAINS
       !! ** Purpose :      passive tracers inventories at initialsation phase
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   Kmm    ! time level index
-      INTEGER             ::  jk, jn  ! dummy loop indices
+      INTEGER             ::  ji, jj, jk, jn  ! dummy loop indices
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) :: zzmsk
       !!----------------------------------------------------------------------
       !
       IF(lwp) WRITE(numout,*)
diff --git a/src/TOP/trcopt.F90 b/src/TOP/trcopt.F90
index d07b4500d..c7a013167 100644
--- a/src/TOP/trcopt.F90
+++ b/src/TOP/trcopt.F90
@@ -58,12 +58,14 @@ CONTAINS
       INTEGER, INTENT(in) ::   kt, knt   ! ocean time step
       INTEGER, INTENT(in) ::   Kbb, Kmm  ! time level indices
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   zchl  ! chlorophyll field
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out) :: ze1, ze2, ze3 ! PAR for individual wavelength
+      REAL(wp), DIMENSION(A2D(0),jpk), INTENT(out) :: ze1, ze2, ze3 ! PAR for individual wavelength
       !
       INTEGER  ::   ji, jj, jk, irgb
       REAL(wp) ::   ztmp
-      REAL(wp), DIMENSION(jpi,jpj    ) :: parsw, zqsr100, zqsr_corr
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze0
+      REAL(wp), DIMENSION(A2D(0)    ) :: parsw, zqsr100, zqsr_corr
+      REAL(wp), DIMENSION(A2D(0),jpk) :: ze0
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zw2d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('trc_opt')
@@ -85,7 +87,7 @@ CONTAINS
 
       !     Attenuation coef. function of Chlorophyll and wavelength (RGB)
       !     --------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          ztmp = ( zchl(ji,jj,jk) + rtrn ) * 1.e6
          ztmp = MIN(  10. , MAX( 0.05, ztmp )  )
          irgb = NINT( 41 + 20.* LOG10( ztmp ) + rtrn )
@@ -99,54 +101,63 @@ CONTAINS
       !     -----------------------------------------------
       IF( ln_qsr_bio ) THEN
          !
-         zqsr_corr(:,:) = parsw(:,:) * qsr(:,:)
+         DO_2D( 0, 0, 0, 0 )
+            zqsr_corr(ji,jj) = parsw(ji,jj) * qsr(ji,jj)
+         END_2D
          !
-         ze0(:,:,1) = (1._wp - 3._wp * parsw(:,:)) * qsr(:,:)  !  ( 1 - 3 * alpha ) * q
+         DO_2D( 0, 0, 0, 0 )
+            ze0(ji,jj,1) = (1._wp - 3._wp * parsw(ji,jj)) * qsr(ji,jj)  !  ( 1 - 3 * alpha ) * q
+         END_2D
          ze1(:,:,1) = zqsr_corr(:,:)
          ze2(:,:,1) = zqsr_corr(:,:)
          ze3(:,:,1) = zqsr_corr(:,:)
          !
-         DO jk = 2, nksrp + 1
-            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-                  ze0(ji,jj,jk) = ze0(ji,jj,jk-1) * EXP( -e3t(ji,jj,jk-1,Kmm) * (1. / rn_si0) )
-                  ze1(ji,jj,jk) = ze1(ji,jj,jk-1) * EXP( -ekb  (ji,jj,jk-1 )        )
-                  ze2(ji,jj,jk) = ze2(ji,jj,jk-1) * EXP( -ekg  (ji,jj,jk-1 )        )
-                  ze3(ji,jj,jk) = ze3(ji,jj,jk-1) * EXP( -ekr  (ji,jj,jk-1 )        )
-            END_2D
-         END DO
+         DO_3D( 0, 0, 0, 0, 2, nksrp + 1 )
+            ze0(ji,jj,jk) = ze0(ji,jj,jk-1) * EXP( -e3t(ji,jj,jk-1,Kmm) * (1. / rn_si0) )
+            ze1(ji,jj,jk) = ze1(ji,jj,jk-1) * EXP( -ekb  (ji,jj,jk-1 )        )
+            ze2(ji,jj,jk) = ze2(ji,jj,jk-1) * EXP( -ekg  (ji,jj,jk-1 )        )
+            ze3(ji,jj,jk) = ze3(ji,jj,jk-1) * EXP( -ekr  (ji,jj,jk-1 )        )
+         END_3D
          !
-         etot3(:,:,1) = qsr(:,:) * tmask(:,:,1)
-         DO jk = 2, nksrp + 1
-            etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
-         END DO
+         DO_2D( 0, 0, 0, 0 )
+            etot3(ji,jj,1) =  qsr(ji,jj) * tmask(ji,jj,1)
+         END_2D
+         DO_3D( 0, 0, 0, 0, 2, nksrp+1 )
+            etot3(ji,jj,jk) =  ( ze0(ji,jj,jk) + ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk) ) * tmask(ji,jj,jk)
+         END_3D
          !                                     !  ------------------------
       ENDIF
 
       !     Photosynthetically Available Radiation (PAR)
       !     --------------------------------------------
-      zqsr_corr(:,:) = parsw(:,:) * qsr(:,:) / ( 1.-fr_i(:,:) + rtrn )
+      DO_2D( 0, 0, 0, 0 )
+         zqsr_corr(ji,jj) = parsw(ji,jj) * qsr(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
+      END_2D
       !
       CALL trc_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )
       !
-      DO jk = 1, nksrp
-         etot (:,:,jk) = ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
-      ENDDO
+      DO_3D( 0, 0, 0, 0, 1, nksr )
+         etot(ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
+      END_3D
 
       ! No Diurnal cycle PAR
       IF( l_trcdm2dc ) THEN
-         zqsr_corr(:,:) = parsw(:,:) * qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn )
+         DO_2D( 0, 0, 0, 0 )
+            zqsr_corr(ji,jj) = parsw(ji,jj) * qsr_mean(ji,jj) / ( 1.-fr_i(ji,jj) + rtrn )
+         END_2D
          !
          CALL trc_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )
-         DO jk = 1, nksrp
-            etot_ndcy(:,:,jk) = ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
-         END DO
+         !
+         DO_3D( 0, 0, 0, 0, 1, nksr )
+            etot_ndcy(ji,jj,jk) = ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk)
+         END_3D
       ELSE
          etot_ndcy(:,:,:) = etot(:,:,:)
       ENDIF
 
       !     Weighted broadband attenuation coefficient
       !     ------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          ztmp = ze1(ji,jj,jk)* ekb(ji,jj,jk) + ze2(ji,jj,jk) * ekg(ji,jj,jk) + ze3(ji,jj,jk) * ekr(ji,jj,jk)
          zeps(ji,jj,jk) = ztmp / e3t(ji,jj,jk,Kmm) / (etot(ji,jj,jk) + rtrn)
       END_3D
@@ -154,26 +165,24 @@ CONTAINS
 
       !     Light at the euphotic depth
       !     ---------------------------
-      zqsr100 = 0.01 * 3. * zqsr_corr(:,:)
+      zqsr100(:,:) = 0.01 * 3. * zqsr_corr(:,:)
 
       !     Euphotic depth and level
       !     ------------------------
-      neln   (:,:) = 1
-      heup   (:,:) = gdepw(:,:,2,Kmm)
-      heup_01(:,:) = gdepw(:,:,2,Kmm)
-      !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksrp )
+      DO_2D( 0, 0, 0, 0 )
+         neln   (ji,jj) = 1
+         heup   (ji,jj) = gdepw(ji,jj,2,Kmm)
+         heup_01(ji,jj) = gdepw(ji,jj,2,Kmm)
+      END_2D
+
+      DO_3D( 0, 0, 0, 0, 2, nksr)
         IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >=  zqsr100(ji,jj) )  THEN
-           ! Euphotic level (1st T-level strictly below Euphotic layer)
-           ! NOTE: ensure compatibility with nmld_trc definition in trdmxl_trc
-           neln(ji,jj) = jk+1
-           !
-           ! Euphotic layer depth
-           heup(ji,jj) = gdepw(ji,jj,jk+1,Kmm)
+           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
+           heup(ji,jj) = gdepw(ji,jj,jk+1,Kmm)     ! Euphotic layer depth
         ENDIF
-        ! Euphotic layer depth (light level definition)
-        IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 )  THEN
-           heup_01(ji,jj) = gdepw(ji,jj,jk+1,Kmm)
+        IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.10 )  THEN
+           heup_01(ji,jj) = gdepw(ji,jj,jk+1,Kmm)  ! Euphotic layer depth (light level definition)
         ENDIF
       END_3D
       !
@@ -181,8 +190,18 @@ CONTAINS
       heup_01(:,:) = MIN( 300., heup_01(:,:) )
       !
       IF( lk_iomput ) THEN
-         CALL iom_put( "xbla" , zeps(:,:,:) * tmask(:,:,:) )
-         CALL iom_put( "Heup" , heup(:,:  ) * tmask(:,:,1) )
+        IF( iom_use( "Heup" ) ) THEN
+           ALLOCATE( zw2d(A2D(0)) )
+           zw2d(A2D(0)) = heup(A2D(0)) * tmask(A2D(0),1)
+           CALL iom_put( "Heup", zw2d )  ! Euphotic layer depth
+           DEALLOCATE( zw2d )
+        ENDIF
+        IF( iom_use( "xbla" ) ) THEN
+           ALLOCATE( zw3d(A2D(0),jpk))   ;    zw3d(A2D(0),jpk) = 0._wp
+           zw3d(A2D(0),1:jpkm1) = zeps(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+           CALL iom_put( "xbla", zw3d )  ! Euphotic layer depth
+           DEALLOCATE( zw3d )
+        ENDIF
       ENDIF
       !
       IF( ln_timing )   CALL timing_stop('trc_opt')
@@ -199,11 +218,11 @@ CONTAINS
       !!
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in)      ::   kt                ! ocean time-step
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)      ::   zqsr              ! real shortwave
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out)     ::   pe1 , pe2 , pe3   ! PAR (R-G-B)
+      REAL(wp), DIMENSION(A2D(0))    , INTENT(in)      ::   zqsr              ! real shortwave
+      REAL(wp), DIMENSION(A2D(0),jpk), INTENT(out)     ::   pe1 , pe2 , pe3   ! PAR (R-G-B)
       !
       INTEGER                       ::   ji, jj, jk        ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj)  ::   we1, we2, we3     ! PAR (R-G-B) at w-level
+      REAL(wp), DIMENSION(A2D(0))  ::   we1, we2, we3     ! PAR (R-G-B) at w-level
       !!----------------------------------------------------------------------
       pe1(:,:,:) = 0. ; pe2(:,:,:) = 0. ; pe3(:,:,:) = 0.
       !
@@ -213,7 +232,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, nksrp )
+         DO_3D( 0, 0, 0, 0, 2, nksrp )
             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) ) )
@@ -225,7 +244,7 @@ CONTAINS
          we2(:,:) = zqsr(:,:)
          we3(:,:) = zqsr(:,:)
          !
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksrp )
+         DO_3D( 0, 0, 0, 0, 1, nksrp )
             ! integrate PAR over current t-level
             pe1(ji,jj,jk) = we1(ji,jj) / (ekb(ji,jj,jk) + rtrn) * (1. - EXP( -ekb(ji,jj,jk) ))
             pe2(ji,jj,jk) = we2(ji,jj) / (ekg(ji,jj,jk) + rtrn) * (1. - EXP( -ekg(ji,jj,jk) ))
@@ -266,7 +285,9 @@ CONTAINS
       IF( ln_varpar ) THEN
          IF( kt == nittrc000 .OR. ( kt /= nittrc000 .AND. ntimes_par > 1 ) ) THEN
             CALL fld_read( kt, 1, sf_par )
-            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0
+            DO_2D( 0, 0, 0, 0 )
+               par_varsw(ji,jj) = ( sf_par(1)%fnow(ji,jj,1) ) / 3.0
+            END_2D
          ENDIF
       ENDIF
       !
@@ -348,8 +369,8 @@ CONTAINS
       !!                     ***  ROUTINE trc_opt_alloc  ***
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( ekb(jpi,jpj,jpk),ekr(jpi,jpj,jpk),  &
-                ekg(jpi,jpj,jpk),zeps(jpi,jpj,jpk),  STAT= trc_opt_alloc  ) 
+      ALLOCATE( ekb(A2D(0),jpk),ekr(A2D(0),jpk),  &
+                ekg(A2D(0),jpk),zeps(A2D(0),jpk),  STAT= trc_opt_alloc  ) 
       !
       IF( trc_opt_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trc_opt_alloc : failed to allocate arrays.' )
       !
-- 
GitLab