diff --git a/src/OCE/IOM/iom.F90 b/src/OCE/IOM/iom.F90
index 47cb14c850b10be7b1e83380f4be98c8e59a75ef..52431da3c8bbecc53ffe4eada6d017d3dc83c53c 100644
--- a/src/OCE/IOM/iom.F90
+++ b/src/OCE/IOM/iom.F90
@@ -1327,6 +1327,7 @@ 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)
@@ -1343,6 +1344,8 @@ CONTAINS
                ELSE
                   CALL ctl_stop( 'iom_get_123d: should have been an impossible case...' )
                ENDIF
+               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)
@@ -1355,11 +1358,6 @@ 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...
@@ -2346,14 +2344,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) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
+            & tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 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) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
+            & tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 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(:))
diff --git a/src/OCE/IOM/iom_nf90.F90 b/src/OCE/IOM/iom_nf90.F90
index c43809493c55fc6f472d1c2fdd029e1377d39c47..67997b9a407e6a19a201079d3187b706460de8a5 100644
--- a/src/OCE/IOM/iom_nf90.F90
+++ b/src/OCE/IOM/iom_nf90.F90
@@ -40,6 +40,8 @@ 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 $
@@ -544,7 +546,7 @@ CONTAINS
       INTEGER               :: idvar                ! variable id
       INTEGER               :: jd                   ! dimension loop counter
       INTEGER               :: ix1, ix2, iy1, iy2   ! subdomain indexes
-      INTEGER, DIMENSION(4) :: idimsz               ! dimensions size
+      INTEGER, DIMENSION(3) :: ishape               ! dimensions size
       INTEGER, DIMENSION(4) :: idimid               ! dimensions id
       CHARACTER(LEN=256)    :: clinfo               ! info character
       INTEGER               :: if90id               ! nf90 file identifier
@@ -627,11 +629,9 @@ 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,23 +673,13 @@ 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(ix1:ix2, iy1:iy2) ), clinfo )
-               CALL iom_nf90_check(    NF90_PUT_VAR( if90id, 2,                            gphit(ix1:ix2, iy1:iy2) ), clinfo )
+               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 )
                SELECT CASE (iom_file(kiomid)%comp)
                CASE ('OCE')
                   CALL iom_nf90_check( NF90_PUT_VAR( if90id, 3,                                           gdept_1d ), clinfo )
@@ -704,6 +694,17 @@ 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
+               ix1 = 1      ;   ix2 = Ni_0   ;   iy1 = 1      ;   iy2 = Nj_0
+            ELSEIF( ishape(1) == jpi  .AND. ishape(2) == jpj  ) THEN
+               ix1 = Nis0   ;   ix2 = Nie0   ;   iy1 = Njs0   ;   iy2 = Nje0
+            ELSE
+               CALL ctl_stop( 'iom_nf90_rp0123d: should have been an impossible case...' )
+            ENDIF
+
          ENDIF
 
          ! write the data
@@ -712,7 +713,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 )