From 4b23ba89db7a194fc0f8e0fef5a5c59677a90b47 Mon Sep 17 00:00:00 2001
From: Renaud Person <renaud.person@locean.ipsl.fr>
Date: Wed, 20 Jul 2022 15:14:25 +0200
Subject: [PATCH] First step of PISCES/P4Z cleanup

---
 cfgs/SHARED/field_def_nemo-pisces.xml | 26 ++++++-------
 src/TOP/PISCES/P4Z/p4zagg.F90         |  4 +-
 src/TOP/PISCES/P4Z/p4zbc.F90          | 16 ++++----
 src/TOP/PISCES/P4Z/p4zche.F90         | 32 ++++++++--------
 src/TOP/PISCES/P4Z/p4zint.F90         |  6 +--
 src/TOP/PISCES/P4Z/p4zlim.F90         | 55 +++++++++++++++++++++------
 src/TOP/PISCES/P4Z/p4zlys.F90         |  2 +-
 src/TOP/PISCES/P4Z/p4zmort.F90        |  4 +-
 src/TOP/PISCES/P4Z/p4zpoc.F90         | 22 +++++------
 src/TOP/PISCES/P4Z/p5zmort.F90        |  6 +--
 10 files changed, 103 insertions(+), 70 deletions(-)

diff --git a/cfgs/SHARED/field_def_nemo-pisces.xml b/cfgs/SHARED/field_def_nemo-pisces.xml
index e15f6d55f..41e65d48f 100644
--- a/cfgs/SHARED/field_def_nemo-pisces.xml
+++ b/cfgs/SHARED/field_def_nemo-pisces.xml
@@ -188,7 +188,7 @@
     <field id="PFeN"        long_name="Primary production of nano iron"         unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="PFeP"        long_name="Primary production of pico iron"         unit="molC/m3/s"  grid_ref="grid_T_3D" />
     <field id="PFeD"        long_name="Primary production of diatoms iron"      unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
-    <field id="xfracal"     long_name="Calcifying fraction"                     unit="1"          grid_ref="grid_T_3D" />
+    <field id="xfracal"     long_name="Calcifying fraction"                     unit="1"          grid_ref="grid_T_3D_inner" />
     <field id="PCAL"        long_name="Calcite production"                      unit="mol/m3/s"   grid_ref="grid_T_3D" />
     <field id="GRAZ1"       long_name="Grazing by microzooplankton"             unit="mol/m3/s"   grid_ref="grid_T_3D" />
     <field id="GRAZ2"       long_name="Grazing by mesozooplankton"              unit="mol/m3/s"   grid_ref="grid_T_3D" />
@@ -204,21 +204,21 @@
     <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="LNnut"       long_name="Nutrient limitation term in Nanophyto"   unit=""           grid_ref="grid_T_3D" />
-    <field id="LPnut"       long_name="Nutrient limitation term in Picophyto"   unit="-"          grid_ref="grid_T_3D" />
-    <field id="LDnut"       long_name="Nutrient limitation term in Diatoms"     unit=""           grid_ref="grid_T_3D" />
-    <field id="LNFe"        long_name="Iron limitation term in Nanophyto"       unit=""           grid_ref="grid_T_3D" />
-    <field id="LPFe"        long_name="Iron limitation term in Picophyto"       unit="-"          grid_ref="grid_T_3D" />
-    <field id="LDFe"        long_name="Iron limitation term in Diatoms"         unit=""           grid_ref="grid_T_3D" />
+    <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" />
+    <field id="LNFe"        long_name="Iron limitation term in Nanophyto"       unit=""           grid_ref="grid_T_3D_inner" />
+    <field id="LPFe"        long_name="Iron limitation term in Picophyto"       unit="-"          grid_ref="grid_T_3D_inner" />
+    <field id="LDFe"        long_name="Iron limitation term in Diatoms"         unit=""           grid_ref="grid_T_3D_inner" />
     <field id="LNlight"     long_name="Light limitation term in Nanophyto"      unit=""           grid_ref="grid_T_3D_inner" />
     <field id="LPlight"     long_name="Light limitation term in Picophyto"      unit="-"          grid_ref="grid_T_3D" />
     <field id="LDlight"     long_name="Light limitation term in Diatoms"        unit=""           grid_ref="grid_T_3D_inner" />
-    <field id="SIZEN"       long_name="Mean relative size of nanophyto."        unit="-"          grid_ref="grid_T_3D" />
-    <field id="SIZEP"       long_name="Mean relative size of picophyto."        unit="-"          grid_ref="grid_T_3D" />
-    <field id="SIZED"       long_name="Mean relative size of diatoms"           unit="-"          grid_ref="grid_T_3D" />
-    <field id="RASSD"       long_name="Size of the protein machinery (Diat.)"   unit="-"          grid_ref="grid_T_3D" />
-    <field id="RASSN"       long_name="Size of the protein machinery (Nano.)"   unit="-"          grid_ref="grid_T_3D" />
-    <field id="RASSP"       long_name="Size of the protein machinery (Pico.)"   unit="-"          grid_ref="grid_T_3D" />
+    <field id="SIZEN"       long_name="Mean relative size of nanophyto."        unit="-"          grid_ref="grid_T_3D_inner" />
+    <field id="SIZEP"       long_name="Mean relative size of picophyto."        unit="-"          grid_ref="grid_T_3D_inner" />
+    <field id="SIZED"       long_name="Mean relative size of diatoms"           unit="-"          grid_ref="grid_T_3D_inner" />
+    <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" />
diff --git a/src/TOP/PISCES/P4Z/p4zagg.F90 b/src/TOP/PISCES/P4Z/p4zagg.F90
index 6eb5208f3..9c16b1647 100644
--- a/src/TOP/PISCES/P4Z/p4zagg.F90
+++ b/src/TOP/PISCES/P4Z/p4zagg.F90
@@ -70,7 +70,7 @@ CONTAINS
       ! PISCES part
       IF( ln_p4z ) THEN
          !
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             !
             zfact = xstep * xdiss(ji,jj,jk)
             ! Part I : Coagulation dependent on turbulence
@@ -117,7 +117,7 @@ CONTAINS
       ELSE    ! ln_p5z
         ! PISCES-QUOTA part
         !
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             !
             zfact = xstep * xdiss(ji,jj,jk)
             !  Part I : Coagulation dependent on turbulence
diff --git a/src/TOP/PISCES/P4Z/p4zbc.F90 b/src/TOP/PISCES/P4Z/p4zbc.F90
index f649f91f7..f0bd7da61 100644
--- a/src/TOP/PISCES/P4Z/p4zbc.F90
+++ b/src/TOP/PISCES/P4Z/p4zbc.F90
@@ -99,7 +99,7 @@ CONTAINS
 
          ! Atmospheric input of Iron dissolves in the water column
          IF ( ln_trc_sbc(jpfer) ) THEN
-            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
                !
                tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zdustdep * mfrac / mMass_Fe 
@@ -116,7 +116,7 @@ CONTAINS
 
          ! Atmospheric input of PO4 dissolves in the water column
          IF ( ln_trc_sbc(jppo4) ) THEN
-            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
                !
                tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zdustdep * 1.e-3 / mMass_P
@@ -125,7 +125,7 @@ CONTAINS
 
          ! Atmospheric input of Si dissolves in the water column
          IF ( ln_trc_sbc(jpsil) ) THEN
-            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
                zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) )
                !
                tr(ji,jj,jk,jpsil,Krhs) = tr(ji,jj,jk,jpsil,Krhs) + zdustdep * 0.269 / mMass_Si
@@ -144,7 +144,7 @@ CONTAINS
       ! ----------------------------------------------------------
       IF( ll_river ) THEN
           jl = n_trc_indcbc(jpno3)
-          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+          DO_2D( 0, 0, 0, 0 )
              DO jk = 1, nk_rnf(ji,jj)
                 zcoef = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_cbc_time ) * tmask(ji,jj,1)
                 zrivdin = rf_trcfac(jl) * sf_trccbc(jl)%fnow(ji,jj,1) * zcoef
@@ -158,14 +158,14 @@ CONTAINS
       IF( ll_ndepo ) THEN
          IF( ln_trc_sbc(jpno3) ) THEN
             jl = n_trc_indsbc(jpno3)
-            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            DO_2D( 0, 0, 0, 0 )
                zndep = rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) / e3t(ji,jj,1,Kmm) / rn_sbc_time
                tr(ji,jj,1,jptal,Krhs) = tr(ji,jj,1,jptal,Krhs) - rno3 * zndep * rfact
             END_2D
          ENDIF
          IF( ln_trc_sbc(jpnh4) ) THEN
             jl = n_trc_indsbc(jpnh4)
-            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            DO_2D( 0, 0, 0, 0 )
                zndep = rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) / e3t(ji,jj,1,Kmm) / rn_sbc_time
                tr(ji,jj,1,jptal,Krhs) = tr(ji,jj,1,jptal,Krhs) + rno3 * zndep * rfact
             END_2D
@@ -183,7 +183,7 @@ CONTAINS
          ! Simple parameterization assuming a fixed constant concentration in
          ! sea-ice (icefeinput)
          ! ------------------------------------------------------------------         
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
             zdep     = rfact / e3t(ji,jj,1,Kmm)
             zwflux   = fmmflx(ji,jj) / 1000._wp
             zironice =  MAX( -0.99 * tr(ji,jj,1,jpfer,Kbb), -zwflux * icefeinput * zdep )
@@ -350,7 +350,7 @@ CONTAINS
          !
          CALL lbc_lnk( 'p4zbc', zcmask , 'T', 1.0_wp )      ! lateral boundary conditions on cmask   (sign unchanged)
          !
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             zexpide   = MIN( 8.,( gdept(ji,jj,jk,Kmm) / 500. )**(-1.5) )
             zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
             zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
diff --git a/src/TOP/PISCES/P4Z/p4zche.F90 b/src/TOP/PISCES/P4Z/p4zche.F90
index a71f9350e..cacce4c2a 100644
--- a/src/TOP/PISCES/P4Z/p4zche.F90
+++ b/src/TOP/PISCES/P4Z/p4zche.F90
@@ -179,7 +179,7 @@ CONTAINS
       ! potential temperature to in situ temperature. The errors is less than 
       ! 0.04°C relative to an exact computation
       ! ---------------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+      DO_3D( 0, 0, 0, 0, 1, jpk )
          zpres = gdept(ji,jj,jk,Kmm) / 1000.
          za1 = 0.04 * ( 1.0 + 0.185 * ts(ji,jj,jk,jp_tem,Kmm) + 0.035 * (salinprac(ji,jj,jk) - 35.0) )
          za2 = 0.0075 * ( 1.0 - ts(ji,jj,jk,jp_tem,Kmm) / 30.0 )
@@ -188,7 +188,7 @@ CONTAINS
       !
       ! CHEMICAL CONSTANTS - SURFACE LAYER
       ! ----------------------------------
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          !                             ! SET ABSOLUTE TEMPERATURE
          ztkel = tempis(ji,jj,1) + 273.15
          zt    = ztkel * 0.01
@@ -216,7 +216,7 @@ CONTAINS
 
       ! OXYGEN SOLUBILITY - DEEP OCEAN
       ! -------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+      DO_3D( 0, 0, 0, 0, 1, jpk )
          ztkel = tempis(ji,jj,jk) + 273.15
          zsal  = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35.
          zsal2 = zsal * zsal
@@ -235,7 +235,7 @@ CONTAINS
 
       ! CHEMICAL CONSTANTS - DEEP OCEAN
       ! -------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+      DO_3D( 0, 0, 0, 0, 1, jpk )
           ! SET PRESSION ACCORDING TO SAUNDER (1980)
           zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) )
           zc1 = 5.92E-3 + zplat**2 * 5.25E-3
@@ -452,7 +452,7 @@ CONTAINS
       !!                    and the 2nd order approximation does not have 
       !!                    a solution
       !!---------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT)  ::  p_hini
+      REAL(wp), DIMENSION(A2D(0),jpk), INTENT(OUT)  ::  p_hini
       INTEGER,                          INTENT(in)   ::  Kbb      ! time level indices
       INTEGER  ::   ji, jj, jk
       REAL(wp)  ::  zca1, zba1
@@ -463,7 +463,7 @@ CONTAINS
 
       IF( ln_timing )  CALL timing_start('ahini_for_at')
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+      DO_3D( 0, 0, 0, 0, 1, jpk )
       zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
       p_alkcb  = tr(ji,jj,jk,jptal,Kbb) * zrhd
       p_dictot = tr(ji,jj,jk,jpdic,Kbb) * zrhd
@@ -512,13 +512,13 @@ CONTAINS
    ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+])
 
    ! Argument variables
-   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf
-   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup
+   REAL(wp), DIMENSION(A2D(0),jpk), INTENT(OUT) :: p_alknw_inf
+   REAL(wp), DIMENSION(A2D(0),jpk), INTENT(OUT) :: p_alknw_sup
    INTEGER,                          INTENT(in)  ::  Kbb      ! time level indices
    INTEGER  ::   ji, jj, jk
    REAL(wp)  ::  zrhd
 
-    DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+    DO_3D( 0, 0, 0, 0, 1, jpk )
       zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
       p_alknw_inf(ji,jj,jk) =  -tr(ji,jj,jk,jppo4,Kbb) * zrhd - sulfat(ji,jj,jk) &
       &              - fluorid(ji,jj,jk)
@@ -536,8 +536,8 @@ CONTAINS
 
    ! Argument variables
    !--------------------
-   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN)   :: p_hini
-   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT)  :: zhi
+   REAL(wp), DIMENSION(A2D(0),jpk), INTENT(IN)   :: p_hini
+   REAL(wp), DIMENSION(A2D(0),jpk), INTENT(OUT)  :: zhi
    INTEGER,                          INTENT(in)   :: Kbb  ! time level indices
 
    ! Local variables
@@ -557,17 +557,17 @@ CONTAINS
    REAL(wp)  ::  zrhd, p_alktot, zdic, zbot, zpt, zst, zft, zsit
    LOGICAL   ::  l_exitnow
    REAL(wp), PARAMETER :: pz_exp_threshold = 1.0
-   REAL(wp), DIMENSION(jpi,jpj,jpk) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin
+   REAL(wp), DIMENSION(A2D(0),jpk) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin
 
    IF( ln_timing )  CALL timing_start('solve_at_general')
 
    CALL anw_infsup( zalknw_inf, zalknw_sup, Kbb )
 
-   rmask(:,:,:) = tmask(:,:,:)
+   rmask(A2D(0),1:jpk) = tmask(A2D(0),1:jpk)
    zhi(:,:,:)   = 0.
 
    ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
-   DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+   DO_3D( 0, 0, 0, 0, 1, jpk )
       IF (rmask(ji,jj,jk) == 1.) THEN
          zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
          p_alktot = tr(ji,jj,jk,jptal,Kbb) * zrhd
@@ -597,7 +597,7 @@ CONTAINS
    zeqn_absmin(:,:,:) = HUGE(1._wp)
 
    DO jn = 1, jp_maxniter_atgen 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+      DO_3D( 0, 0, 0, 0, 1, jpk )
       IF (rmask(ji,jj,jk) == 1.) THEN
          zrhd = 1._wp / ( rhd(ji,jj,jk) + 1. )
          p_alktot = tr(ji,jj,jk,jptal,Kbb) * zrhd
@@ -799,7 +799,7 @@ CONTAINS
 
       ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), STAT=ierr(1) )
 
-      ALLOCATE( akb3(jpi,jpj,jpk)     , tempis(jpi, jpj, jpk),       &
+      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)   ,       &
diff --git a/src/TOP/PISCES/P4Z/p4zint.F90 b/src/TOP/PISCES/P4Z/p4zint.F90
index 9b95c2e49..3b19ae86d 100644
--- a/src/TOP/PISCES/P4Z/p4zint.F90
+++ b/src/TOP/PISCES/P4Z/p4zint.F90
@@ -53,7 +53,7 @@ CONTAINS
       ! Computation of the silicon dependant half saturation  constant for silica uptake
       ! This is based on an old study by Pondaven et al. (1998)
       ! --------------------------------------------------------------------------------
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          zvar = tr(ji,jj,1,jpsil,Kbb) * tr(ji,jj,1,jpsil,Kbb)
          xksimax(ji,jj) = MAX( xksimax(ji,jj), ( 1.+ 7.* zvar / ( xksilim * xksilim + zvar ) ) * 1e-6 )
       END_2D
@@ -74,13 +74,13 @@ CONTAINS
 
       ! day length in hours
       strn(:,:) = 0.
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
          zargu = MAX( -1., MIN(  1., zargu ) )
          strn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
       END_2D
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
         ! denitrification factor computed from O2 levels
          ! This factor diagnoses below which level of O2 denitrification
          ! is active
diff --git a/src/TOP/PISCES/P4Z/p4zlim.F90 b/src/TOP/PISCES/P4Z/p4zlim.F90
index 4e227808a..93a988c71 100644
--- a/src/TOP/PISCES/P4Z/p4zlim.F90
+++ b/src/TOP/PISCES/P4Z/p4zlim.F90
@@ -98,13 +98,14 @@ CONTAINS
       REAL(wp) ::   zconc1d, zconc1dnh4, zconc0n, zconc0nnh4   
       REAL(wp) ::   fananof, fadiatf, znutlim, zfalim
       REAL(wp) ::   znutlimtot, zlimno3, zlimnh4, zbiron
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_lim')
       !
       sizena(:,:,:) = 1.0  ;  sizeda(:,:,:) = 1.0
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          
          ! Computation of a variable Ks for iron on diatoms taking into account
          ! that increasing biomass is made of generally bigger cells
@@ -219,7 +220,7 @@ CONTAINS
       ! Size estimation of phytoplankton based on total biomass
       ! Assumes that larger biomass implies addition of larger cells
       ! ------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          zcoef = tr(ji,jj,jk,jpphy,Kbb) - MIN(xsizephy, tr(ji,jj,jk,jpphy,Kbb) )
          sizena(ji,jj,jk) = 1. + ( xsizern -1.0 ) * zcoef / ( xsizephy + zcoef )
          zcoef = tr(ji,jj,jk,jpdia,Kbb) - MIN(xsizedia, tr(ji,jj,jk,jpdia,Kbb) )
@@ -231,7 +232,7 @@ CONTAINS
       ! This is a purely adhoc formulation described in Aumont et al. (2015)
       ! This fraction depends on nutrient limitation, light, temperature
       ! --------------------------------------------------------------------
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          zlim1  = xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) 
          zlim2  = tr(ji,jj,jk,jppo4,Kbb) / ( tr(ji,jj,jk,jppo4,Kbb) + concnnh4 )
          zlim3  = tr(ji,jj,jk,jpfer,Kbb) / ( tr(ji,jj,jk,jpfer,Kbb) +  6.E-11   )
@@ -250,7 +251,7 @@ CONTAINS
          xfracal(ji,jj,jk) = MAX( 0.02, xfracal(ji,jj,jk) )
       END_3D
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          ! denitrification factor computed from O2 levels
          nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - tr(ji,jj,jk,jpoxy,Kbb) )    &
             &                                / ( oxymin + tr(ji,jj,jk,jpoxy,Kbb) )  )
@@ -263,13 +264,45 @@ CONTAINS
       END_3D
       !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics
-        CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! euphotic layer deptht
-        CALL iom_put( "LNnut"  , xlimphy(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
-        CALL iom_put( "LDnut"  , xlimdia(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
-        CALL iom_put( "LNFe"   , xlimnfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "LDFe"   , xlimdfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "SIZEN"  , sizen  (:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-        CALL iom_put( "SIZED"  , sized  (:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        !
+        ALLOCATE( zw3d(A2D(0),jpk) )  ;  zw3d(A2D(0),jpk) = 0._wp
+        !
+        IF( iom_use ( "xfracal" ) ) THEN   ! euphotic layer deptht
+          zw3d(A2D(0),1:jpkm1) = xfracal(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "xfracal",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "LNnut" ) ) THEN   ! Nutrient limitation term
+          zw3d(A2D(0),1:jpkm1) = xlimphy(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "LNnut",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "LDnut" ) ) THEN   ! Nutrient limitation term
+          zw3d(A2D(0),1:jpkm1) = xlimdia(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "LDnut",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "LNFe" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = xlimnfe(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "LNFe",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "LDFe" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = xlimdfe(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "LDFe",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "SIZEN" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = sizen(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "SIZEN",  zw3d)
+        ENDIF
+        !
+        IF( iom_use ( "SIZED" ) ) THEN   ! Iron limitation term
+          zw3d(A2D(0),1:jpkm1) = sized(A2D(0),1:jpkm1) * tmask(A2D(0),1:jpkm1)
+          CALL iom_put( "SIZED",  zw3d)
+        ENDIF
+        !
+        DEALLOCATE( zw3d )
       ENDIF
       !
       IF( ln_timing )   CALL timing_stop('p4z_lim')
diff --git a/src/TOP/PISCES/P4Z/p4zlys.F90 b/src/TOP/PISCES/P4Z/p4zlys.F90
index 24d573977..9450e39f3 100644
--- a/src/TOP/PISCES/P4Z/p4zlys.F90
+++ b/src/TOP/PISCES/P4Z/p4zlys.F90
@@ -129,7 +129,7 @@ CONTAINS
       !
 
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-             ALLOCATE( zw3d(A2D(0),jpk) )    ;    zw3d(A2D(0),jpk) = 0._wp
+          ALLOCATE( zw3d(A2D(0),jpk) )    ;    zw3d(A2D(0),jpk) = 0._wp
           IF( iom_use( "PH" ) ) THEN
              DO_3D( 0, 0, 0, 0, 1, jpkm1)
                 zw3d(ji,jj,jk) = -1. * LOG10( MAX( hi(ji,jj,jk), rtrn ) ) * tmask(ji,jj,jk) 
diff --git a/src/TOP/PISCES/P4Z/p4zmort.F90 b/src/TOP/PISCES/P4Z/p4zmort.F90
index 5a32dd99e..3c3903e74 100644
--- a/src/TOP/PISCES/P4Z/p4zmort.F90
+++ b/src/TOP/PISCES/P4Z/p4zmort.F90
@@ -74,7 +74,7 @@ CONTAINS
       IF( ln_timing )   CALL timing_start('p4z_mort_nano')
       !
       prodcal(:,:,:) = 0._wp   ! calcite production variable set to zero
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          zcompaph = MAX( ( tr(ji,jj,jk,jpphy,Kbb) - 1e-9 ), 0.e0 )
 
          ! Quadratic mortality of nano due to aggregation during
@@ -152,7 +152,7 @@ CONTAINS
       ! This is due to the production of EPS by stressed cells
       ! -------------------------------------------------------------
 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
 
          zcompadi = MAX( ( tr(ji,jj,jk,jpdia,Kbb) - 1e-9), 0. )
 
diff --git a/src/TOP/PISCES/P4Z/p4zpoc.F90 b/src/TOP/PISCES/P4Z/p4zpoc.F90
index 54b8f202e..d729eff71 100644
--- a/src/TOP/PISCES/P4Z/p4zpoc.F90
+++ b/src/TOP/PISCES/P4Z/p4zpoc.F90
@@ -74,9 +74,9 @@ CONTAINS
       REAL(wp) ::   zofer2, zofer3
       REAL(wp) ::   zrfact2
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj  )   :: totprod, totthick, totcons 
-      REAL(wp), DIMENSION(jpi,jpj,jpk)   :: zremipoc, zremigoc, zorem3, ztremint, zfolimi
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jcpoc) :: alphag
+      REAL(wp), DIMENSION(A2D(0)  )   :: totprod, totthick, totcons 
+      REAL(wp), DIMENSION(A2D(0),jpk)   :: zremipoc, zremigoc, zorem3, ztremint, zfolimi
+      REAL(wp), DIMENSION(A2D(0),jpk,jcpoc) :: alphag
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )  CALL timing_start('p4z_poc')
@@ -118,7 +118,7 @@ CONTAINS
      ! a standard parameterisation with a constant lability
      ! -----------------------------------------------------------------------
      ztremint(:,:,:) = zremigoc(:,:,:)
-     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.) THEN
           zdep = hmld(ji,jj)
           !
@@ -204,7 +204,7 @@ CONTAINS
 
       IF( ln_p4z ) THEN
          ! The standard PISCES part
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
             ! POC degradation by bacterial activity. It is a function 
             ! of the mean lability and of temperature. This also includes
             ! shrinking of particles due to the bacterial activity
@@ -226,7 +226,7 @@ CONTAINS
             zfolimi(ji,jj,jk)   = zofer2
          END_3D
       ELSE
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1)
             ! POC degradation by bacterial activity. It is a function 
             ! of the mean lability and of temperature. This also includes
             ! shrinking of particles due to the bacterial activity
@@ -274,7 +274,7 @@ CONTAINS
 
      ! intregrated production and consumption of POC in the mixed layer
      ! ----------------------------------------------------------------
-     DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+     DO_3D( 0, 0, 0, 0, 1, jpkm1)
         zdep = hmld(ji,jj)
         IF (tmask(ji,jj,jk) == 1. .AND. gdept(ji,jj,jk,Kmm) <= zdep ) THEN
           totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * e3t(ji,jj,jk,Kmm) * rday/ rfact2
@@ -290,7 +290,7 @@ CONTAINS
      ! mixing.
      ! ---------------------------------------------------------------------
      ztremint(:,:,:) = zremipoc(:,:,:)
-     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.) THEN
           zdep = hmld(ji,jj)
           alphat = 0.0
@@ -323,7 +323,7 @@ CONTAINS
      ! that since we need the lability spectrum of GOC, GOC spectrum 
      ! should be determined before.
      ! -----------------------------------------------------------------------
-     DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1)
+     DO_3D( 0, 0, 0, 0, 2, jpkm1)
         IF (tmask(ji,jj,jk) == 1.) THEN
           zdep = hmld(ji,jj)
           IF( gdept(ji,jj,jk,Kmm) > zdep ) THEN
@@ -397,7 +397,7 @@ CONTAINS
      ENDIF
 
      IF( ln_p4z ) THEN
-         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.) THEN
               ! POC disaggregation by turbulence and bacterial activity.It is a function
               ! of the mean lability and of temperature  
@@ -416,7 +416,7 @@ CONTAINS
             ENDIF
          END_3D
      ELSE
-       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+       DO_3D( 0, 0, 0, 0, 1, jpkm1)
           ! POC disaggregation by turbulence and bacterial activity.It is a function
           ! of the mean lability and of temperature  
           !--------------------------------------------------------
diff --git a/src/TOP/PISCES/P4Z/p5zmort.F90 b/src/TOP/PISCES/P4Z/p5zmort.F90
index 396be16f7..4fb33dac3 100644
--- a/src/TOP/PISCES/P4Z/p5zmort.F90
+++ b/src/TOP/PISCES/P4Z/p5zmort.F90
@@ -80,7 +80,7 @@ CONTAINS
       IF( ln_timing )   CALL timing_start('p5z_mort_nano')
       !
       prodcal(:,:,:) = 0.  !: calcite production variable set to zero
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          zcompaph = MAX( ( tr(ji,jj,jk,jpphy,Kbb) - 1e-9 ), 0.e0 )
 
          ! Quadratic mortality of nano due to aggregation during
@@ -151,7 +151,7 @@ CONTAINS
       !
       IF( ln_timing )   CALL timing_start('p5z_mort_pico')
       !
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          zcompaph = MAX( ( tr(ji,jj,jk,jppic,Kbb) - 1e-9 ), 0.e0 )
 
          ! Quadratic mortality of pico due to aggregation during
@@ -215,7 +215,7 @@ CONTAINS
       IF( ln_timing )   CALL timing_start('p5z_mort_diat')
       !
 
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
 
          zcompadi = MAX( ( tr(ji,jj,jk,jpdia,Kbb) - 1E-9), 0. )
 
-- 
GitLab