diff --git a/cfgs/SHARED/field_def_nemo-oce.xml b/cfgs/SHARED/field_def_nemo-oce.xml
index 838cd2834cb3d8daad80c923b1f658e41fcd36e1..7b618ff67965e56b470d15decce18d951977bc48 100644
--- a/cfgs/SHARED/field_def_nemo-oce.xml
+++ b/cfgs/SHARED/field_def_nemo-oce.xml
@@ -615,12 +615,12 @@ that are available in the tidal-forcing implementation (see
     <field id="psiu_mle"     long_name="MLE streamfunction along i-axis"   unit="m3/s"   grid_ref="grid_U_3D" />
 
     <!-- uoce_eiv: available EIV (ln_ldfeiv=T and ln_ldfeiv_dia=T) -->
-    <field id="uoce_eiv"      long_name="EIV ocean current along i-axis"                                  standard_name="bolus_sea_water_x_velocity"                     unit="m/s"   grid_ref="grid_U_3D" />
-    <field id="ueiv_masstr"   long_name="EIV Ocean Mass X Transport"                                      standard_name="bolus_ocean_mass_x_transport"                   unit="kg/s"  grid_ref="grid_U_3D" />
-    <field id="ueiv_heattr"   long_name="ocean bolus heat transport along i-axis"                         standard_name="ocean_heat_x_transport_due_to_bolus_advection"  unit="W"                         />
-    <field id="ueiv_salttr"   long_name="ocean bolus salt transport along i-axis"                         standard_name="ocean_salt_x_transport_due_to_bolus_advection"  unit="Kg"                        />
-    <field id="ueiv_heattr3d" long_name="ocean bolus heat transport along i-axis"                         standard_name="ocean_heat_x_transport_due_to_bolus_advection"  unit="W"     grid_ref="grid_U_3D" />
-    <field id="ueiv_salttr3d" long_name="ocean bolus salt transport along i-axis"                         standard_name="ocean_salt_x_transport_due_to_bolus_advection"  unit="kg"    grid_ref="grid_U_3D" />
+    <field id="uoce_eiv"      long_name="EIV ocean current along i-axis"                                  standard_name="bolus_sea_water_x_velocity"                     unit="m/s"   grid_ref="grid_U_3D_inner" />
+    <field id="ueiv_masstr"   long_name="EIV Ocean Mass X Transport"                                      standard_name="bolus_ocean_mass_x_transport"                   unit="kg/s"  grid_ref="grid_U_3D_inner" />
+    <field id="ueiv_heattr"   long_name="ocean bolus heat transport along i-axis"                         standard_name="ocean_heat_x_transport_due_to_bolus_advection"  unit="W" grid_ref="grid_U_2D_inner"                        />
+    <field id="ueiv_salttr"   long_name="ocean bolus salt transport along i-axis"                         standard_name="ocean_salt_x_transport_due_to_bolus_advection"  unit="Kg" grid_ref="grid_U_2D_inner"                       />
+    <field id="ueiv_heattr3d" long_name="ocean bolus heat transport along i-axis"                         standard_name="ocean_heat_x_transport_due_to_bolus_advection"  unit="W"     grid_ref="grid_U_3D_inner" />
+    <field id="ueiv_salttr3d" long_name="ocean bolus salt transport along i-axis"                         standard_name="ocean_salt_x_transport_due_to_bolus_advection"  unit="kg"    grid_ref="grid_U_3D_inner" />
 
     <!-- uoce_bbl: available with ln_trabbl=T and nn_bbl_adv=1 -->
     <field id="uoce_bbl"     long_name="BBL ocean current along i-axis"    unit="m/s"  />
@@ -674,12 +674,12 @@ that are available in the tidal-forcing implementation (see
     <field id="psiv_mle"     long_name="MLE streamfunction along j-axis"   unit="m3/s"   grid_ref="grid_V_3D" />
 
     <!-- voce_eiv: available EIV (ln_ldfeiv=T and ln_ldfeiv_dia=T)  -->
-    <field id="voce_eiv"     long_name="EIV ocean current along j-axis"  standard_name="bolus_sea_water_y_velocity"     unit="m/s"   grid_ref="grid_V_3D" />
-    <field id="veiv_masstr"  long_name="EIV Ocean Mass Y Transport"      standard_name="bolus_ocean_mass_y_transport"   unit="kg/s"  grid_ref="grid_V_3D" />
-    <field id="veiv_heattr"   long_name="ocean bolus heat transport along j-axis"       standard_name="ocean_heat_y_transport_due_to_bolus_advection"   unit="W"                         />
-    <field id="veiv_salttr"   long_name="ocean bolus salt transport along j-axis"       standard_name="ocean_salt_x_transport_due_to_bolus_advection"   unit="Kg"                        />
-    <field id="veiv_heattr3d" long_name="ocean bolus heat transport along j-axis"       standard_name="ocean_heat_y_transport_due_to_bolus_advection"   unit="W"    grid_ref="grid_V_3D" />
-    <field id="veiv_salttr3d" long_name="ocean bolus salt transport along j-axis"       standard_name="ocean_salt_y_transport_due_to_bolus_advection"   unit="kg"   grid_ref="grid_V_3D" />
+    <field id="voce_eiv"     long_name="EIV ocean current along j-axis"  standard_name="bolus_sea_water_y_velocity"     unit="m/s"   grid_ref="grid_V_3D_inner" />
+    <field id="veiv_masstr"  long_name="EIV Ocean Mass Y Transport"      standard_name="bolus_ocean_mass_y_transport"   unit="kg/s"  grid_ref="grid_V_3D_inner" />
+    <field id="veiv_heattr"   long_name="ocean bolus heat transport along j-axis"       standard_name="ocean_heat_y_transport_due_to_bolus_advection"   unit="W"  grid_ref="grid_V_2D_inner"   />
+    <field id="veiv_salttr"   long_name="ocean bolus salt transport along j-axis"       standard_name="ocean_salt_x_transport_due_to_bolus_advection"   unit="Kg" grid_ref="grid_V_2D_inner"   />
+    <field id="veiv_heattr3d" long_name="ocean bolus heat transport along j-axis"       standard_name="ocean_heat_y_transport_due_to_bolus_advection"   unit="W"    grid_ref="grid_V_3D_inner" />
+    <field id="veiv_salttr3d" long_name="ocean bolus salt transport along j-axis"       standard_name="ocean_salt_y_transport_due_to_bolus_advection"   unit="kg"   grid_ref="grid_V_3D_inner" />
 
 
     <!-- voce_bbl: available with ln_trabbl=T and nn_bbl_adv=1 -->
@@ -713,8 +713,8 @@ that are available in the tidal-forcing implementation (see
     <field id="wocetr_eff"   long_name="effective ocean vertical transport"                                                                   unit="m3/s" />
 
     <!-- woce_eiv: available with EIV  (ln_ldfeiv=T and ln_ldfeiv_dia=T)  -->
-    <field id="woce_eiv"     long_name="EIV ocean vertical velocity"                    standard_name="bolus_upward_sea_water_velocity"       unit="m/s"  />
-    <field id="weiv_masstr"  long_name="EIV Upward Ocean Mass Transport"  standard_name="bolus_upward_ocean_mass_transport"             unit="kg/s"   />
+    <field id="woce_eiv"     long_name="EIV ocean vertical velocity"                    standard_name="bolus_upward_sea_water_velocity"       unit="m/s" grid_ref="grid_W_3D_inner"  />
+    <field id="weiv_masstr"  long_name="EIV Upward Ocean Mass Transport"  standard_name="bolus_upward_ocean_mass_transport"             unit="kg/s" grid_ref="grid_W_3D_inner"  />
     <field id="weiv_heattr3d" long_name="ocean bolus heat transport"    standard_name="ocean_heat_z_transport_due_to_bolus_advection"   unit="W"    />
     <field id="weiv_salttr3d" long_name="ocean bolus salt transport"    standard_name="ocean_salt_z_transport_due_to_bolus_advection"   unit="kg"   />
 
diff --git a/src/OCE/LDF/ldfdyn.F90 b/src/OCE/LDF/ldfdyn.F90
index 92623ba08b2e2271515a09c0a65ee12dd7ce2111..8f1e38f26ed1c596f2cce40c3e9d298036c53602 100644
--- a/src/OCE/LDF/ldfdyn.F90
+++ b/src/OCE/LDF/ldfdyn.F90
@@ -321,10 +321,10 @@ CONTAINS
             l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
             !
             !                          ! allocate arrays used in ldf_dyn. 
-            ALLOCATE( dtensq(jpi,jpj,jpk) , dshesq(jpi,jpj,jpk) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr )
+            ALLOCATE( dtensq(A2D(1),jpk) , dshesq(A2D(1),jpk) , esqt(A2D(0)) , esqf(A2D(0)) , STAT=ierr )
             IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays')
             !
-            DO_2D( 1, 1, 1, 1 )        ! Set local gridscale values
+            DO_2D( 0, 0, 0, 0 )        ! Set local gridscale values
                esqt(ji,jj) = ( 2._wp * e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 
                esqf(ji,jj) = ( 2._wp * e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 
             END_2D
@@ -419,7 +419,7 @@ CONTAINS
             !                                                                       ! of |U|L^3/16 in blp case
             DO jk = 1, jpkm1
                !
-               DO_2D( 0, 0, 0, 0 )
+               DO_2D( 0, 1, 0, 1 )
                   zdb =    ( uu(ji,jj,jk,Kbb) * r1_e2u(ji,jj) -  uu(ji-1,jj,jk,Kbb) * r1_e2u(ji-1,jj) )  &
                        &                      * r1_e1t(ji,jj) * e2t(ji,jj)                           &
                        & - ( vv(ji,jj,jk,Kbb) * r1_e1v(ji,jj) -  vv(ji,jj-1,jk,Kbb) * r1_e1v(ji,jj-1) )  &
@@ -437,8 +437,6 @@ CONTAINS
                !
             END DO
             !
-            CALL lbc_lnk( 'ldfdyn', dtensq, 'T', 1.0_wp )  ! lbc_lnk on dshesq not needed
-            !
             DO jk = 1, jpkm1
               !
                DO_2D( 0, 0, 0, 0 )                                   ! T-point value
@@ -455,7 +453,7 @@ CONTAINS
                   !
                END_2D
                !
-               DO_2D( 1, 0, 1, 0 )                                   ! F-point value
+               DO_2D( 0, 0, 0, 0 )                                   ! F-point value
                   !
                   zu2pv2_ij_p1 = uu(ji  ,jj+1,jk, kbb) * uu(ji  ,jj+1,jk, kbb) + vv(ji+1,jj  ,jk, kbb) * vv(ji+1,jj  ,jk, kbb)
                   zu2pv2_ij    = uu(ji  ,jj  ,jk, kbb) * uu(ji  ,jj  ,jk, kbb) + vv(ji  ,jj  ,jk, kbb) * vv(ji  ,jj  ,jk, kbb)
diff --git a/src/OCE/LDF/ldfslp.F90 b/src/OCE/LDF/ldfslp.F90
index d4804673554de52a04ebd0b8dd6de52161314407..dacc75fc0fae325933e745d80ef679fc01b8c266 100644
--- a/src/OCE/LDF/ldfslp.F90
+++ b/src/OCE/LDF/ldfslp.F90
@@ -113,14 +113,12 @@ CONTAINS
       REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
       !!
       INTEGER  ::   ji , jj , jk    ! dummy loop indices
-      INTEGER  ::   ii0, ii1        ! temporary integer
-      INTEGER  ::   ij0, ij1        ! temporary integer
       REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw, z1_slpmax ! local scalars
       REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
       REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
       REAL(wp) ::   zck, zfk,      zbw             !   -      -
       REAL(wp) ::   zdepu, zdepv                   !   -      -
-      REAL(wp), DIMENSION(jpi,jpj)     ::  zslpml_hmlpu, zslpml_hmlpv
+      REAL(wp), DIMENSION(A2D(0))     ::  zslpml_hmlpu, zslpml_hmlpv
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zgru, zwz, zdzr
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zgrv, zww
       !!----------------------------------------------------------------------
@@ -593,8 +591,6 @@ CONTAINS
       !
       uslpml (1,:) = 0._wp      ;      uslpml (jpi,:) = 0._wp
       vslpml (1,:) = 0._wp      ;      vslpml (jpi,:) = 0._wp
-      wslpiml(1,:) = 0._wp      ;      wslpiml(jpi,:) = 0._wp
-      wslpjml(1,:) = 0._wp      ;      wslpjml(jpi,:) = 0._wp
       !
       !                                            !==   surface mixed layer mask   !
       DO_3D( 1, 1, 1, 1, 1, jpk )                  ! =1 inside the mixed layer, =0 otherwise
@@ -657,8 +653,6 @@ CONTAINS
          wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik)
          wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik)
       END_2D
-      !!gm this lbc_lnk should be useless....
-      CALL lbc_lnk( 'ldfslp', uslpml , 'U', -1.0_wp , vslpml , 'V', -1.0_wp , wslpiml, 'W', -1.0_wp , wslpjml, 'W', -1.0_wp ) 
       !
    END SUBROUTINE ldf_slp_mxl
 
@@ -700,8 +694,8 @@ CONTAINS
       ELSE                             ! Madec operator : slopes at u-, v-, and w-points
          IF(lwp) WRITE(numout,*) '   ==>>>   iso operator (Madec)'
          ALLOCATE( omlmask(jpi,jpj,jpk) ,                                                                        &
-            &      uslp(jpi,jpj,jpk) , uslpml(jpi,jpj) , wslpi(jpi,jpj,jpk) , wslpiml(jpi,jpj) ,     &
-            &      vslp(jpi,jpj,jpk) , vslpml(jpi,jpj) , wslpj(jpi,jpj,jpk) , wslpjml(jpi,jpj) , STAT=ierr )
+            &      uslp(jpi,jpj,jpk) , uslpml(A2D(0)) , wslpi(jpi,jpj,jpk) , wslpiml(A2D(0)) ,     &
+            &      vslp(jpi,jpj,jpk) , vslpml(A2D(0)) , wslpj(jpi,jpj,jpk) , wslpjml(A2D(0)) , STAT=ierr )
          IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' )
 
          ! Direction of lateral diffusion (tracers and/or momentum)
diff --git a/src/OCE/LDF/ldftra.F90 b/src/OCE/LDF/ldftra.F90
index 87e1d746557589acb22e01749b03156bbc5172cb..dea21729609f0d6bd405cbbcc9e8014202f7f362 100644
--- a/src/OCE/LDF/ldftra.F90
+++ b/src/OCE/LDF/ldftra.F90
@@ -794,8 +794,8 @@ CONTAINS
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
       REAL(wp) ::   zztmp   ! local scalar
-      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zw2d   ! 2D workspace
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zw3d   ! 3D workspace
+      REAL(wp), DIMENSION(A2D(0))     ::   zw2d   ! 2D workspace
+      REAL(wp), DIMENSION(A2D(0),jpk) ::   zw3d   ! 3D workspace
       !!----------------------------------------------------------------------
       !
 !!gm I don't like this routine....   Crazy  way of doing things, not optimal at all...
@@ -828,27 +828,27 @@ CONTAINS
          DO_2D( 0, 0, 0, 0 )
             zw2d(ji,jj) = rho0 * e1e2t(ji,jj)
          END_2D
-         DO jk = 1, jpk
-            zw3d(:,:,jk) = zw3d(:,:,jk) * zw2d(:,:)
-         END DO
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * zw2d(ji,jj) 
+         END_3D
          CALL iom_put( "weiv_masstr" , zw3d )
       ENDIF
       !
       IF( iom_use('ueiv_masstr') ) THEN
-         zw3d(:,:,:) = 0.e0
-         DO jk = 1, jpkm1
-            zw3d(:,:,jk) = rho0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) )
-         END DO
+         zw3d(:,:,jpk) = 0.e0
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zw3d(ji,jj,jk) = rho0 * ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) )
+         END_3D
          CALL iom_put( "ueiv_masstr", zw3d )                  ! mass transport in i-direction
       ENDIF
       !
       zztmp = 0.5_wp * rho0 * rcp
       IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
         zw2d(:,:)   = 0._wp
-        zw3d(:,:,:) = 0._wp
+        zw3d(:,:,jpk) = 0._wp
         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-           zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
-              &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji+1,jj,jk,jp_tem,Kmm) )
+           zw3d(ji,jj,jk) = ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
+              &           * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji+1,jj,jk,jp_tem,Kmm) )
            zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
         END_3D
         CALL iom_put( "ueiv_heattr"  , zztmp * zw2d )                  ! heat transport in i-direction
@@ -856,18 +856,18 @@ CONTAINS
       ENDIF
       !
       IF( iom_use('veiv_masstr') ) THEN
-         zw3d(:,:,:) = 0.e0
-         DO jk = 1, jpkm1
-            zw3d(:,:,jk) = rho0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) )
-         END DO
+         zw3d(:,:,jpk) = 0.e0
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zw3d(ji,jj,jk) = rho0 * ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) )
+         END_3D
          CALL iom_put( "veiv_masstr", zw3d )                  ! mass transport in i-direction
       ENDIF
       !
       zw2d(:,:)   = 0._wp
-      zw3d(:,:,:) = 0._wp
+      zw3d(:,:,jpk) = 0._wp
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-         zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
-            &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji,jj+1,jk,jp_tem,Kmm) )
+         zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
+            &           * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji,jj+1,jk,jp_tem,Kmm) )
          zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
       END_3D
       CALL iom_put( "veiv_heattr"  , zztmp * zw2d )                  !  heat transport in j-direction
@@ -878,20 +878,20 @@ CONTAINS
       zztmp = 0.5_wp * 0.5
       IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN
         zw2d(:,:) = 0._wp
-        zw3d(:,:,:) = 0._wp
+        zw3d(:,:,jpk) = 0._wp
         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-           zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
-              &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji+1,jj,jk,jp_sal,Kmm) )
+           zw3d(ji,jj,jk) = ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
+              &           * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji+1,jj,jk,jp_sal,Kmm) )
            zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
         END_3D
         CALL iom_put( "ueiv_salttr", zztmp * zw2d )                  ! salt transport in i-direction
         CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                ! salt transport in i-direction
       ENDIF
       zw2d(:,:) = 0._wp
-      zw3d(:,:,:) = 0._wp
+      zw3d(:,:,jpk) = 0._wp
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-         zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
-            &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji,jj+1,jk,jp_sal,Kmm) )
+         zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
+            &           * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji,jj+1,jk,jp_sal,Kmm) )
          zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
       END_3D
       CALL iom_put( "veiv_salttr"  , zztmp * zw2d )                  !  salt transport in j-direction