diff --git a/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg b/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
index 4b17c10299ec87c10936d1ed936390a783d0f3e3..c6e94339d09e11830e49e55ebf1aafcd7761c391 100644
--- a/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
+++ b/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
@@ -118,6 +118,17 @@
 &namsbc_abl    !   Atmospheric Boundary Layer formulation           (ln_abl = T)
 !-----------------------------------------------------------------------
    cn_dom           = 'dom_cfg_abl_L25Z10'
+
+   nn_dyn_restore = 2         ! restoring option for dynamical ABL variables: = 0 no restoring
+                              !                                               = 1 equatorial restoring
+                              !                                               = 2 global restoring
+                              !                                               = 1 equatorial restoring
+                              !                                               = 2 global restoring
+   rn_ldyn_min   = 24.        ! dynamics nudging magnitude inside the ABL [hour] (~3 rn_Dt)
+   rn_ldyn_max   =  6.        ! dynamics nudging magnitude above  the ABL [hour] (~1 rn_Dt)
+   rn_ltra_min   = 24.        ! tracers  nudging magnitude inside the ABL [hour] (~3 rn_Dt)
+   rn_ltra_max   =  6.        ! tracers  nudging magnitude above  the ABL [hour] (~1 rn_Dt)
+   rn_vfac       =  1.
 /
 !-----------------------------------------------------------------------
 &namtra_qsr    !   penetrative solar radiation                          (ln_traqsr =T)
diff --git a/cfgs/SHARED/grid_def_nemo.xml b/cfgs/SHARED/grid_def_nemo.xml
index 374f67c6544360e11ac7660a6bd654850f59c3ed..0120bb3cce47eaf6f49af406da56a241a64b791e 100644
--- a/cfgs/SHARED/grid_def_nemo.xml
+++ b/cfgs/SHARED/grid_def_nemo.xml
@@ -318,14 +318,14 @@
 
   <!-- ABL grid definition -->
   <grid id="grid_TA_2D">
-    <domain domain_ref="grid_T" />
+    <domain domain_ref="grid_T_inner" />
   </grid>
   <grid id="grid_TA_3D">
-    <domain domain_ref="grid_T" />
+    <domain domain_ref="grid_T_inner" />
     <axis  id="ght_abl" />
   </grid>
   <grid id="grid_WA_3D">
-    <domain domain_ref="grid_T" />
+    <domain domain_ref="grid_T_inner" />
     <axis  id="ghw_abl" />
   </grid>
   <!--  -->
diff --git a/cfgs/SHARED/namelist_ref b/cfgs/SHARED/namelist_ref
index 28a67e4d54a5df21cc92367afbd896ee85a385ab..00688501d60bce2d58b1bdd7985b60939225587f 100644
--- a/cfgs/SHARED/namelist_ref
+++ b/cfgs/SHARED/namelist_ref
@@ -312,7 +312,7 @@
 &namsbc_abl    !   Atmospheric Boundary Layer formulation           (ln_abl = T)
 !-----------------------------------------------------------------------
    cn_dir           = './'      !  root directory for the location of the ABL grid file
-   cn_dom           = 'dom_cfg_abl.nc'
+   cn_dom           = 'dom_cfg_abl'
 
    cn_ablrst_in     = "restart_abl"   !  suffix of abl restart name (input)
    cn_ablrst_out    = "restart_abl"   !  suffix of abl restart name (output)
@@ -326,11 +326,11 @@
    nn_dyn_restore = 0         ! restoring option for dynamical ABL variables: = 0 no restoring
                               !                                               = 1 equatorial restoring
                               !                                               = 2 global restoring
-   rn_vfac       = 0.
-   rn_ldyn_min   =  4.5       ! dynamics nudging magnitude inside the ABL [hour] (~3 rn_Dt)
-   rn_ldyn_max   =  1.5       ! dynamics nudging magnitude above  the ABL [hour] (~1 rn_Dt)
-   rn_ltra_min   =  4.5       ! tracers  nudging magnitude inside the ABL [hour] (~3 rn_Dt)
-   rn_ltra_max   =  1.5       ! tracers  nudging magnitude above  the ABL [hour] (~1 rn_Dt)
+   rn_ldyn_min   =  0.        ! dynamics nudging magnitude inside the ABL [hour] (~3 rn_Dt)
+   rn_ldyn_max   =  0.        ! dynamics nudging magnitude above  the ABL [hour] (~1 rn_Dt)
+   rn_ltra_min   =  0.        ! tracers  nudging magnitude inside the ABL [hour] (~3 rn_Dt)
+   rn_ltra_max   =  0.        ! tracers  nudging magnitude above  the ABL [hour] (~1 rn_Dt)
+   rn_vfac       =  0.
    nn_amxl       =  0         ! mixing length: = 0 Deardorff 80 length-scale
                               !                = 1 length-scale based on the distance to the PBL height
                               !                = 2 Bougeault & Lacarrere 89 length-scale
diff --git a/makenemo b/makenemo
index 8f429bd2db3ce33385783bcae5658d7ff39d89a1..a7c8bd9068c7671682912b7d11f0088f77698616 100755
--- a/makenemo
+++ b/makenemo
@@ -223,8 +223,8 @@ done
 
 if [ -n "$x_name" ]   # is the configuration existing in cfgs or tests? or is it a new conf?
 then
-    incfg=$( find $CFGS_DIR   -type d -name $x_name | wc -l )   # this configuration exists in CFGS_DIR
-    intst=$( find $TESTS_DIR  -type d -name $x_name | wc -l )   # this configuration exists in TESTS_DIR
+    incfg=$( find $CFGS_DIR  -type d -name $x_name | wc -l )   # this configuration exists in CFGS_DIR
+    intst=$( find $TESTS_DIR -type d -name $x_name | wc -l )   # this configuration exists in TESTS_DIR
     if [ $(( $incfg + $intst )) -eq 2 ] ; then
 	echo -e "\033[0;31m\nERROR: the $x_name directory is found twice: in $CFGS_DIR and in $TESTS_DIR\033[0m\n"
 	exit 4
diff --git a/sette/BATCH_TEMPLATE/batch-X64_AA_INTEL_OMPI b/sette/BATCH_TEMPLATE/batch-X64_AA_INTEL_OMPI
index bdd075c6c25cfc14359dcc9eced3b355eed74968..3eca6d7c613a4911733c6577b9a0bc0a05442565 100644
--- a/sette/BATCH_TEMPLATE/batch-X64_AA_INTEL_OMPI
+++ b/sette/BATCH_TEMPLATE/batch-X64_AA_INTEL_OMPI
@@ -98,7 +98,7 @@ EOF
 #
 # Run SPMD case
 #
-    time ./nemo
+    $MPIRUN --ntasks=TOTAL_NPROCS ./nemo
   fi
 #
 
diff --git a/sette/all_functions.sh b/sette/all_functions.sh
index 8bc6e975c89614583364b1080507f0794b3303e2..8737aa4bef0b7e57b99981927aa4d6e383640beb 100755
--- a/sette/all_functions.sh
+++ b/sette/all_functions.sh
@@ -172,7 +172,7 @@ clean_config() {
 # define validation dir
 set_valid_dir () {
     if [ ${DETACHED_HEAD} == "no" ] ; then
-      branchname=$( git branch --show-current )
+      branchname=$( git -C ${MAIN_DIR} branch --show-current )
       REVISION_NB=$( git -C ${MAIN_DIR} rev-list --abbrev-commit origin/$branchname | head -1l )
     else
       REVISION_NB=${DETACHED_CMIT}
diff --git a/sette/sette.sh b/sette/sette.sh
index 7537af5e849258a035b630a4dd8d78f7b3f0a852..cc2846f44373a5cb081c662fb08088aa3a031427 100755
--- a/sette/sette.sh
+++ b/sette/sette.sh
@@ -40,10 +40,10 @@ export USER_INPUT='yes'        # Default: yes => request user input on decisions
 #
 # Check that git branch is usable
 export DETACHED_HEAD="no"
-git branch --show-current >& /dev/null
+git -C ${MAIN_DIR} branch --show-current >& /dev/null
 if [[ $? == 0 ]] ; then
   # subdirectory below NEMO_VALIDATION_DIR defaults to branchname
-  export SETTE_SUB_VAL="$(git branch --show-current)"
+  export SETTE_SUB_VAL="$(git -C ${MAIN_DIR} branch --show-current)"
   if [ -z $SETTE_SUB_VAL ] ; then
    # Probabably on a detached HEAD (possibly testing an old commit).
    # Verify this and try to recover original commit
@@ -76,6 +76,7 @@ if [ $# -gt 0 ]; then
         y) dry_run=1
            echo "";;
         b) NEMO_DEBUG="-b"
+           echo "-b: Nemo will be compiled with DEBUG options"
            echo "";;
         r) NO_REPORT=1
            echo "";;
diff --git a/src/ABL/abl.F90 b/src/ABL/abl.F90
index ccffa088d6e3caddc643458ab076c440454fccf5..c89f0c1b6e4875be1d00fc64bd0dbcf99acbbe98 100644
--- a/src/ABL/abl.F90
+++ b/src/ABL/abl.F90
@@ -41,6 +41,8 @@ MODULE abl
    !
    INTEGER , PUBLIC :: nt_n, nt_a       !: now / after indices (equal 1 or 2)
    ! 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
    !! $Id: abl.F90 4990 2014-12-15 16:42:49Z timgraham $ 
@@ -55,18 +57,18 @@ CONTAINS
       INTEGER :: ierr
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( u_abl   (1:jpi,1:jpj,1:jpka,jptime     ), &
-         &      v_abl   (1:jpi,1:jpj,1:jpka,jptime     ), &
-         &      tq_abl  (1:jpi,1:jpj,1:jpka,jptime,jptq), &
-         &      tke_abl (1:jpi,1:jpj,1:jpka,jptime     ), &
-         &      avm_abl (1:jpi,1:jpj,1:jpka            ), &
-         &      avt_abl (1:jpi,1:jpj,1:jpka            ), &
-         &      mxld_abl(1:jpi,1:jpj,1:jpka            ), &
-         &      mxlm_abl(1:jpi,1:jpj,1:jpka            ), &
-         &      fft_abl (1:jpi,1:jpj                   ), &
-         &      pblh    (1:jpi,1:jpj                   ), &
-         &      msk_abl (1:jpi,1:jpj                   ), &
-         &      rest_eq (1:jpi,1:jpj                   ), &
+      ALLOCATE( u_abl   (A2D(0),1:jpka,jptime     ), &
+         &      v_abl   (A2D(0),1:jpka,jptime     ), &
+         &      tq_abl  (A2D(0),1:jpka,jptime,jptq), &
+         &      tke_abl (A2D(0),1:jpka,jptime     ), &
+         &      avm_abl (A2D(0),1:jpka            ), &
+         &      avt_abl (A2D(0),1:jpka            ), &
+         &      mxld_abl(A2D(0),1:jpka            ), &
+         &      mxlm_abl(A2D(0),1:jpka            ), &
+         &      fft_abl (A2D(0)                   ), &
+         &      pblh    (A2D(1)                   ), &   ! needed for optional smoothing
+         &      msk_abl (A2D(0)                   ), &
+         &      rest_eq (A2D(0)                   ), &
          &      e3t_abl (1:jpka), e3w_abl(1:jpka)       , &
          &      ght_abl (1:jpka), ghw_abl(1:jpka)       , STAT=ierr )
          !
diff --git a/src/ABL/ablmod.F90 b/src/ABL/ablmod.F90
index ffaddbbcd36d7635404b1ef1d294790d4846ea31..84a2a99d389cdd865aa3c9de3dc9d91c53ba637b 100644
--- a/src/ABL/ablmod.F90
+++ b/src/ABL/ablmod.F90
@@ -67,54 +67,50 @@ CONTAINS
       !!              - Finalize flux computation in psen, pevp, pwndm, ptaui, ptauj, ptaum
       !!
       !!----------------------------------------------------------------------
-      INTEGER  , INTENT(in   )                   ::   kt         ! time step index
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psst       ! sea-surface temperature [Celsius]
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu       ! sea-surface u (U-point)
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv       ! sea-surface v (V-point)
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq       ! sea-surface humidity
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pu_dta     ! large-scale windi
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pv_dta     ! large-scale windj
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pgu_dta    ! large-scale hpgi
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pgv_dta    ! large-scale hpgj
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pt_dta     ! large-scale pot. temp.
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:,:) ::   pq_dta     ! large-scale humidity
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pslp_dta   ! sea-level pressure
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du     ! Cd x Du (T-point)
-      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   psen       ! Ch x Du
-      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pevp       ! Ce x Du
-      REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pwndm      ! ||uwnd - uoce||
-      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   plat       ! latent heat flux
-      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui      ! taux
-      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj      ! tauy
-      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaum      ! ||tau||
-      !
+
+      INTEGER  , INTENT(in   )                           ::   kt         ! time step index
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0))        ::   psst       ! sea-surface temperature [Celsius]
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(nn_hls))   ::   pssu       ! sea-surface u (U-point)
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(nn_hls))   ::   pssv       ! sea-surface v (V-point)
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0))        ::   pssq       ! sea-surface humidity
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0),1:jpka) ::   pu_dta     ! large-scale windi
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0),1:jpka) ::   pv_dta     ! large-scale windj
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0),1:jpka) ::   pgu_dta    ! large-scale hpgi
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0),1:jpka) ::   pgv_dta    ! large-scale hpgj
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0),1:jpka) ::   pt_dta     ! large-scale pot. temp.
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0),1:jpka) ::   pq_dta     ! large-scale humidity
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0))        ::   pslp_dta   ! sea-level pressure
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0))        ::   pcd_du     ! Cd x Du (T-point)
+      REAL(wp) , INTENT(inout), DIMENSION(A2D(0))        ::   psen       ! Ch x Du
+      REAL(wp) , INTENT(inout), DIMENSION(A2D(0))        ::   pevp       ! Ce x Du
+      REAL(wp) , INTENT(inout), DIMENSION(A2D(0))        ::   pwndm      ! ||uwnd - uoce||
+      REAL(wp) , INTENT(  out), DIMENSION(A2D(0))        ::   plat       ! latent heat flux
+      REAL(wp) , INTENT(  out), DIMENSION(A2D(nn_hls))   ::   ptaui      ! taux
+      REAL(wp) , INTENT(  out), DIMENSION(A2D(nn_hls))   ::   ptauj      ! tauy
+      REAL(wp) , INTENT(  out), DIMENSION(A2D(0))        ::   ptaum      ! ||tau||
 #if defined key_si3
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptm_su       ! ice-surface temperature [K]
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu_ice     ! ice-surface u (U-point)
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv_ice     ! ice-surface v (V-point)
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq_ice     ! ice-surface humidity
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du_ice   ! Cd x Du over ice (T-point)
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psen_ice     ! Ch x Du over ice (T-point)
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pevp_ice     ! Ce x Du over ice (T-point)
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndm_ice    ! ||uwnd - uice||
-      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pfrac_oce    ! ocean fraction
-      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui_ice    ! ice-surface taux stress (T-point)
-      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (T-point)
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0))        ::   ptm_su     ! ice-surface temperature [K]
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(nn_hls))   ::   pssu_ice   ! ice-surface u (U-point)
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(nn_hls))   ::   pssv_ice   ! ice-surface v (V-point)
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0))        ::   pssq_ice   ! ice-surface humidity
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0))        ::   pcd_du_ice ! Cd x Du over ice (T-point)
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0))        ::   psen_ice   ! Ch x Du over ice (T-point)
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0))        ::   pevp_ice   ! Ce x Du over ice (T-point)
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0))        ::   pwndm_ice  ! ||uwnd - uice||
+      REAL(wp) , INTENT(in   ), DIMENSION(A2D(0))        ::   pfrac_oce  ! ocean fraction
+      REAL(wp) , INTENT(  out), DIMENSION(A2D(nn_hls))   ::   ptaui_ice  ! ice-surface taux stress (T-point)
+      REAL(wp) , INTENT(  out), DIMENSION(A2D(nn_hls))   ::   ptauj_ice  ! ice-surface tauy stress (T-point)
 #endif
       !
-      REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   zwnd_i, zwnd_j
-      REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   zsspt
-      REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   ztabs, zpre
-      REAL(wp), DIMENSION(1:jpi      ,2:jpka)    ::   zCF
-      !
-      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_a
-      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_b
-      REAL(wp), DIMENSION(1:jpi      ,1:jpka)    ::   z_elem_c
+      REAL(wp), DIMENSION(A2D(0))         ::   zwnd_i, zwnd_j
+      REAL(wp), DIMENSION(A2D(0))         ::   zsspt, ztabs, zpre
+      REAL(wp), DIMENSION(A1Di(0),2:jpka) ::   zCF
+      REAL(wp), DIMENSION(A1Di(0),1:jpka) ::   z_elem_a, z_elem_b, z_elem_c
       !
-      INTEGER             ::   ji, jj, jk, jtra, jbak               ! dummy loop indices
-      REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2
-      REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce
-      LOGICAL             ::   SemiImp_Cor = .TRUE.
+      INTEGER  ::   ji, jj, jk, jtra, jbak               ! dummy loop indices
+      REAL(wp) ::   zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2
+      REAL(wp) ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce
+      LOGICAL  ::   SemiImp_Cor = .TRUE.
       !
       !!---------------------------------------------------------------------
       !
@@ -124,12 +120,12 @@ CONTAINS
          WRITE(numout,*) '~~~~~~'
       ENDIF
       !
-      IF( kt == nit000 ) ALLOCATE ( ustar2( 1:jpi, 1:jpj ) )
-      IF( kt == nit000 ) ALLOCATE ( zrough( 1:jpi, 1:jpj ) )
+      IF( kt == nit000 ) ALLOCATE ( ustar2( A2D(0) ) )
+      IF( kt == nit000 ) ALLOCATE ( zrough( A2D(0) ) )
       !! Compute ustar squared as Cd || Uatm-Uoce ||^2
       !! needed for surface boundary condition of TKE
       !! pwndm contains | U10m - U_oce | (see blk_oce_1 in sbcblk)
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          zzoce         = pCd_du    (ji,jj) * pwndm    (ji,jj)
 #if defined key_si3
          zzice         = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj)
@@ -137,39 +133,40 @@ CONTAINS
 #else
          ustar2(ji,jj) = zzoce
 #endif
-         !#LB: sorry Cdn_oce is gone:
-         !zrough(ji,jj) = ght_abl(2) * EXP( - vkarmn / SQRT( MAX( Cdn_oce(ji,jj), 1.e-4 ) ) ) !<-- recover the value of z0 from Cdn_oce
+         zsspt(ji,jj) = theta_exner( psst(ji,jj)+rt0, pslp_dta(ji,jj) )   ! potential SST [K]
       END_2D
 
+      !#LB: sorry Cdn_oce is gone:
+      !zrough(:,:) = ght_abl(2) * EXP( - vkarmn / SQRT( MAX( Cdn_oce(:,:), 1.e-4 ) ) ) !<-- recover the value of z0 from Cdn_oce
       zrough(:,:) = z0_from_Cd( ght_abl(2), pCd_du(:,:) / MAX( pwndm(:,:), 0.5_wp ) ) ! #LB: z0_from_Cd is define in sbc_phy.F90...
 
       ! sea surface potential temperature [K]
-      zsspt(:,:) = theta_exner( psst(:,:)+rt0, pslp_dta(:,:) )
+      !zsspt(:,:) = theta_exner( psst(A2D(0))+rt0, pslp_dta(:,:) )
+
 
-      !
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  1 *** Advance TKE to time n+1 and compute Avm_abl, Avt_abl, PBLh
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-      CALL abl_zdf_tke( )                       !--> Avm_abl, Avt_abl, pblh defined on (1,jpi) x (1,jpj)
+      CALL abl_zdf_tke( )                       !--> Avm_abl, Avt_abl, pblh defined on (A1Di(0)) x (A1Dj(0))
 
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  2 *** Advance tracers to time n+1
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
       !-------------
-      DO jj = 1, jpj    ! outer loop             !--> tq_abl computed on (1:jpi) x (1:jpj)
+      DO_1Dj(0,0)   ! outer loop                !--> tq_abl computed on (A1Di(0)) x (A1Dj(0))
       !-------------
          ! Compute matrix elements for interior points
          DO jk = 3, jpkam1
-            DO ji = 1, jpi   ! vector opt.
+            DO_1Di(0,0)
                z_elem_a( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
                z_elem_c( ji, jk ) = - rDt_abl * Avt_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal
                z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   !       diagonal
-            END DO
+            END_1D
          END DO
          ! Boundary conditions
-         DO ji = 1, jpi   ! vector opt.
+         DO_1Di(0,0)
             ! Neumann at the bottom
             z_elem_a( ji,    2 ) = 0._wp
             z_elem_c( ji,    2 ) = - rDt_abl * Avt_abl( ji, jj,    2 ) / e3w_abl(    2 )
@@ -177,19 +174,18 @@ CONTAINS
             z_elem_a( ji, jpka ) = - rDt_abl * Avt_abl( ji, jj, jpka ) / e3w_abl( jpka )
             z_elem_c( ji, jpka ) = 0._wp
             z_elem_b( ji, jpka ) = e3t_abl( jpka ) - z_elem_a( ji, jpka )
-         END DO
+         END_1D
 
          DO jtra = 1,jptq  ! loop on active tracers
 
             DO jk = 3, jpkam1
-               !DO ji = 2, jpim1
-               DO ji = 1,jpi  !!GS: to be checked if needed
+               DO_1Di(0,0)
                   tq_abl( ji, jj, jk, nt_a, jtra ) = e3t_abl(jk) * tq_abl( ji, jj, jk, nt_n, jtra )   ! initialize right-hand-side
-               END DO
+               END_1D
             END DO
 
             IF(jtra == jp_ta) THEN
-               DO ji = 1,jpi  ! surface boundary condition for temperature
+               DO_1Di(0,0)   ! surface boundary condition for temperature
                   zztmp1 = psen(ji, jj)
                   zztmp2 = psen(ji, jj) * zsspt(ji, jj)
 #if defined key_si3
@@ -199,9 +195,9 @@ CONTAINS
                   z_elem_b( ji,        2             ) = e3t_abl(    2 ) - z_elem_c( ji,        2             ) + rDt_abl * zztmp1
                   tq_abl  ( ji, jj,    2, nt_a, jtra ) = e3t_abl(    2 ) * tq_abl  ( ji, jj,    2, nt_n, jtra ) + rDt_abl * zztmp2
                   tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
-               END DO
+               END_1D
             ELSE ! jp_qa
-               DO ji = 1,jpi  ! surface boundary condition for humidity
+               DO_1Di(0,0)   ! surface boundary condition for humidity
                   zztmp1 = pevp(ji, jj)
                   zztmp2 = pevp(ji, jj) * pssq(ji, jj)
 #if defined key_si3
@@ -211,31 +207,31 @@ CONTAINS
                   z_elem_b( ji,     2                ) = e3t_abl(    2 ) - z_elem_c( ji,        2             ) + rDt_abl * zztmp1
                   tq_abl  ( ji, jj, 2   , nt_a, jtra ) = e3t_abl(    2 ) * tq_abl  ( ji, jj,    2, nt_n, jtra ) + rDt_abl * zztmp2
                   tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra )
-               END DO
+               END_1D
             END IF
             !!
             !! Matrix inversion
             !! ----------------------------------------------------------
-            DO ji = 1,jpi
+            DO_1Di(0,0)
                zcff                            =  1._wp / z_elem_b( ji, 2 )
                zCF   ( ji,     2             ) = - zcff * z_elem_c( ji, 2 )
                tq_abl( ji, jj, 2, nt_a, jtra ) =   zcff * tq_abl( ji, jj, 2, nt_a, jtra )
-            END DO
+            END_1D
 
             DO jk = 3, jpka
-               DO ji = 1,jpi
+               DO_1Di(0,0)
                   zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF( ji, jk-1 ) )
                   zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
                   tq_abl(ji,jj,jk,nt_a,jtra) = zcff * ( tq_abl(ji,jj,jk  ,nt_a,jtra)   &
                      &           - z_elem_a(ji, jk) *   tq_abl(ji,jj,jk-1,nt_a,jtra) )
-               END DO
+               END_1D
             END DO
             !!FL at this point we could check positivity of tq_abl(:,:,:,nt_a,jp_qa) ... test to do ...
             DO jk = jpkam1,2,-1
-               DO ji = 1,jpi
+               DO_1Di(0,0)
                   tq_abl(ji,jj,jk,nt_a,jtra) = tq_abl(ji,jj,jk,nt_a,jtra) +    &
                      &                        zCF(ji,jk) * tq_abl(ji,jj,jk+1,nt_a,jtra)
-               END DO
+               END_1D
             END DO
 
          END DO   !<-- loop on tracers
@@ -254,7 +250,7 @@ CONTAINS
          !-------------
             !
             ! Advance u_abl & v_abl to time n+1
-            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            DO_2D( 0, 0, 0, 0 )
                zcff = ( fft_abl(ji,jj) * rDt_abl )*( fft_abl(ji,jj) * rDt_abl )  ! (f dt)**2
 
                u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *(                                          &
@@ -269,13 +265,13 @@ CONTAINS
             END_2D
             !
          !-------------
-         END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
+         END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (A1Di(0)) x (A1Dj(0))
          !-------------
          !
          IF( ln_geos_winds ) THEN
-            DO jj = 1, jpj    ! outer loop
+            DO_1Dj( 0, 0 )   ! outer loop
                DO jk = 1, jpka
-                  DO ji = 1, jpi
+                  DO_1Di( 0, 0 )
                      u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a )   &
                         &                      - rDt_abl * e3t_abl(jk) * fft_abl(ji  , jj) * pgv_dta(ji  ,jj  ,jk)
                      v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a )   &
@@ -286,14 +282,14 @@ CONTAINS
          END IF
          !
          IF( ln_hpgls_frc ) THEN
-            DO jj = 1, jpj    ! outer loop
+            DO_1Dj( 0, 0 )    ! outer loop
                DO jk = 1, jpka
-                  DO ji = 1, jpi
+                  DO_1Di( 0, 0 )
                      u_abl( ji, jj, jk, nt_a ) = u_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgu_dta(ji,jj,jk)
                      v_abl( ji, jj, jk, nt_a ) = v_abl( ji, jj, jk, nt_a ) - rDt_abl * e3t_abl(jk) * pgv_dta(ji,jj,jk)
-                  ENDDO
+                  END_1D
                ENDDO
-            ENDDO
+            END_1D
          END IF
 
       ELSE ! SemiImp_Cor = .FALSE.
@@ -306,29 +302,25 @@ CONTAINS
                !
                IF( MOD( kt, 2 ) == 0 ) then
                   ! Advance u_abl & v_abl to time n+1
-                  DO jj = 1, jpj
-                     DO ji = 1, jpi
-                        zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_n ) - pgv_dta(ji  ,jj  ,jk)  )
-                        u_abl( ji, jj, jk, nt_a ) =                u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff
-                        zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_a ) - pgu_dta(ji  ,jj  ,jk)  )
-                        v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff )
-                        u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  u_abl( ji, jj, jk, nt_a )
-                     END DO
-                  END DO
+                  DO_2D( 0, 0, 0, 0 )
+                     zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_n ) - pgv_dta(ji  ,jj  ,jk)  )
+                     u_abl( ji, jj, jk, nt_a ) =                u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff
+                     zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_a ) - pgu_dta(ji  ,jj  ,jk)  )
+                     v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff )
+                     u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  u_abl( ji, jj, jk, nt_a )
+                  END_2D
                ELSE
-                  DO jj = 1, jpj
-                     DO ji = 1, jpi
-                        zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_n ) - pgu_dta(ji  ,jj  ,jk)  )
-                        v_abl( ji, jj, jk, nt_a ) =                v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff
-                        zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_a ) - pgv_dta(ji  ,jj  ,jk)  )
-                        u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff )
-                        v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  v_abl( ji, jj, jk, nt_a )
-                     END DO
-                  END DO
+                  DO_2D( 0, 0, 0, 0 )
+                     zcff = fft_abl(ji,jj) * ( u_abl ( ji , jj  , jk, nt_n ) - pgu_dta(ji  ,jj  ,jk)  )
+                     v_abl( ji, jj, jk, nt_a ) =                v_abl( ji, jj, jk, nt_n ) - rDt_abl * zcff
+                     zcff = fft_abl(ji,jj) * ( v_abl ( ji , jj  , jk, nt_a ) - pgv_dta(ji  ,jj  ,jk)  )
+                     u_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *( u_abl( ji, jj, jk, nt_n ) + rDt_abl * zcff )
+                     v_abl( ji, jj, jk, nt_a ) = e3t_abl(jk) *  v_abl( ji, jj, jk, nt_a )
+                  END_2D
                END IF
                !
             !-------------
-            END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (1:jpi) x (1:jpj)
+            END DO             ! end outer loop  !<-- u_abl and v_abl are properly updated on (A1Di(0)) x (A1Dj(0))
             !-------------
 
          ENDIF ! ln_geos_winds
@@ -340,18 +332,19 @@ CONTAINS
       !
       !  Vertical diffusion for u_abl
       !-------------
-      DO jj = 1, jpj    ! outer loop
+      DO_1Dj(0,0)   ! outer loop
       !-------------
 
          DO jk = 3, jpkam1
-            DO ji = 1, jpi
+            DO_1Di(0,0)
                z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
                z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal
                z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   !       diagonal
-            END DO
+            END_1D
          END DO
 
-         DO ji = 2, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at ji=jpi)
+         DO_1Di(0,0)   ! boundary conditions
+
             !++ Surface boundary condition
             z_elem_a( ji, 2 ) = 0._wp
             z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
@@ -381,29 +374,29 @@ CONTAINS
                u_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pu_dta(ji,jj,jk)
             !ENDIF
 
-         END DO
+         END_1D
          !!
          !! Matrix inversion
          !! ----------------------------------------------------------
-         DO ji = 1, jpi  !DO ji = 2, jpi !!GS: TBI
+         DO_1Di(0,0)
             zcff                 =   1._wp / z_elem_b( ji, 2 )
             zCF   (ji,   2     ) =  - zcff * z_elem_c( ji,     2       )
             u_abl (ji,jj,2,nt_a) =    zcff * u_abl   ( ji, jj, 2, nt_a )
-         END DO
+         END_1D
 
          DO jk = 3, jpka
-            DO ji = 2, jpi
+            DO_1Di(0,0)
                zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
                zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
                u_abl(ji,jj,jk,nt_a) = zcff * ( u_abl(ji,jj,jk  ,nt_a)   &
                &          - z_elem_a(ji, jk) * u_abl(ji,jj,jk-1,nt_a) )
-            END DO
+            END_1D
          END DO
 
          DO jk = jpkam1,2,-1
-            DO ji = 2, jpi
+            DO_1Di(0,0)
                u_abl(ji,jj,jk,nt_a) = u_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * u_abl(ji,jj,jk+1,nt_a)
-            END DO
+            END_1D
          END DO
 
       !-------------
@@ -413,18 +406,19 @@ CONTAINS
       !
       !  Vertical diffusion for v_abl
       !-------------
-      DO jj = 2, jpj   ! outer loop
+      DO_1Dj(0,0)   ! outer loop
       !-------------
          !
          DO jk = 3, jpkam1
-            DO ji = 1, jpi
+            DO_1Di(0,0)
                z_elem_a( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk-1 ) / e3w_abl( jk-1 )   ! lower-diagonal
                z_elem_c( ji, jk ) = - rDt_abl * Avm_abl( ji, jj, jk   ) / e3w_abl( jk   )   ! upper-diagonal
                z_elem_b( ji, jk ) = e3t_abl(jk) - z_elem_a( ji, jk ) - z_elem_c( ji, jk )   !       diagonal
-            END DO
+            END_1D
          END DO
 
-         DO ji = 1, jpi   ! boundary conditions (Avm_abl and pcd_du must be available at jj=jpj)
+         DO_1Di(0,0)   ! boundary conditions
+
             !++ Surface boundary condition
             z_elem_a( ji, 2 ) = 0._wp
             z_elem_c( ji, 2 ) = - rDt_abl * Avm_abl( ji, jj, 2 ) / e3w_abl( 2 )
@@ -454,29 +448,29 @@ CONTAINS
                v_abl   ( ji, jj, jpka, nt_a ) = e3t_abl( jpka ) * pv_dta(ji,jj,jk)
             !ENDIF
 
-         END DO
+         END_1D
          !!
          !! Matrix inversion
          !! ----------------------------------------------------------
-         DO ji = 1, jpi
+         DO_1Di(0,0)
             zcff                 =   1._wp / z_elem_b( ji, 2 )
             zCF   (ji,   2     ) =  - zcff * z_elem_c( ji,     2       )
             v_abl (ji,jj,2,nt_a) =    zcff * v_abl   ( ji, jj, 2, nt_a )
-         END DO
+         END_1D
 
          DO jk = 3, jpka
-            DO ji = 1, jpi
+            DO_1Di(0,0)
                zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF   (ji, jk-1 ) )
                zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
                v_abl(ji,jj,jk,nt_a) = zcff * ( v_abl(ji,jj,jk  ,nt_a)   &
                &          - z_elem_a(ji, jk) * v_abl(ji,jj,jk-1,nt_a) )
-            END DO
+            END_1D
          END DO
 
          DO jk = jpkam1,2,-1
-            DO ji = 1, jpi
+            DO_1Di(0,0)
                v_abl(ji,jj,jk,nt_a) = v_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * v_abl(ji,jj,jk+1,nt_a)
-            END DO
+            END_1D
          END DO
          !
       !-------------
@@ -491,7 +485,7 @@ CONTAINS
          !-------------
          DO jk = 2, jpka    ! outer loop
          !-------------
-            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+            DO_2D( 0, 0, 0, 0)
                zcff1 = pblh( ji, jj )
                zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )
                zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )
@@ -514,7 +508,7 @@ CONTAINS
       !-------------
       DO jk = 2, jpka    ! outer loop
       !-------------
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0)
             zcff1 = pblh( ji, jj )
             zsig  = ght_abl(jk) / MAX( jp_pblh_min,  MIN(  jp_pblh_max, zcff1  ) )
             zsig  =               MIN( jp_bmax    ,  MAX(         zsig, jp_bmin) )
@@ -528,21 +522,17 @@ CONTAINS
 
             tq_abl( ji, jj, jk, nt_a, jp_qa ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_qa )   &
                &                                       + zcff   * pq_dta( ji, jj, jk )
-
          END_2D
       !-------------
       END DO             ! end outer loop
       !-------------
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                            !  6 *** MPI exchanges & IOM outputs
+      !                            !  6 *** IOM outputs
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      CALL lbc_lnk( 'ablmod',  u_abl(:,:,:,nt_a      ), 'T', -1._wp,  v_abl(:,:,:,nt_a)      , 'T', -1._wp                            )
-      !CALL lbc_lnk( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T',  1._wp, tq_abl(:,:,:,nt_a,jp_qa), 'T',  1._wp, kfillmode = jpfillnothing )   ! ++++ this should not be needed...
-      !
 #if defined key_xios
       ! 2D & first ABL level
-      IF ( iom_use("pblh"   ) ) CALL iom_put (    "pblh",    pblh(:,:             ) )
+      IF ( iom_use("pblh"   ) ) CALL iom_put (    "pblh",    pblh(A2D(0)          ) )
       IF ( iom_use("uz1_abl") ) CALL iom_put ( "uz1_abl",   u_abl(:,:,2,nt_a      ) )
       IF ( iom_use("vz1_abl") ) CALL iom_put ( "vz1_abl",   v_abl(:,:,2,nt_a      ) )
       IF ( iom_use("tz1_abl") ) CALL iom_put ( "tz1_abl",  tq_abl(:,:,2,nt_a,jp_ta) )
@@ -588,7 +578,7 @@ CONTAINS
       !                            !  7 *** Finalize flux computation
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          ztemp          = tq_abl( ji, jj, 2, nt_a, jp_ta )
          zhumi          = tq_abl( ji, jj, 2, nt_a, jp_qa )
          zpre( ji, jj ) = pres_temp( zhumi, pslp_dta(ji,jj), ght_abl(2), ptpot=ztemp, pta=ztabs( ji, jj ) )
@@ -599,15 +589,13 @@ CONTAINS
          rhoa( ji, jj ) = zcff
       END_2D
 
-      DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          zwnd_i(ji,jj) = u_abl(ji  ,jj,2,nt_a) - 0.5_wp * ( pssu(ji  ,jj) + pssu(ji-1,jj) ) * rn_vfac
          zwnd_j(ji,jj) = v_abl(ji,jj  ,2,nt_a) - 0.5_wp * ( pssv(ji,jj  ) + pssv(ji,jj-1) ) * rn_vfac
       END_2D
       !
-      CALL lbc_lnk( 'ablmod', zwnd_i(:,:) , 'T', -1.0_wp, zwnd_j(:,:) , 'T', -1.0_wp )
-      !
       ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          zcff          = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   &
             &                 + zwnd_j(ji,jj) * zwnd_j(ji,jj) )   ! * msk_abl(ji,jj)
          zztmp         = rhoa(ji,jj) * pcd_du(ji,jj)
@@ -617,13 +605,14 @@ CONTAINS
          zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
          zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
       END_2D
+      CALL iom_put( "taum_oce", ptaum )
+
       ! ... utau, vtau at T-points
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         ptaui(ji,jj) = zwnd_i(ji,jj) * msk_abl(ji,jj) !!clem tau: check
-         ptauj(ji,jj) = zwnd_j(ji,jj) * msk_abl(ji,jj) !!clem tau: check
+      DO_2D( 0, 0, 0, 0 )
+         ptaui(ji,jj) = zwnd_i(ji,jj) * msk_abl(ji,jj)
+         ptauj(ji,jj) = zwnd_j(ji,jj) * msk_abl(ji,jj)
       END_2D
-
-      CALL iom_put( "taum_oce", ptaum )
+      CALL lbc_lnk( 'ablmod', ptaui(:,:) , 'T', -1.0_wp, ptauj(:,:) , 'T', -1.0_wp )
 
       IF(sn_cfctl%l_prtctl) THEN
          CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau   : ', mask1=tmask,   &
@@ -635,24 +624,12 @@ CONTAINS
       ! ------------------------------------------------------------ !
       !    Wind stress relative to the moving ice ( U10m - U_ice )   !
       ! ------------------------------------------------------------ !
-      !DO_2D( 0, 0, 0, 0 )
-      !   ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   &
-      !      &                      * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) )
-      !   ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   &
-      !      &                      * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) )
-      !END_2D
-      !CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp )
-      !!
-      !IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   &
-      !   &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' )
-
-      ! ------------------------------------------------------------ !
-      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
-      ! ------------------------------------------------------------ !
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         ptaui_ice(ji,jj) = rhoa(ji,jj) * pCd_du_ice(ji,jj) * ( u_abl(ji,jj,2,nt_a) - pssu_ice(ji,jj) * rn_vfac )
-         ptauj_ice(ji,jj) = rhoa(ji,jj) * pCd_du_ice(ji,jj) * ( v_abl(ji,jj,2,nt_a) - pssv_ice(ji,jj) * rn_vfac )
+      DO_2D(0,0,0,0)
+         zztmp            = rhoa(ji,jj) * pCd_du_ice(ji,jj)
+         ptaui_ice(ji,jj) = zztmp * ( u_abl(ji,jj,2,nt_a) - 0.5_wp * ( pssu_ice(ji,jj) + pssu_ice(ji-1,jj  ) ) * rn_vfac )
+         ptauj_ice(ji,jj) = zztmp * ( v_abl(ji,jj,2,nt_a) - 0.5_wp * ( pssv_ice(ji,jj) + pssv_ice(ji  ,jj-1) ) * rn_vfac )
       END_2D
+      CALL lbc_lnk( 'ablmod', ptaui_ice(:,:) , 'T', -1.0_wp, ptauj_ice(:,:) , 'T', -1.0_wp )
       !
       IF(sn_cfctl%l_prtctl) THEN
          CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=tmask,   &
@@ -690,21 +667,22 @@ CONTAINS
       !!              - avmu, avmv : production of TKE by shear at u and v-points
       !!                (= Kz dz[Ub] * dz[Un] )
       !! ---------------------------------------------------------------------
-      INTEGER                                 ::   ji, jj, jk, tind, jbak, jkup, jkdwn
-      INTEGER, DIMENSION(1:jpi          )     ::   ikbl
-      REAL(wp)                                ::   zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
-      REAL(wp)                                ::   zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2
-      REAL(wp)                                ::   zdU2, zdV2, zbuoy1, zbuoy2    ! zbuoy for BL89
-      REAL(wp)                                ::   zwndi, zwndj
-      REAL(wp), DIMENSION(1:jpi,      1:jpka) ::   zsh2
-      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) ::   zbn2
-      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   zFC, zRH, zCF
-      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_a
-      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_b
-      REAL(wp), DIMENSION(1:jpi,1:jpka  )     ::   z_elem_c
-      LOGICAL                                 ::   ln_Patankar    = .FALSE.
-      LOGICAL                                 ::   ln_dumpvar     = .FALSE.
-      LOGICAL , DIMENSION(1:jpi         )     ::   ln_foundl
+
+      INTEGER                             ::   ji, jj, jk, tind, jbak, jkup, jkdwn
+      INTEGER, DIMENSION(A1Di(0))         ::   ikbl
+      REAL(wp)                            ::   zcff, zcff2, ztken, zesrf, zetop, ziRic, ztv
+      REAL(wp)                            ::   zdU , zdV , zcff1, zshear, zbuoy, zsig, zustar2
+      REAL(wp)                            ::   zdU2, zdV2, zbuoy1, zbuoy2    ! zbuoy for BL89
+      REAL(wp)                            ::   zwndi, zwndj
+      REAL(wp), DIMENSION(A1Di(0),1:jpka) ::   zsh2
+      REAL(wp), DIMENSION(A2D(0) ,1:jpka) ::   zbn2
+      REAL(wp), DIMENSION(A1Di(0),1:jpka) ::   zFC, zRH, zCF
+      REAL(wp), DIMENSION(A1Di(0),1:jpka) ::   z_elem_a
+      REAL(wp), DIMENSION(A1Di(0),1:jpka) ::   z_elem_b
+      REAL(wp), DIMENSION(A1Di(0),1:jpka) ::   z_elem_c
+      LOGICAL                             ::   ln_Patankar    = .FALSE.
+      LOGICAL                             ::   ln_dumpvar     = .FALSE.
+      LOGICAL , DIMENSION(A1Di(0))        ::   ln_foundl
       !
       tind  = nt_n
       ziRic = 1._wp / rn_Ric
@@ -713,33 +691,33 @@ CONTAINS
       !                            !  Advance TKE equation to time n+1
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !-------------
-      DO jj = 1, jpj    ! outer loop
+      DO_1Dj( 0, 0 )   ! outer loop
       !-------------
          !
          ! Compute vertical shear
          DO jk = 2, jpkam1
-            DO ji = 1, jpi
+            DO_1Di( 0, 0 )
                zcff        = 1.0_wp / e3w_abl( jk )**2
                zdU         = zcff* Avm_abl(ji,jj,jk) * (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2
                zdV         = zcff* Avm_abl(ji,jj,jk) * (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
                zsh2(ji,jk) = zdU+zdV   !<-- zsh2 = Km ( ( du/dz )^2 + ( dv/dz )^2 )
-            END DO
+            END_1D
          END DO
          !
          ! Compute brunt-vaisala frequency
          DO jk = 2, jpkam1
-            DO ji = 1,jpi
+            DO_1Di( 0, 0 )
                zcff  = grav * itvref / e3w_abl( jk )
                zcff1 =  tq_abl( ji, jj, jk+1, tind, jp_ta) - tq_abl( ji, jj, jk  , tind, jp_ta)
                zcff2 =  tq_abl( ji, jj, jk+1, tind, jp_ta) * tq_abl( ji, jj, jk+1, tind, jp_qa)        &
                   &   - tq_abl( ji, jj, jk  , tind, jp_ta) * tq_abl( ji, jj, jk  , tind, jp_qa)
                zbn2(ji,jj,jk) = zcff * ( zcff1 + rctv0 * zcff2 )  !<-- zbn2 defined on (2,jpi)
-            END DO
+            END_1D
          END DO
          !
          ! Terms for the tridiagonal problem
          DO jk = 2, jpkam1
-            DO ji = 1, jpi
+            DO_1Di( 0, 0 )
                zshear      = zsh2( ji, jk )                           ! zsh2 is already multiplied by Avm_abl at this point
                zsh2(ji,jk) = zsh2( ji, jk ) / Avm_abl( ji, jj, jk )   ! reformulate zsh2 as a 'true' vertical shear for PBLH computation
                zbuoy       = - Avt_abl( ji, jj, jk ) * zbn2( ji, jj, jk )
@@ -756,10 +734,10 @@ CONTAINS
                      &               - e3w_abl(jk) * rDt_abl * zbuoy
                   tke_abl( ji, jj, jk, nt_a ) = e3w_abl(jk) * ( tke_abl( ji, jj, jk, nt_n ) + rDt_abl *  zshear )                    ! right-hand-side
                END IF
-            END DO
+            END_1D
          END DO
 
-         DO ji = 1,jpi    ! vector opt.
+         DO_1Di( 0, 0 )
             zesrf = MAX( rn_Esfc * ustar2(ji,jj), tke_min )
             zetop = tke_min
 
@@ -788,44 +766,41 @@ CONTAINS
          !!
          !! Matrix inversion
          !! ----------------------------------------------------------
-         DO ji = 1,jpi
+         DO_1Di( 0, 0 )
             zcff                  =  1._wp / z_elem_b( ji, 1 )
             zCF    (ji,   1     ) = - zcff * z_elem_c( ji,     1       )
             tke_abl(ji,jj,1,nt_a) =   zcff * tke_abl ( ji, jj, 1, nt_a )
-         END DO
+         END_1D
 
          DO jk = 2, jpka
-            DO ji = 1,jpi
+            DO_1Di( 0, 0 )
                zcff = 1._wp / ( z_elem_b( ji, jk ) + z_elem_a( ji, jk ) * zCF(ji, jk-1 ) )
                zCF(ji,jk) = - zcff * z_elem_c( ji, jk )
                tke_abl(ji,jj,jk,nt_a) =   zcff * ( tke_abl(ji,jj,jk  ,nt_a)   &
                &          - z_elem_a(ji, jk) * tke_abl(ji,jj,jk-1,nt_a) )
-            END DO
+            END_1D
          END DO
 
          DO jk = jpkam1,1,-1
-            DO ji = 1,jpi
+            DO_1Di( 0, 0 )
                tke_abl(ji,jj,jk,nt_a) = tke_abl(ji,jj,jk,nt_a) + zCF(ji,jk) * tke_abl(ji,jj,jk+1,nt_a)
-            END DO
+            END_1D
          END DO
 
-!!FL should not be needed because of Patankar procedure
-         tke_abl(2:jpi,jj,1:jpka,nt_a) = MAX( tke_abl(2:jpi,jj,1:jpka,nt_a), tke_min )
+         !!FL should not be needed because of Patankar procedure
+         tke_abl(A2D(0),1:jpka,nt_a) = MAX( tke_abl(A2D(0),1:jpka,nt_a), tke_min )
 
          !!
          !! Diagnose PBL height
          !! ----------------------------------------------------------
-
-
-         !
          ! arrays zRH, zFC and zCF are available at this point
          ! and zFC(:, 1 ) = 0.
          ! diagnose PBL height based on zsh2 and zbn2
          zFC (  :  ,1) = 0._wp
-         ikbl( 1:jpi ) = 0
+         ikbl( A1Di(0) ) = 0
 
          DO jk = 2,jpka
-            DO ji = 1, jpi
+            DO_1Di( 0, 0 )
                zcff  = ghw_abl( jk-1 )
                zcff1 = zcff / ( zcff + rn_epssfc * pblh ( ji, jj ) )
                zcff  = ghw_abl( jk   )
@@ -837,11 +812,11 @@ CONTAINS
                            - rn_Cek  * ( fft_abl( ji, jj  ) * fft_abl( ji, jj ) ) ) &
                            &                                                 )
                IF( ikbl(ji) == 0 .and. zFC( ji, jk ).lt.0._wp ) ikbl(ji)=jk
-            END DO
+            END_1D
          END DO
          !
          ! finalize the computation of the PBL height
-         DO ji = 1, jpi
+         DO_1Di( 0, 0 )
             jk = ikbl(ji)
             IF( jk > 2 ) THEN ! linear interpolation to get subgrid value of pblh
                pblh( ji, jj ) =  (  ghw_abl( jk-1 ) * zFC( ji, jk   )       &
@@ -852,16 +827,16 @@ CONTAINS
             ELSE
                pblh( ji, jj ) = ghw_abl(jpka)
             END IF
-         END DO
+         END_1D
       !-------------
       END DO
       !-------------
       !
       ! Optional : could add pblh smoothing if pblh is noisy horizontally ...
       IF(ln_smth_pblh) THEN
-         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing)
-         CALL smooth_pblh( pblh, msk_abl )
-         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp) !, kfillmode = jpfillnothing)
+         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp )
+         CALL smooth_pblh( pblh, tmask(A2D(0),1) )
+         CALL lbc_lnk( 'ablmod', pblh, 'T', 1.0_wp )
       ENDIF
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  Diagnostic mixing length computation
@@ -872,55 +847,55 @@ CONTAINS
       CASE ( 0 )           ! Deardroff 80 length-scale bounded by the distance to surface and bottom
 #   define zlup zRH
 #   define zldw zFC
-         DO jj = 1, jpj     ! outer loop
+         DO_1Dj(0,0)   ! outer loop
             !
-            DO ji = 1, jpi
+            DO_1Di(0,0)
                mxld_abl( ji, jj,    1 ) = mxl_min
                mxld_abl( ji, jj, jpka ) = mxl_min
                mxlm_abl( ji, jj,    1 ) = mxl_min
                mxlm_abl( ji, jj, jpka ) = mxl_min
                zldw    ( ji,        1 ) = zrough(ji,jj) * rn_Lsfc
                zlup    ( ji,     jpka ) = mxl_min
-            END DO
+            END_1D
             !
             DO jk = 2, jpkam1
-               DO ji = 1, jpi
+               DO_1Di(0,0)
                   zbuoy = MAX( zbn2(ji, jj, jk), rsmall )
                   mxlm_abl( ji, jj, jk ) = MAX( mxl_min,  &
                      &               SQRT( 2._wp * tke_abl( ji, jj, jk, nt_a ) / zbuoy ) )
-               END DO
+               END_1D
             END DO
             !
             ! Limit mxl
             DO jk = jpkam1,1,-1
-               DO ji = 1, jpi
+               DO_1Di(0,0)
                   zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
-               END DO
+               END_1D
             END DO
             !
             DO jk = 2, jpka
-               DO ji = 1, jpi
+               DO_1Di(0,0)
                   zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
-               END DO
+               END_1D
             END DO
             !
 !            DO jk = 1, jpka
-!               DO ji = 1, jpi
+!               DO_1Di(0,0)
 !                  mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
 !                  mxld_abl( ji, jj, jk ) = MIN ( zldw( ji, jk ),  zlup( ji, jk ) )
-!               END DO
+!               END_1D
 !            END DO
             !
             DO jk = 1, jpka
-               DO ji = 1, jpi
+               DO_1Di(0,0)
 !                  zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
                   zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
                   mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
                   mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
-               END DO
+               END_1D
             END DO
             !
-         END DO
+         END_1D   ! outer loop
 #   undef zlup
 #   undef zldw
          !
@@ -928,18 +903,18 @@ CONTAINS
       CASE ( 1 )           ! Modified Deardroff 80 length-scale bounded by the distance to surface and bottom
 #   define zlup zRH
 #   define zldw zFC
-         DO jj = 1, jpj     ! outer loop
+         DO_1Dj( 0, 0 )   ! outer loop
             !
             DO jk = 2, jpkam1
-               DO ji = 1,jpi
-                              zcff        = 1.0_wp / e3w_abl( jk )**2
+               DO_1Di( 0, 0 )
+                  zcff        = 1.0_wp / e3w_abl( jk )**2
                   zdU         = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2
                   zdV         = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
                   zsh2(ji,jk) = SQRT(zdU+zdV)   !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
-                           ENDDO
-                        ENDDO
-                        !
-            DO ji = 1, jpi
+               END_1D
+            END DO
+            !
+            DO_1Di( 0, 0 )
                zcff                      = zrough(ji,jj) * rn_Lsfc
                mxld_abl ( ji, jj,    1 ) = zcff
                mxld_abl ( ji, jj, jpka ) = mxl_min
@@ -947,42 +922,42 @@ CONTAINS
                mxlm_abl ( ji, jj, jpka ) = mxl_min
                zldw     ( ji,        1 ) = zcff
                zlup     ( ji,     jpka ) = mxl_min
-            END DO
+            END_1D
             !
             DO jk = 2, jpkam1
-               DO ji = 1, jpi
+               DO_1Di( 0, 0 )
                   zbuoy    = MAX( zbn2(ji, jj, jk), rsmall )
                   zcff     = 2.0_wp*SQRT(tke_abl( ji, jj, jk, nt_a )) / ( rn_Rod*zsh2(ji,jk) &
                                 &             + SQRT(rn_Rod*rn_Rod*zsh2(ji,jk)*zsh2(ji,jk)+2.0_wp*zbuoy ) )
                                   mxlm_abl( ji, jj, jk ) = MAX( mxl_min, zcff )
-               END DO
+               END_1D
             END DO
             !
             ! Limit mxl
             DO jk = jpkam1,1,-1
-               DO ji = 1, jpi
+               DO_1Di( 0, 0 )
                   zlup(ji,jk) = MIN( zlup(ji,jk+1) + (ghw_abl(jk+1)-ghw_abl(jk)) , mxlm_abl(ji, jj, jk) )
-               END DO
+               END_1D
             END DO
             !
             DO jk = 2, jpka
-               DO ji = 1, jpi
+               DO_1Di( 0, 0 )
                   zldw(ji,jk) = MIN( zldw(ji,jk-1) + (ghw_abl(jk)-ghw_abl(jk-1)) , mxlm_abl(ji, jj, jk) )
-               END DO
+               END_1D
             END DO
             !
             DO jk = 1, jpka
-               DO ji = 1, jpi
+               DO_1Di( 0, 0 )
                   !mxlm_abl( ji, jj, jk ) = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
                   !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
                   zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
                   mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
                   !mxld_abl( ji, jj, jk ) = MIN( zldw( ji, jk ), zlup( ji, jk ) )
                   mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
-               END DO
+               END_1D
             END DO
- !
-         END DO
+            !
+         END_1D   ! outer loop
 #   undef zlup
 #   undef zldw
          !
@@ -992,35 +967,35 @@ CONTAINS
 #   define zldw zFC
 ! zCF is used for matrix inversion
 !
-       DO jj = 1, jpj      ! outer loop
+       DO_1Dj( 0, 0 )   ! outer loop
 
-         DO ji = 1, jpi
+         DO_1Di( 0, 0 )
             zcff             = zrough(ji,jj) * rn_Lsfc
             zlup( ji,    1 ) = zcff
             zldw( ji,    1 ) = zcff
             zlup( ji, jpka ) = mxl_min
             zldw( ji, jpka ) = mxl_min
-         END DO
+         END_1D
 
          DO jk = 2,jpka-1
-            DO ji = 1, jpi
+            DO_1Di( 0, 0 )
                zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
                zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)
-            END DO
+            END_1D
          END DO
          !!
          !! BL89 search for lup
          !! ----------------------------------------------------------
          DO jk=2,jpka-1
             !
-            DO ji = 1, jpi
+            DO_1Di( 0, 0 )
                zCF(ji,1:jpka) = 0._wp
                zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
                ln_foundl(ji ) = .false.
-            END DO
+            END_1D
             !
             DO jkup=jk+1,jpka-1
-               DO ji = 1, jpi
+               DO_1Di( 0, 0 )
                   zbuoy1 = MAX( zbn2(ji,jj,jkup  ), rsmall )
                   zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall )
                   zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
@@ -1035,7 +1010,7 @@ CONTAINS
                      zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk)
                      ln_foundl(ji) = .true.
                   END IF
-               END DO
+               END_1D
             END DO
             !
          END DO
@@ -1044,14 +1019,14 @@ CONTAINS
          !! ----------------------------------------------------------
          DO jk=2,jpka-1
             !
-            DO ji = 1, jpi
+            DO_1Di( 0, 0 )
                zCF(ji,1:jpka) = 0._wp
                zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
                ln_foundl(ji ) = .false.
-            END DO
+            END_1D
             !
             DO jkdwn=jk-1,1,-1
-               DO ji = 1, jpi
+               DO_1Di( 0, 0 )
                   zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
                   zbuoy2 = MAX( zbn2(ji,jj,jkdwn  ), rsmall )
                   zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
@@ -1066,21 +1041,21 @@ CONTAINS
                      zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  )
                      ln_foundl(ji) = .true.
                   END IF
-               END DO
+               END_1D
             END DO
             !
          END DO
 
          DO jk = 1, jpka
-            DO ji = 1, jpi
+            DO_1Di( 0, 0 )
                !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
                zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
                mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
                mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
-            END DO
+            END_1D
          END DO
 
-      END DO
+      END_1D   ! outer loop
 #   undef zlup
 #   undef zldw
          !
@@ -1090,46 +1065,46 @@ CONTAINS
 #   define zldw zFC
 ! zCF is used for matrix inversion
 !
-       DO jj = 1, jpj      ! outer loop
+       DO_1Dj( 0, 0 )   ! outer loop
           !
           DO jk = 2, jpkam1
-             DO ji = 1,jpi
+             DO_1Di( 0, 0 )
                             zcff        = 1.0_wp / e3w_abl( jk )**2
                 zdU         = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2
                 zdV         = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
                 zsh2(ji,jk) = SQRT(zdU+zdV)   !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
-                         ENDDO
-                  ENDDO
+             END_1D
+          ENDDO
           zsh2(:,      1) = zsh2( :,      2)
           zsh2(:,   jpka) = zsh2( :, jpkam1)
 
-                 DO ji = 1, jpi
-            zcff              = zrough(ji,jj) * rn_Lsfc
-                        zlup( ji,    1 )  = zcff
-            zldw( ji,    1 )  = zcff
-            zlup( ji, jpka ) = mxl_min
-            zldw( ji, jpka ) = mxl_min
-         END DO
+          DO_1Di( 0, 0 )
+             zcff              = zrough(ji,jj) * rn_Lsfc
+             zlup( ji,    1 )  = zcff
+             zldw( ji,    1 )  = zcff
+             zlup( ji, jpka ) = mxl_min
+             zldw( ji, jpka ) = mxl_min
+          END_1D
 
          DO jk = 2,jpka-1
-            DO ji = 1, jpi
+            DO_1Di( 0, 0 )
                zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
                zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)
-            END DO
+            END_1D
          END DO
          !!
          !! BL89 search for lup
          !! ----------------------------------------------------------
          DO jk=2,jpka-1
             !
-            DO ji = 1, jpi
+            DO_1Di( 0, 0 )
                zCF(ji,1:jpka) = 0._wp
                zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
                ln_foundl(ji ) = .false.
-            END DO
+            END_1D
             !
             DO jkup=jk+1,jpka-1
-               DO ji = 1, jpi
+               DO_1Di( 0, 0 )
                   zbuoy1             = MAX( zbn2(ji,jj,jkup  ), rsmall )
                   zbuoy2             = MAX( zbn2(ji,jj,jkup-1), rsmall )
                   zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) *                 &
@@ -1148,7 +1123,7 @@ CONTAINS
                      zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk)
                      ln_foundl(ji) = .true.
                   END IF
-               END DO
+               END_1D
             END DO
             !
          END DO
@@ -1157,14 +1132,14 @@ CONTAINS
          !! ----------------------------------------------------------
          DO jk=2,jpka-1
             !
-            DO ji = 1, jpi
+            DO_1Di( 0, 0 )
                zCF(ji,1:jpka) = 0._wp
                zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
                ln_foundl(ji ) = .false.
-            END DO
+            END_1D
             !
             DO jkdwn=jk-1,1,-1
-               DO ji = 1, jpi
+               DO_1Di( 0, 0 )
                   zbuoy1             = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
                   zbuoy2             = MAX( zbn2(ji,jj,jkdwn  ), rsmall )
                   zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
@@ -1183,18 +1158,18 @@ CONTAINS
                      zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  )
                      ln_foundl(ji) = .true.
                   END IF
-               END DO
+               END_1D
             END DO
             !
          END DO
 
          DO jk = 1, jpka
-            DO ji = 1, jpi
+            DO_1Di( 0, 0 )
                !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
                zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
                mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
                mxld_abl( ji, jj, jk ) = MAX(  MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
-            END DO
+            END_1D
          END DO
 
       END DO
@@ -1208,10 +1183,10 @@ CONTAINS
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
       !-------------
-      DO jj = 1, jpj     ! outer loop
+      DO_1Dj( 0, 0 )   ! outer loop
       !-------------
          DO jk = 1, jpka
-            DO ji = 1, jpi  ! vector opt.
+            DO_1Di( 0, 0 )
                zcff  = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk )  &
                &     * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) )
                zcff2 =  1. / ( 1. + zcff )   !<-- phi_z(z)
@@ -1219,10 +1194,10 @@ CONTAINS
                !!FL: MAX function probably useless because of the definition of mxl_min
                Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff         , avm_bak   )
                Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   )
-            END DO
+            END_1D
          END DO
       !-------------
-      END DO
+      END_1D   ! outer loop
       !-------------
 
 !---------------------------------------------------------------------------------------------------
@@ -1240,42 +1215,43 @@ CONTAINS
       !! ** Purpose :   2D Hanning filter on atmospheric PBL height
       !!
       !! ---------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: msk
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
-      INTEGER                                     :: ji,jj
-      REAL(wp)                                    :: smth_a, smth_b
-      REAL(wp), DIMENSION(jpi,jpj)                :: zdX,zdY,zFX,zFY
-      REAL(wp)                                    :: zumsk,zvmsk
-      !!
+
+      REAL(wp), DIMENSION(A2D(1)), INTENT(in   ) :: msk
+      REAL(wp), DIMENSION(A2D(1)), INTENT(inout) :: pvar2d
+      INTEGER                                    :: ji,jj
+      REAL(wp)                                   :: smth_a, smth_b
+      REAL(wp), DIMENSION(A2D(1))                :: zdX,zdY,zFX,zFY
+      REAL(wp)                                   :: zumsk,zvmsk
+
       !!=========================================================
       !!
       !! Hanning filter
       smth_a = 1._wp / 8._wp
       smth_b = 1._wp / 4._wp
       !
-      DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls )
+      DO_2D( 1, 0, 1, 1 )
          zumsk = msk(ji,jj) * msk(ji+1,jj)
          zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji  ,jj ) ) * zumsk
       END_2D
 
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 )
+      DO_2D( 1, 1, 1, 0 )
          zvmsk = msk(ji,jj) * msk(ji,jj+1)
          zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji  ,jj ) ) * zvmsk
       END_2D
 
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )
+      DO_2D( 0, 0, 1, 0 )
          zFY ( ji, jj  ) =   zdY ( ji, jj   )                        &
             & +  smth_a*  ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 ))   &
             &            -  (zdX ( ji, jj   ) - zdX( ji-1, jj   ))  )
       END_2D
 
-      DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 1, 0, 0, 0 )
          zFX( ji, jj  ) =    zdX( ji, jj   )                         &
            &    + smth_a*(  (zdY( ji+1, jj ) - zdY( ji+1, jj-1))     &
            &             -  (zdY( ji  , jj ) - zdY( ji  , jj-1)) )
       END_2D
 
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D( 0, 0, 0, 0 )
          pvar2d( ji  ,jj ) = pvar2d( ji  ,jj )              &
   &         + msk(ji,jj) * smth_b * (                       &
   &                  zFX( ji, jj ) - zFX( ji-1, jj )        &
diff --git a/src/ABL/ablrst.F90 b/src/ABL/ablrst.F90
index 925234b3e73663d1aeeba26e750649f20f3ab062..de4172974cef6032d7050d3aca32e655e31caf07 100644
--- a/src/ABL/ablrst.F90
+++ b/src/ABL/ablrst.F90
@@ -202,12 +202,12 @@ CONTAINS
       ! --- mandatory fields --- ! 
       CALL iom_get( numrar, jpdom_auto,   'u_abl',   u_abl(:,:,:,nt_n      ), cd_type = 'T', psgn = -1._wp )
       CALL iom_get( numrar, jpdom_auto,   'v_abl',   v_abl(:,:,:,nt_n      ), cd_type = 'T', psgn = -1._wp )
-      CALL iom_get( numrar, jpdom_auto,   't_abl',  tq_abl(:,:,:,nt_n,jp_ta), kfill = jpfillcopy )
+      CALL iom_get( numrar, jpdom_auto,   't_abl',  tq_abl(:,:,:,nt_n,jp_ta) ) !, kfill = jpfillcopy )
       CALL iom_get( numrar, jpdom_auto,   'q_abl',  tq_abl(:,:,:,nt_n,jp_qa) )
       CALL iom_get( numrar, jpdom_auto, 'tke_abl', tke_abl(:,:,:,nt_n      ) )
-      CALL iom_get( numrar, jpdom_auto, 'avm_abl', avm_abl(:,:,:           ), kfill = jpfillcopy )
+      CALL iom_get( numrar, jpdom_auto, 'avm_abl', avm_abl(:,:,:           ) ) !, kfill = jpfillcopy )
       CALL iom_get( numrar, jpdom_auto, 'avt_abl', avt_abl(:,:,:           ) )
-      CALL iom_get( numrar, jpdom_auto,'mxld_abl',mxld_abl(:,:,:           ), kfill = jpfillcopy )
+      CALL iom_get( numrar, jpdom_auto,'mxld_abl',mxld_abl(:,:,:           ) ) !, kfill = jpfillcopy )
       CALL iom_get( numrar, jpdom_auto,    'pblh',    pblh(:,:             ), kfill = jpfillcopy )
 
       IF(.NOT.lrxios) CALL iom_delay_rst( 'READ', 'ABL', numrar )   ! read only abl delayed global communication variables
diff --git a/src/ABL/sbcabl.F90 b/src/ABL/sbcabl.F90
index d0eac5483d1a9caeda3b34add824131d3b7d6b2e..cc34c4a19804ecd48a00a353665bcd36fa897969 100644
--- a/src/ABL/sbcabl.F90
+++ b/src/ABL/sbcabl.F90
@@ -261,7 +261,7 @@ CONTAINS
          rest_eq(:,:) = 1._wp
       END IF
       ! T-mask
-      msk_abl(:,:) = tmask(:,:,1)
+      msk_abl(:,:) = tmask(A2D(0),1)
 
       !!-------------------------------------------------------------------------------------------
 
@@ -320,9 +320,11 @@ CONTAINS
       !!---------------------------------------------------------------------
       INTEGER ,         INTENT(in) ::   kt   ! ocean time step
       !!
-      REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zlat, zevp
+      !REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zlat, zevp
+      REAL(wp), DIMENSION(A2D(0))  ::   zssq, zcd_du, zsen, zlat, zevp
 #if defined key_si3
-      REAL(wp), DIMENSION(jpi,jpj) ::   zssqi, zcd_dui, zseni, zevpi
+      !REAL(wp), DIMENSION(jpi,jpj) ::   zssqi, zcd_dui, zseni, zevpi
+      REAL(wp), DIMENSION(A2D(0))  ::   zssqi, zcd_dui, zseni, zevpi
 #endif
       INTEGER                      ::   jbak, jbak_dta, ji, jj
       !!---------------------------------------------------------------------
@@ -357,7 +359,7 @@ CONTAINS
          !! 3 - Advance ABL variables from now (n) to after (n+1)
          !!-------------------------------------------------------------------------------------------
 
-         CALL abl_stp( kt, tsk_m, ssu_m, ssv_m, zssq,                          &   !   <<= in
+         CALL abl_stp( kt, tsk_m(A2D(0)), ssu_m, ssv_m, zssq,                  &   !   <<= in
             &              sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:),   &   !   <<= in
             &              sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:),   &   !   <<= in
             &              sf(jp_slp )%fnow(:,:,1),                            &   !   <<= in
@@ -366,7 +368,7 @@ CONTAINS
             &              zlat, wndm, utau, vtau, taum                        &   !   =>> out
 #if defined key_si3
             &            , tm_su, u_ice, v_ice, zssqi, zcd_dui                 &   !   <<= in
-            &            , zseni, zevpi, wndm_ice, ato_i                       &   !   <<= in
+            &            , zseni, zevpi, wndm_ice, ato_i(A2D(0))               &   !   <<= in
             &            , utau_ice, vtau_ice                                  &   !   =>> out
 #endif
             &                                                                  )
diff --git a/src/OCE/SBC/sbcblk.F90 b/src/OCE/SBC/sbcblk.F90
index ba5255c909b32dffe6939c279d6c29c49abae526..118951c9ab13ab29657ce8562826dd8b1a6451e2 100644
--- a/src/OCE/SBC/sbcblk.F90
+++ b/src/OCE/SBC/sbcblk.F90
@@ -559,6 +559,7 @@ CONTAINS
          ! Potential temperature of air at z=rn_zqt (most reanalysis products provide absolute temp., not potential temp.)
          IF( ln_tair_pot ) THEN
             ! temperature read into file is already potential temperature, do nothing...
+            IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => air temperature already converted to POTENTIAL!'
             theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1)
          ELSE
             ! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...)
diff --git a/src/OCE/ZDF/zdfgls.F90 b/src/OCE/ZDF/zdfgls.F90
index f44ccc814faf6b6e84f425af23e6f0dfbb449c5c..798fdaf49ef0ee75ff4f5515d28b330d94a2a0ce 100644
--- a/src/OCE/ZDF/zdfgls.F90
+++ b/src/OCE/ZDF/zdfgls.F90
@@ -973,7 +973,8 @@ CONTAINS
                CALL ctl_stop( 'zdf_gls_init: wrong value for nn_mxlice, should be 0,1,2,3 ')
          END SELECT
          IF     ( (nn_mxlice>0).AND.(nn_ice==0) ) THEN
-            CALL ctl_stop( 'zdf_gls_init: with no ice at all, nn_mxlice must be 0 ') 
+            CALL ctl_warn( 'zdf_gls_init: with no ice at all, nn_mxlice is set to 0 ') 
+            nn_mxlice = 0
          ELSEIF ( (nn_mxlice>1).AND.(nn_ice==1) ) THEN
             CALL ctl_stop( 'zdf_gls_init: with no ice model, nn_mxlice must be 0 or 1')
          ENDIF
diff --git a/src/OCE/ZDF/zdftke.F90 b/src/OCE/ZDF/zdftke.F90
index 60371fbcd60ecb0cfba7b84e8af95601aaa84334..0661ea22d0158b4cca4a54b391ae90232e3e25ad 100644
--- a/src/OCE/ZDF/zdftke.F90
+++ b/src/OCE/ZDF/zdftke.F90
@@ -761,7 +761,8 @@ CONTAINS
                CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4')
             END SELECT
             IF     ( (nn_mxlice>0).AND.(nn_ice==0) ) THEN
-               CALL ctl_stop( 'zdf_tke_init: with no ice at all, nn_mxlice must be 0 ') 
+               CALL ctl_warn( 'zdf_tke_init: with no ice at all, nn_mxlice is set to 0 ')
+               nn_mxlice = 0
             ELSEIF ( (nn_mxlice>1).AND.(nn_ice==1) ) THEN
                CALL ctl_stop( 'zdf_tke_init: with no ice model, nn_mxlice must be 0 or 1')
             ENDIF