From f09c3227c396d4809301dec4f6ab7db5c2229436 Mon Sep 17 00:00:00 2001
From: cetlod <Christian.Ethe@ipsl.fr>
Date: Thu, 14 Jul 2022 12:39:56 +0200
Subject: [PATCH] TOP/PISCES cleaning : continuation

---
 cfgs/SHARED/field_def_nemo-pisces.xml | 20 ++++-----
 src/OCE/IOM/iom.F90                   | 59 +++++++++++----------------
 src/OCE/IOM/iom_nf90.F90              | 39 ++++++++----------
 src/TOP/PISCES/P4Z/p4zlys.F90         | 14 +++----
 src/TOP/PISCES/P4Z/p4zprod.F90        | 11 +++--
 src/TOP/PISCES/sms_pisces.F90         | 44 ++++++++++----------
 src/TOP/trc.F90                       | 12 +++---
 src/TOP/trcwri.F90                    |  4 +-
 8 files changed, 93 insertions(+), 110 deletions(-)

diff --git a/cfgs/SHARED/field_def_nemo-pisces.xml b/cfgs/SHARED/field_def_nemo-pisces.xml
index 6b63d355e..e15f6d55f 100644
--- a/cfgs/SHARED/field_def_nemo-pisces.xml
+++ b/cfgs/SHARED/field_def_nemo-pisces.xml
@@ -176,20 +176,20 @@
     <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="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="PPPHYP"      long_name="Primary production of picophyto"         unit="molC/m3/s"  grid_ref="grid_T_3D" />
     <field id="PPPHYD"      long_name="Primary production of diatoms"           unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="PPNEWN"      long_name="New Primary production of nanophyto"     unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="PPNEWP"      long_name="New Primary production of picophyto"     unit="molC/m3/s"  grid_ref="grid_T_3D" />
     <field id="PPNEWD"      long_name="New Primary production of diatoms"       unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="PBSi"        long_name="Primary production of Si diatoms"        unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="PFeN"        long_name="Primary production of nano iron"         unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
-    <field id="PFeP"        long_name="Primary production of pico iron"         unit="molC/m3/s"  grid_ref="grid_T_3D_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="PCAL"        long_name="Calcite production"                      unit="mol/m3/s"   grid_ref="grid_T_3D" />
-    <field id="DCAL"        long_name="Calcite dissolution"                     unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
     <field id="GRAZ1"       long_name="Grazing by microzooplankton"             unit="mol/m3/s"   grid_ref="grid_T_3D" />
     <field id="GRAZ2"       long_name="Grazing by mesozooplankton"              unit="mol/m3/s"   grid_ref="grid_T_3D" />
     <field id="REMIN"       long_name="Oxic remineralization of OM"             unit="mol/m3/s"   grid_ref="grid_T_3D" />
@@ -199,19 +199,19 @@
     <field id="Nfix"        long_name="Nitrogen fixation"                       unit="mol/m3/s"   grid_ref="grid_T_3D" />
     <field id="Mumax"       long_name="Maximum growth rate"                     unit="s-1"        grid_ref="grid_T_3D_inner" />
     <field id="MuN"         long_name="Realized growth rate for nanophyto"      unit="s-1"        grid_ref="grid_T_3D_inner" />
-    <field id="MuP"         long_name="Realized growth rate for picophyto"      unit="s-1"        grid_ref="grid_T_3D_inner" />
+    <field id="MuP"         long_name="Realized growth rate for picophyto"      unit="s-1"        grid_ref="grid_T_3D" />
     <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="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="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="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_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" />
@@ -257,10 +257,10 @@
     <field id="FEBACT"      long_name="Bacterial uptake of Fe"                  unit="molFe/m3/s"  grid_ref="grid_T_3D" />
     <field id="FEPREC"      long_name="Precipitation of Fe"                     unit="molFe/m3/s"  grid_ref="grid_T_3D" />
     <field id="LPRODR"      long_name="OM remineralisation ligand production rate" unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
-    <field id="LPRODP"      long_name="phytoplankton ligand production rate"    unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
+    <field id="LPRODP"      long_name="phytoplankton ligand production rate"    unit="nmol-L/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="LIGREM"      long_name="Remineralisation loss of ligands"        unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
     <field id="LIGPR"       long_name="Photochemical loss of ligands"           unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
-    <field id="LDETP"       long_name="Ligand destruction during phytoplankton uptake" unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
+    <field id="LDETP"       long_name="Ligand destruction during phytoplankton uptake" unit="nmol-L/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="LPRODZ2"     long_name="mesozooplankton ligand production rate"  unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
     <field id="LPRODZ"      long_name="microzooplankton ligand production rate" unit="nmol-L/m3/s"  grid_ref="grid_T_3D" />
     <field id="FEZOO"       long_name="microzooplankton iron recycling rate"    unit="nmol-FeL/m3/s"  grid_ref="grid_T_3D" />
diff --git a/src/OCE/IOM/iom.F90 b/src/OCE/IOM/iom.F90
index 9a1f1c859..d66b13d19 100644
--- a/src/OCE/IOM/iom.F90
+++ b/src/OCE/IOM/iom.F90
@@ -1202,7 +1202,6 @@ CONTAINS
       CHARACTER(LEN=1)               ::   cl_type     ! local value of cd_type
       LOGICAL                        ::   ll_only3rd  ! T => if kstart, kcount present then *only* use values for 3rd spatial dimension.
       INTEGER                        ::   inlev       ! number of levels for 3D data
-      INTEGER                        ::   ihls        ! local value of the halo size
       REAL(dp)                       ::   gma, gmi
       !---------------------------------------------------------------------
       CHARACTER(LEN=lc)                               ::   context
@@ -1328,29 +1327,14 @@ CONTAINS
                IF( irankpv == 1 )        ishape(1:1) = SHAPE(pv_r1d)
                IF( irankpv == 2 )        ishape(1:2) = SHAPE(pv_r2d)
                IF( irankpv == 3 )        ishape(1:3) = SHAPE(pv_r3d)
-               ix1 = 1      ;   ix2 = icnt(1)   ;   iy1 = 1   ;   iy2 = icnt(2)         ! index of the array to be read
                ctmp1 = 'd'
             ELSE
-               IF( irankpv == 2 )   ishape(1:2) = SHAPE(pv_r2d)
-               IF( irankpv == 3 )   ishape(1:3) = SHAPE(pv_r3d)
-               IF(     ishape(1) == Ni_0   .AND. ishape(2) == Nj_0   ) THEN           ! array with 0 halo
-                  ihls = 0
-                  ix1 = 1      ;   ix2 = Ni_0     ;   iy1 = 1      ;   iy2 = Nj_0     ! index of the array to be read
-                  ctmp1 = 'd(:,:'
-               ELSEIF( ishape(1) == jpi    .AND. ishape(2) == jpj    ) THEN           ! array with nn_hls halos
-                  ihls = nn_hls
-                  ix1 = Nis0   ;   ix2 = Nie0     ;   iy1 = Njs0   ;   iy2 = Nje0     ! index of the array to be read
-                  ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0'
-               ELSEIF( ishape(1) == Ni_0+1 .AND. ishape(2) == Nj_0+1 ) THEN           ! nn_hls = 2 and array with 1 halo
-                  ihls = 1
-                  ix1 = 2      ;   ix2 = Ni_0+1   ;   iy1 = 2      ;   iy2 = Nj_0+1   ! index of the array to be read
-                  ctmp1 = 'd(2:Ni_0+1,2:Ni_0+1'
-               ELSE
-                  CALL ctl_stop( 'iom_get_123d: should have been an impossible case...' )
+               IF( irankpv == 2 ) THEN
+                  ishape(1:2) = SHAPE(pv_r2d(Nis0:Nie0,Njs0:Nje0  ))   ;   ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0)'
+               ENDIF
+               IF( irankpv == 3 ) THEN
+                  ishape(1:3) = SHAPE(pv_r3d(Nis0:Nie0,Njs0:Nje0,:))   ;   ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0,:)'
                ENDIF
-               ishape(1:2) = (/ Ni_0, Nj_0 /)   ! update and force ishape to match the inner domain
-               IF( irankpv == 3 )   ctmp1 = TRIM(ctmp1)//',:'
-               ctmp1 = TRIM(ctmp1)//')'
             ENDIF
             DO jl = 1, irankpv
                WRITE( ctmp2, FMT="(', ', i1,'): ', i5,' /= icnt(', i1,'):', i5)" ) jl, ishape(jl), jl, icnt(jl)
@@ -1363,22 +1347,25 @@ CONTAINS
          !-
          IF( idvar > 0 .AND. istop == nstop ) THEN   ! no additional errors until this point...
             !
+            ! find the right index of the array to be read
+            IF( idom /= jpdom_unknown ) THEN   ;   ix1 = Nis0   ;   ix2 = Nie0      ;   iy1 = Njs0   ;   iy2 = Nje0
+            ELSE                               ;   ix1 = 1      ;   ix2 = icnt(1)   ;   iy1 = 1      ;   iy2 = icnt(2)
+            ENDIF
+
             CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, pv_r3d )
 
             IF( istop == nstop ) THEN   ! no additional errors until this point...
                IF(lwp) WRITE(numout,"(10x,' read ',a,' (rec: ',i6,') in ',a,' ok')") TRIM(cdvar), itime, TRIM(iom_file(kiomid)%name)
 
-               !--- halos and NP folding (NP folding to be done even if no halos)
-               IF( idom /= jpdom_unknown .AND. cl_type /= 'Z' .AND. ( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) ) THEN
-                  cl_type = 'T'
-                  IF( PRESENT(cd_type) )   cl_type = cd_type
-                  zsgn = 1._wp
-                  IF( PRESENT(psgn   ) )   zsgn    = psgn
-                  IF(     PRESENT(pv_r2d) .AND. llok ) THEN
-                     CALL lbc_lnk( 'iom', pv_r2d, cl_type, zsgn, kfillmode = kfill, khls = ihls )
-                  ELSEIF( PRESENT(pv_r3d) .AND. llok ) THEN
-                     CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill, khls = ihls )
-                  ENDIF
+               cl_type = 'T'
+               IF( PRESENT(cd_type) )   cl_type = cd_type
+               zsgn = 1._wp
+               IF( PRESENT(psgn   ) )   zsgn    = psgn
+               !--- overlap areas and extra hallows (mpp)
+               IF(     PRESENT(pv_r2d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+                  CALL lbc_lnk( 'iom', pv_r2d, cl_type, zsgn, kfillmode = kfill )
+               ELSEIF( PRESENT(pv_r3d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+                  CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill )
                ENDIF
                !
             ELSE
@@ -2349,14 +2336,14 @@ CONTAINS
             idb(jn) = -nn_hls                         ! Tile data offset (halo size)
          END DO
 
+         ! Tile_[ij]begin are defined with respect to the processor data domain, so data_[ij]begin is added
          CALL iom_set_domain_attr("grid_"//cdgrd, ntiles=nijtile,                                     &
-            & tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
+            & tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
             & tile_ni=ini(:), tile_nj=inj(:),                                                         &
             & tile_data_ibegin=idb(:), tile_data_jbegin=idb(:),                                       &
             & tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
-         idb(:) = 0
          CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", ntiles=nijtile,                           &
-            & tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
+            & tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
             & tile_ni=ini(:), tile_nj=inj(:),                                                         &
             & tile_data_ibegin=idb(:), tile_data_jbegin=idb(:),                                       &
             & tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
@@ -2466,7 +2453,7 @@ CONTAINS
 !      CALL dom_ngb( -168.53_wp, 65.03_wp, ix, iy, 'T' ) !  i-line that passes through Bering Strait: Reference latitude (used in plots)
       CALL dom_ngb( 180.0_wp, 90.0_wp, ix, iy, 'T' ) !  i-line that passes near the North Pole : Reference latitude (used in plots)
       CALL iom_set_domain_attr("gznl", ni_glo=Ni0glo, nj_glo=Nj0glo, ibegin=mig0(Nis0)-1, jbegin=mjg0(Njs0)-1, ni=Ni_0, nj=Nj_0)
-      CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin=0, data_ni=Ni_0, data_jbegin=0, data_nj=Nj_0)
+      CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj)
       CALL iom_set_domain_attr("gznl", lonvalue = real(zlon, dp),   &
          &                             latvalue = real(RESHAPE(plat(Nis0:Nie0, Njs0:Nje0),(/ Ni_0*Nj_0 /)),dp))
       CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj0glo)
diff --git a/src/OCE/IOM/iom_nf90.F90 b/src/OCE/IOM/iom_nf90.F90
index 171ef53bf..c49d9538e 100644
--- a/src/OCE/IOM/iom_nf90.F90
+++ b/src/OCE/IOM/iom_nf90.F90
@@ -40,8 +40,6 @@ MODULE iom_nf90
       MODULE PROCEDURE iom_nf90_rp0123d_dp
    END INTERFACE
 
-   !! * Substitutions
-#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
    !! $Id: iom_nf90.F90 14433 2021-02-11 08:06:49Z smasson $
@@ -546,7 +544,7 @@ CONTAINS
       INTEGER               :: idvar                ! variable id
       INTEGER               :: jd                   ! dimension loop counter
       INTEGER               :: ix1, ix2, iy1, iy2   ! subdomain indexes
-      INTEGER, DIMENSION(3) :: ishape               ! dimensions size
+      INTEGER, DIMENSION(4) :: idimsz               ! dimensions size
       INTEGER, DIMENSION(4) :: idimid               ! dimensions id
       CHARACTER(LEN=256)    :: clinfo               ! info character
       INTEGER               :: if90id               ! nf90 file identifier
@@ -629,9 +627,11 @@ CONTAINS
             itype = NF90_DOUBLE
          ENDIF
          IF( PRESENT(pv_r0d) ) THEN
-            CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype,                  iom_file(kiomid)%nvid(idvar) ), clinfo )
+            CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype,                    &
+               &                              iom_file(kiomid)%nvid(idvar) ), clinfo )
          ELSE
-            CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims), iom_file(kiomid)%nvid(idvar) ), clinfo )
+            CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims),   &
+               &                              iom_file(kiomid)%nvid(idvar) ), clinfo )
          ENDIF
          lchunk = .false.
          IF( snc4set%luse .AND. idims == 4 )   lchunk = .true.
@@ -673,13 +673,23 @@ CONTAINS
          ENDIF
          ! on what kind of domain must the data be written?
          IF( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) THEN
+            idimsz(1:2) = iom_file(kiomid)%dimsz(1:2,idvar)
+            IF(     idimsz(1) == Ni_0 .AND. idimsz(2) == Nj_0 ) THEN
+               ix1 = Nis0   ;   ix2 = Nie0   ;   iy1 = Njs0   ;   iy2 = Nje0
+            ELSEIF( idimsz(1) == jpi  .AND. idimsz(2) == jpj  ) THEN
+               ix1 = 1      ;   ix2 = jpi    ;   iy1 = 1      ;   iy2 = jpj
+            ELSEIF( idimsz(1) == jpi  .AND. idimsz(2) == jpj  ) THEN
+               ix1 = 1      ;   ix2 = jpi    ;   iy1 = 1      ;   iy2 = jpj
+            ELSE
+               CALL ctl_stop( 'iom_nf90_rp0123d: should have been an impossible case...' )
+            ENDIF
 
             ! write dimension variables if it is not already done
             ! =============
             ! trick: is defined to 0 => dimension variable are defined but not yet written
             IF( iom_file(kiomid)%dimsz(1, 4) == 0 ) THEN   ! time_counter = 0
-               CALL iom_nf90_check(    NF90_PUT_VAR( if90id, 1,                                      glamt(A2D(0)) ), clinfo )
-               CALL iom_nf90_check(    NF90_PUT_VAR( if90id, 2,                                      gphit(A2D(0)) ), clinfo )
+               CALL iom_nf90_check(    NF90_PUT_VAR( if90id, 1,                            glamt(ix1:ix2, iy1:iy2) ), clinfo )
+               CALL iom_nf90_check(    NF90_PUT_VAR( if90id, 2,                            gphit(ix1:ix2, iy1:iy2) ), clinfo )
                SELECT CASE (iom_file(kiomid)%comp)
                CASE ('OCE')
                   CALL iom_nf90_check( NF90_PUT_VAR( if90id, 3,                                           gdept_1d ), clinfo )
@@ -694,19 +704,6 @@ CONTAINS
                iom_file(kiomid)%dimsz(1, 4) = 1   ! so we don't enter this IF case any more...
                IF(lwp) WRITE(numout,*) TRIM(clinfo)//' write dimension variables done'
             ENDIF
-
-            IF( PRESENT(pv_r2d) )  ishape(1:2) = SHAPE(pv_r2d)
-            IF( PRESENT(pv_r3d) )  ishape(1:3) = SHAPE(pv_r3d)
-            IF(     ishape(1) == Ni_0   .AND. ishape(2) == Nj_0   ) THEN           ! array with 0 halo
-               ix1 = 1      ;   ix2 = Ni_0     ;   iy1 = 1      ;   iy2 = Nj_0
-            ELSEIF( ishape(1) == jpi    .AND. ishape(2) == jpj    ) THEN           ! array with nn_hls halos
-               ix1 = Nis0   ;   ix2 = Nie0     ;   iy1 = Njs0   ;   iy2 = Nje0
-            ELSEIF( ishape(1) == Ni_0+1 .AND. ishape(2) == Nj_0+1 ) THEN           ! nn_hls = 2 and array with 1 halo
-               ix1 = 2      ;   ix2 = Ni_0+1   ;   iy1 = 2      ;   iy2 = Nj_0+1
-            ELSE
-               CALL ctl_stop( 'iom_nf90_rp0123d: should have been an impossible case...' )
-            ENDIF
-
          ENDIF
 
          ! write the data
@@ -715,7 +712,7 @@ CONTAINS
             CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r0d                    ), clinfo )
          ELSEIF( PRESENT(pv_r1d) ) THEN
             CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r1d(:)                 ), clinfo )
-         ELSEIF( PRESENT(pv_r2d) ) THEN     
+         ELSEIF( PRESENT(pv_r2d) ) THEN
             CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r2d(ix1:ix2,iy1:iy2)   ), clinfo )
          ELSEIF( PRESENT(pv_r3d) ) THEN
             CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r3d(ix1:ix2,iy1:iy2,:) ), clinfo )
diff --git a/src/TOP/PISCES/P4Z/p4zlys.F90 b/src/TOP/PISCES/P4Z/p4zlys.F90
index 729daab9d..baf433510 100644
--- a/src/TOP/PISCES/P4Z/p4zlys.F90
+++ b/src/TOP/PISCES/P4Z/p4zlys.F90
@@ -32,7 +32,7 @@ MODULE p4zlys
    REAL(wp), PUBLIC ::   kdca   !: diss. rate constant calcite
    REAL(wp), PUBLIC ::   nca    !: order of reaction for calcite dissolution
 
-   INTEGER  ::   rmtss              ! number of seconds per month 
+   INTEGER  :: rmtss              ! number of seconds per month 
 
    !! * Module variables
    REAL(wp) :: calcon = 1.03E-2           !: mean calcite concentration [Ca2+]  in sea water [mole/kg solution]
@@ -67,11 +67,11 @@ CONTAINS
       REAL(wp) ::   zomegaca, zexcess, zexcess0, zkd
       CHARACTER (len=25) ::   charout
       REAL(wp), DIMENSION(A2D(0),jpk) ::   zco3, zcaldiss, zhinit, zhi, zco3sat
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)  :: zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )  CALL timing_start('p4z_lys')
       !
-
       zhinit  (:,:,:) = hi(:,:,:) / ( rhd(:,:,:) + 1._wp )
       !
       !     -------------------------------------------
@@ -126,10 +126,7 @@ CONTAINS
       !
 
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-          IF( l_dia ) THEN
-             ALLOCATE( zw3d(A2D(0),jpk) ) 
-             zw3d(A2D(0),jpk) = 0._wp
-          ENDIF
+             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) 
@@ -154,7 +151,7 @@ CONTAINS
              END_3D
              CALL iom_put( "DCAL", zw3d ) 
          ENDIF              
-         IF( l_dia ) DEALLOCATE( zw3d )
+         DEALLOCATE( zw3d )
       ENDIF
       !
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
@@ -205,6 +202,9 @@ CONTAINS
       ! Number of seconds per month 
       rmtss =  nyear_len(1) * rday / raamo
       !
+      ! CE not really needed ; tempory, shoub be removed when quotan( A2D(0),jpk )
+      excess(:,:,:) = 0._wp
+      !
    END SUBROUTINE p4z_lys_init
 
    !!======================================================================
diff --git a/src/TOP/PISCES/P4Z/p4zprod.F90 b/src/TOP/PISCES/P4Z/p4zprod.F90
index 00bbcfe8c..a48a6242a 100644
--- a/src/TOP/PISCES/P4Z/p4zprod.F90
+++ b/src/TOP/PISCES/P4Z/p4zprod.F90
@@ -83,7 +83,7 @@ CONTAINS
       REAL(wp), DIMENSION(A2D(0),jpk) :: zprorcan, zprorcad, zprofed, zprofen
       REAL(wp), DIMENSION(A2D(0),jpk) :: zpronewn, zpronewd
       REAL(wp), DIMENSION(A2D(0),jpk) :: zmxl_fac, zmxl_chl
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zpligprod
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zpligprod, zw3d
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_prod')
@@ -383,10 +383,6 @@ CONTAINS
        DEALLOCATE ( zw3d ) 
     ENDIF
     
-    IF( lk_iomput .AND.  knt == nrdttrc ) THEN
-       zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
-       !
-
     IF( lk_iomput .AND.  knt == nrdttrc ) THEN
        zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
        !
@@ -546,6 +542,9 @@ CONTAINS
       texcretn  = 1._wp - excretn
       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 
       !
    END SUBROUTINE p4z_prod_init
 
@@ -554,7 +553,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_prod_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( quotan(A2D(0),jpk), quotad(A2D(0),jpk), STAT = p4z_prod_alloc )
+      ALLOCATE( quotan(jpi,jpj,jpk), quotad(jpi,jpj,jpk), STAT = p4z_prod_alloc )
       !
       IF( p4z_prod_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_prod_alloc : failed to allocate arrays.' )
       !
diff --git a/src/TOP/PISCES/sms_pisces.F90 b/src/TOP/PISCES/sms_pisces.F90
index 58679aaeb..baf78a5d7 100644
--- a/src/TOP/PISCES/sms_pisces.F90
+++ b/src/TOP/PISCES/sms_pisces.F90
@@ -147,47 +147,47 @@ CONTAINS
       IF( ln_p4z .OR. ln_p5z ) THEN
 
          !* Optics
-         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) )
+         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) )
 
          !* Biological SMS
-         ALLOCATE( xksimax(jpi,jpj)  , biron(A2D(0),jpk)      ,  STAT=ierr(3) )
+         ALLOCATE( xksimax(jpi,jpj)  , biron(jpi,jpj,jpk)      ,  STAT=ierr(3) )
 
          ! Biological SMS
-         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) )
+         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) )
 
          !* Carbonate chemistry
-         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) )
+         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) )
          !
          !* Temperature dependency of SMS terms
-         ALLOCATE( tgfunc (A2D(0),jpk) , tgfunc2(A2D(0),jpk),   STAT=ierr(6) )
+         ALLOCATE( tgfunc (jpi,jpj,jpk) , tgfunc2(jpi,jpj,jpk),   STAT=ierr(6) )
          !
          !* Sinking speed
-         ALLOCATE( wsbio3 (A2D(0),jpk) , wsbio4 (A2D(0),jpk),   STAT=ierr(7) )
+         ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4 (jpi,jpj,jpk),   STAT=ierr(7) )
 
          !*  Size of phytoplankton cells
-         ALLOCATE( sizen (A2D(0),jpk), sized (A2D(0),jpk),        &
-           &       sizena(A2D(0),jpk), sizeda(A2D(0),jpk),      STAT=ierr(8) )
+         ALLOCATE( sizen (jpi,jpj,jpk), sized (jpi,jpj,jpk),        &
+           &       sizena(jpi,jpj,jpk), sizeda(jpi,jpj,jpk),      STAT=ierr(8) )
          ! 
-         ALLOCATE( plig(A2D(0),jpk)  ,                           STAT=ierr(9) )
+         ALLOCATE( plig(jpi,jpj,jpk)  ,                           STAT=ierr(9) )
       ENDIF
       !
       IF( ln_p5z ) THEN
          ! PISCES-QUOTA specific part      
-         ALLOCATE( epico(A2D(0),jpk)   , epicom(A2D(0),jpk) ,   STAT=ierr(10) ) 
+         ALLOCATE( epico(jpi,jpj,jpk)   , epicom(jpi,jpj,jpk) ,   STAT=ierr(10) ) 
 
          !*  Size of phytoplankton cells
-         ALLOCATE( sizep(A2D(0),jpk), sizepa(A2D(0),jpk),       STAT=ierr(11) )
+         ALLOCATE( sizep(jpi,jpj,jpk), sizepa(jpi,jpj,jpk),       STAT=ierr(11) )
       ENDIF
       !
       sms_pisces_alloc = MAXVAL( ierr )
diff --git a/src/TOP/trc.F90 b/src/TOP/trc.F90
index 10c40818b..6958508cc 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(A2D(0),jptra)  , trc_o(A2D(0),jptra)                        ,       &
+         &      trc_i(jpi,jpj,jptra)  , trc_o(jpi,jpj,jptra)                        ,       &
          &      trc_ice_ratio(jptra) , trc_ice_prescr(jptra) , cn_trc_o(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)                    ,       &  
+         &      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)                    ,       &  
          &      cvol(jpi,jpj,jpk)    , trai(jptra)                                ,       &
          &      ctrcnm(jptra)        , ctrcln(jptra)         , ctrcun(jptra)      ,       &
          &      ln_trc_ini(jptra)    ,                                                    &
@@ -175,9 +175,9 @@ CONTAINS
       !
       IF( ln_bdy       )   ALLOCATE( trcdta_bdy(jptra, jp_bdy)  , STAT = ierr(2) )
       !
-      IF (jp_dia3d > 0 )   ALLOCATE( trc3d(A2D(0),jpk,jp_dia3d), STAT = ierr(3) )
+      IF (jp_dia3d > 0 )   ALLOCATE( trc3d(jpi,jpj,jpk,jp_dia3d), STAT = ierr(3) )
       !
-      IF (jp_dia2d > 0 )   ALLOCATE( trc2d(A2D(0),jp_dia2d)    , STAT = ierr(4) )
+      IF (jp_dia2d > 0 )   ALLOCATE( trc2d(jpi,jpj,jp_dia2d)    , STAT = ierr(4) )
       ! 
       trc_alloc = MAXVAL( ierr )
       IF( trc_alloc /= 0 )   CALL ctl_stop( 'STOP', 'trc_alloc: failed to allocate arrays' )
diff --git a/src/TOP/trcwri.F90 b/src/TOP/trcwri.F90
index 8849c2256..abd64e4b3 100644
--- a/src/TOP/trcwri.F90
+++ b/src/TOP/trcwri.F90
@@ -58,8 +58,8 @@ CONTAINS
            WRITE(inum,*) clhstnam
            CLOSE(inum)
          ENDIF
-         
-         ALLOCATE( z3d(A2D(0),jpk) )  ;   z3d(A2D(0),:) = 0._wp
+
+         ALLOCATE( z3d(jpi,jpj,jpk) )  ;   z3d(:,:,:) = 0._wp
 
          ! Output of initial vertical scale factor
          CALL iom_put( "e3t_0", e3t_0(:,:,:) )
-- 
GitLab