Index: NEMO/trunk/src/OCE/ZDF/zdfgls.F90
===================================================================
--- NEMO/trunk/src/OCE/ZDF/zdfgls.F90	(revision 14156)
+++ NEMO/trunk/src/OCE/ZDF/zdfgls.F90	(revision 14834)
@@ -137,8 +137,8 @@
       USE zdf_oce , ONLY : en, avtb, avmb   ! ocean vertical physics
       !!
-      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step
-      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
+      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time step
+      INTEGER                             , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   p_sh2          ! shear production term
+      REAL(wp), DIMENSION(:,:,:)          , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop arguments
@@ -151,33 +151,35 @@
       REAL(wp) ::   gh, gm, shr, dif, zsqen, zavt, zavm !   -      -
       REAL(wp) ::   zmsku, zmskv                        !   -      -
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zdep
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zkar
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs       ! Turbulence fluxed induced by internal waves
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro       ! Surface roughness (surface waves)
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra    ! Tapering of wave breaking under sea ice
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb          ! tke at time before
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hmxl_b      ! mixing length at time before
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eps         ! dissipation rate
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi)
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   psi         ! psi at time now
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zd_lw, zd_up, zdiag   ! lower, upper  and diagonal of the matrix
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zstt, zstm  ! stability function on tracer and momentum
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zdep
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zkar
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zflxs                 ! Turbulence fluxed induced by internal waves
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zhsro                 ! Surface roughness (surface waves)
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zice_fra              ! Tapering of wave breaking under sea ice
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   eb                    ! tke at time before
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   hmxl_b                ! mixing length at time before
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   eps                   ! dissipation rate
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zwall_psi             ! Wall function use in the wb case (ln_sigpsi)
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   psi                   ! psi at time now
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zd_lw, zd_up, zdiag   ! lower, upper  and diagonal of the matrix
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zstt, zstm            ! stability function on tracer and momentum
       !!--------------------------------------------------------------------
       !
       ! Preliminary computing
-
-      ustar2_surf(:,:) = 0._wp   ;         psi(:,:,:) = 0._wp
-      ustar2_top (:,:) = 0._wp   ;   zwall_psi(:,:,:) = 0._wp
-      ustar2_bot (:,:) = 0._wp
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         ustar2_surf(ji,jj) = 0._wp   ;   ustar2_top(ji,jj) = 0._wp   ;   ustar2_bot(ji,jj) = 0._wp
+      END_2D
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+         psi(ji,jj,jk) = 0._wp   ;   zwall_psi(ji,jj,jk) = 0._wp
+      END_3D
 
       SELECT CASE ( nn_z0_ice )
       CASE( 0 )   ;   zice_fra(:,:) = 0._wp
-      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp )
-      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:)
-      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp )
+      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(A2D(nn_hls)) * 10._wp )
+      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(A2D(nn_hls))
+      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp )
       END SELECT
 
       ! Compute surface, top and bottom friction at T-points
-      DO_2D( 0, 0, 0, 0 )          !==  surface ocean friction
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          !==  surface ocean friction
          ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1)   ! surface friction
       END_2D
@@ -186,5 +188,5 @@
       !
       IF( .NOT.ln_drg_OFF ) THEN     !== top/bottom friction   (explicit before friction)
-         DO_2D( 0, 0, 0, 0 )         ! bottom friction (explicit before friction)
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )          ! bottom friction (explicit before friction)
             zmsku = 0.5_wp * ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
             zmskv = 0.5_wp * ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
@@ -193,5 +195,5 @@
          END_2D
          IF( ln_isfcav ) THEN
-            DO_2D( 0, 0, 0, 0 )      ! top friction
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )      ! top friction
                zmsku = 0.5_wp * ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
                zmskv = 0.5_wp * ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
@@ -206,28 +208,35 @@
          zhsro(:,:) = rn_hsro
       CASE ( 1 )             ! Standard Charnock formula
-         zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(:,:) , rn_hsro )
+         zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(A2D(nn_hls)) , rn_hsro )
       CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
 !!gm faster coding : the 2 comment lines should be used
 !!gm         zcof = 2._wp * 0.6_wp / 28._wp
 !!gm         zdep(:,:)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) )  )       ! Wave age (eq. 10)
-         zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) )         ! Wave age (eq. 10)
-         zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro)          ! zhsro = rn_frac_hs * Hsw (eq. 11)
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            zcof = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(ji,jj),rsmall))) )          ! Wave age (eq. 10)
+            zhsro(ji,jj) = MAX(rsbc_zs2 * ustar2_surf(ji,jj) * zcof**1.5, rn_hsro)        ! zhsro = rn_frac_hs * Hsw (eq. 11)
+         END_2D
       CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)
-         zhsro(:,:) = MAX(rn_frac_hs * hsw(:,:), rn_hsro)   ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 )
+         zhsro(:,:) = MAX(rn_frac_hs * hsw(A2D(nn_hls)), rn_hsro)   ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 )
       END SELECT
       !
       ! adapt roughness where there is sea ice
-      zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1)  + (1._wp - tmask(:,:,1))*rn_hsro
-      !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )  !==  Compute dissipation rate  ==!
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         zhsro(ji,jj) = ( (1._wp-zice_fra(ji,jj)) * zhsro(ji,jj) + zice_fra(ji,jj) * rn_hsri )*tmask(ji,jj,1)  + &
+            &           (1._wp - tmask(ji,jj,1))*rn_hsro
+      END_2D
+      !
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )  !==  Compute dissipation rate  ==!
          eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
       END_3D
 
       ! Save tke at before time step
-      eb    (:,:,:) = en    (:,:,:)
-      hmxl_b(:,:,:) = hmxl_n(:,:,:)
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+         eb    (ji,jj,jk) = en    (ji,jj,jk)
+         hmxl_b(ji,jj,jk) = hmxl_n(ji,jj,jk)
+      END_3D
 
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             zup   = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)
             zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) )
@@ -250,5 +259,5 @@
       ! Warning : after this step, en : right hand side of the matrix
 
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
          !
          buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction
@@ -303,39 +312,43 @@
       !
       CASE ( 0 )             ! Dirichlet boundary condition (set e at k=1 & 2)
-      ! First level
-      en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3  )
-      zd_lw(:,:,1) = en(:,:,1)
-      zd_up(:,:,1) = 0._wp
-      zdiag(:,:,1) = 1._wp
-      !
-      ! One level below
-      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) &
-         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp) , rn_emin   )
-      zd_lw(:,:,2) = 0._wp
-      zd_up(:,:,2) = 0._wp
-      zdiag(:,:,2) = 1._wp
-      !
-      !
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            ! First level
+            en   (ji,jj,1) = MAX(  rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3  )
+            zd_lw(ji,jj,1) = en(ji,jj,1)
+            zd_up(ji,jj,1) = 0._wp
+            zdiag(ji,jj,1) = 1._wp
+            !
+            ! One level below
+            en   (ji,jj,2) =  MAX( rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1          &
+               &                             * ((zhsro(ji,jj)+gdepw(ji,jj,2,Kmm)) / zhsro(ji,jj) )**(1.5_wp*ra_sf)  )**r2_3 )
+            zd_lw(ji,jj,2) = 0._wp
+            zd_up(ji,jj,2) = 0._wp
+            zdiag(ji,jj,2) = 1._wp
+         END_2D
+         !
+         !
       CASE ( 1 )             ! Neumann boundary condition (set d(e)/dz)
-      !
-      ! Dirichlet conditions at k=1
-      en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 , rn_emin  )
-      zd_lw(:,:,1) = en(:,:,1)
-      zd_up(:,:,1) = 0._wp
-      zdiag(:,:,1) = 1._wp
-      !
-      ! at k=2, set de/dz=Fw
-      !cbr
-      DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo
-         zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
-         zd_lw(ji,jj,2) = 0._wp
-      END_2D
-      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) ))
-      zflxs(:,:)   = rsbc_tke2 * (1._wp-zice_fra(:,:)) * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &
-          &                    * (  ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:)  )**(1.5_wp*ra_sf)
+         !
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            ! Dirichlet conditions at k=1
+            en   (ji,jj,1) = MAX(  rn_emin , rc02r * ustar2_surf(ji,jj) * (1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1)**r2_3  )
+            zd_lw(ji,jj,1) = en(ji,jj,1)
+            zd_up(ji,jj,1) = 0._wp
+            zdiag(ji,jj,1) = 1._wp
+            !
+            ! at k=2, set de/dz=Fw
+            !cbr
+            ! zdiag zd_lw not defined/used on the halo
+            zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
+            zd_lw(ji,jj,2) = 0._wp
+            !
+            zkar (ji,jj)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj)) ))
+            zflxs(ji,jj)   = rsbc_tke2 * (1._wp-zice_fra(ji,jj)) * ustar2_surf(ji,jj)**1.5_wp * zkar(ji,jj) &
+                &                    * (  ( zhsro(ji,jj)+gdept(ji,jj,1,Kmm) ) / zhsro(ji,jj)  )**(1.5_wp*ra_sf)
 !!gm why not   :                        * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf)
-      en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)
-      !
-      !
+            en(ji,jj,2) = en(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm)
+         END_2D
+         !
+         !
       END SELECT
 
@@ -348,5 +361,5 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
 !!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ????
 !!   With thick deep ocean level thickness, this may be quite large, no ???
@@ -365,6 +378,7 @@
          END_2D
          !
+         ! NOTE: ctl_stop with ln_isfcav when using GLS
          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                itop   = mikt(ji,jj)       ! k   top w-point
                itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
@@ -384,5 +398,5 @@
       CASE ( 1 )             ! Neumman boundary condition
          !
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -398,6 +412,7 @@
             en   (ji,jj,ibot) = z_en
          END_2D
+         ! NOTE: ctl_stop with ln_isfcav when using GLS
          IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                itop   = mikt(ji,jj)       ! k   top w-point
                itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
@@ -420,15 +435,17 @@
       ! ----------------------------------------------------------
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
       END_3D
-      DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
       !                                            ! set the minimum value of tke
-      en(:,:,:) = MAX( en(:,:,:), rn_emin )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin )
+      END_3D
 
       !!----------------------------------------!!
@@ -441,20 +458,20 @@
       !
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
          END_3D
          !
       CASE( 1 )               ! k-eps
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             psi(ji,jj,jk)  = eps(ji,jj,jk)
          END_3D
          !
       CASE( 2 )               ! k-w
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
          END_3D
          !
       CASE( 3 )               ! generic
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn
          END_3D
@@ -469,5 +486,5 @@
       ! Warning : after this step, en : right hand side of the matrix
 
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
          !
          ! psi / k
@@ -516,43 +533,45 @@
       CASE ( 0 )             ! Dirichlet boundary conditions
          !
-         ! Surface value
-         zdep    (:,:)   = zhsro(:,:) * rl_sf ! Cosmetic
-         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-         zd_lw(:,:,1) = psi(:,:,1)
-         zd_up(:,:,1) = 0._wp
-         zdiag(:,:,1) = 1._wp
-         !
-         ! One level below
-         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(:,:,2,Kmm)/zhsro(:,:) )))
-         zdep    (:,:)   = (zhsro(:,:) + gdepw(:,:,2,Kmm)) * zkar(:,:)
-         psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-         zd_lw(:,:,2) = 0._wp
-         zd_up(:,:,2) = 0._wp
-         zdiag(:,:,2) = 1._wp
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            ! Surface value
+            zdep    (ji,jj)   = zhsro(ji,jj) * rl_sf ! Cosmetic
+            psi     (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
+            zd_lw(ji,jj,1) = psi(ji,jj,1)
+            zd_up(ji,jj,1) = 0._wp
+            zdiag(ji,jj,1) = 1._wp
+            !
+            ! One level below
+            zkar    (ji,jj)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(ji,jj,2,Kmm)/zhsro(ji,jj) )))
+            zdep    (ji,jj)   = (zhsro(ji,jj) + gdepw(ji,jj,2,Kmm)) * zkar(ji,jj)
+            psi     (ji,jj,2) = rc0**rpp * en(ji,jj,2)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
+            zd_lw(ji,jj,2) = 0._wp
+            zd_up(ji,jj,2) = 0._wp
+            zdiag(ji,jj,2) = 1._wp
+         END_2D
          !
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
          !
-         ! Surface value: Dirichlet
-         zdep    (:,:)   = zhsro(:,:) * rl_sf
-         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-         zd_lw(:,:,1) = psi(:,:,1)
-         zd_up(:,:,1) = 0._wp
-         zdiag(:,:,1) = 1._wp
-         !
-         ! Neumann condition at k=2
-         DO_2D( 0, 0, 0, 0 )   ! zdiag zd_lw not defined/used on the halo
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            ! Surface value: Dirichlet
+            zdep    (ji,jj)   = zhsro(ji,jj) * rl_sf
+            psi     (ji,jj,1) = rc0**rpp * en(ji,jj,1)**rmm * zdep(ji,jj)**rnn * tmask(ji,jj,1)
+            zd_lw(ji,jj,1) = psi(ji,jj,1)
+            zd_up(ji,jj,1) = 0._wp
+            zdiag(ji,jj,1) = 1._wp
+            !
+            ! Neumann condition at k=2, zdiag zd_lw not defined/used on the halo
             zdiag(ji,jj,2) = zdiag(ji,jj,2) +  zd_lw(ji,jj,2) ! Remove zd_lw from zdiag
             zd_lw(ji,jj,2) = 0._wp
+            !
+            ! Set psi vertical flux at the surface:
+            zkar (ji,jj)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(ji,jj,1,Kmm)/zhsro(ji,jj) )) ! Lengh scale slope
+            zdep (ji,jj)   = ((zhsro(ji,jj) + gdept(ji,jj,1,Kmm)) / zhsro(ji,jj))**(rmm*ra_sf)
+            zflxs(ji,jj)   = (rnn + (1._wp-zice_fra(ji,jj))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(ji,jj)) &
+               &           *(1._wp + (1._wp-zice_fra(ji,jj))*rsbc_tke1*zdep(ji,jj))**(2._wp*rmm/3._wp-1_wp)
+            zdep (ji,jj)   = rsbc_psi1 * (zwall_psi(ji,jj,1)*p_avm(ji,jj,1)+zwall_psi(ji,jj,2)*p_avm(ji,jj,2)) * &
+               &           ustar2_surf(ji,jj)**rmm * zkar(ji,jj)**rnn * (zhsro(ji,jj) + gdept(ji,jj,1,Kmm))**(rnn-1.)
+            zflxs(ji,jj)   = zdep(ji,jj) * zflxs(ji,jj)
+            psi  (ji,jj,2) = psi(ji,jj,2) + zflxs(ji,jj) / e3w(ji,jj,2,Kmm)
          END_2D
-         !
-         ! Set psi vertical flux at the surface:
-         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope
-         zdep (:,:)   = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf)
-         zflxs(:,:)   = (rnn + (1._wp-zice_fra(:,:))*rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:)) &
-            &           *(1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
-         zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * &
-            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.)
-         zflxs(:,:)   = zdep(:,:) * zflxs(:,:)
-         psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm)
          !
       END SELECT
@@ -569,5 +588,5 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot
          !                      ! Balance between the production and the dissipation terms
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -588,5 +607,5 @@
       CASE ( 1 )             ! Neumman boundary condition
          !
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
             ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -616,11 +635,11 @@
       ! ----------------
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
       END_3D
-      DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+      DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
@@ -632,15 +651,15 @@
       !
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin)
          END_3D
          !
       CASE( 1 )               ! k-eps
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             eps(ji,jj,jk) = psi(ji,jj,jk)
          END_3D
          !
       CASE( 2 )               ! k-w
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)
          END_3D
@@ -650,5 +669,5 @@
          zex1  =      ( 1.5_wp + rmm/rnn )
          zex2  = -1._wp / rnn
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
             eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
          END_3D
@@ -658,12 +677,15 @@
       ! Limit dissipation rate under stable stratification
       ! --------------------------------------------------
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
          ! limitation
          eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
          hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
-         ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
-         zrn2 = MAX( rn2(ji,jj,jk), rsmall )
-         IF( ln_length_lim )   hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) )
-      END_3D
+      END_3D
+      IF( ln_length_lim ) THEN        ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
+            zrn2 = MAX( rn2(ji,jj,jk), rsmall )
+            hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) )
+         END_3D
+      ENDIF
 
       !
@@ -674,5 +696,5 @@
       !
       CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             ! zcof =  l²/q²
             zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
@@ -691,5 +713,5 @@
          !
       CASE ( 2, 3 )               ! Canuto stability functions
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
             ! zcof =  l²/q²
             zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
@@ -723,10 +745,10 @@
       ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0)
       zstm(:,:,jpk) = 0.
-      DO_2D( 0, 0, 0, 0 )             ! update bottom with good values
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )             ! update bottom with good values
          zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
       END_2D
 
-      zstt(:,:,  1) = wmask(:,:,  1)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
-      zstt(:,:,jpk) = wmask(:,:,jpk)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
+      zstt(:,:,  1) = wmask(A2D(nn_hls),  1)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
+      zstt(:,:,jpk) = wmask(A2D(nn_hls),jpk)  ! default value not needed but avoid a bug when looking for undefined values (-fpe0)
 
 !!gm should be done for ISF (top boundary cond.)
@@ -738,5 +760,5 @@
       !     later overwritten by surface/bottom boundaries conditions, so we don't really care of p_avm(:,:1) and p_avm(:,:jpk)
       !     for zd_lw and zd_up but they have to be defined to avoid a bug when looking for undefined values (-fpe0)
-      DO_3D( 0, 0, 0, 0, 1, jpk )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
          zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
          zavt  = zsqen * zstt(ji,jj,jk)
@@ -745,5 +767,5 @@
          p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
       END_3D
-      p_avt(:,:,1) = 0._wp
+      p_avt(A2D(nn_hls),1) = 0._wp
       !
       IF(sn_cfctl%l_prtctl) THEN
