diff --git a/cfgs/AGRIF_DEMO/EXPREF/1_context_nemo.xml b/cfgs/AGRIF_DEMO/EXPREF/1_context_nemo.xml
index 87b0d3b2c41627901934387958845a9ac578a9ab..ee852b6043d294582fe869594d56dc18367495c5 100644
--- a/cfgs/AGRIF_DEMO/EXPREF/1_context_nemo.xml
+++ b/cfgs/AGRIF_DEMO/EXPREF/1_context_nemo.xml
@@ -27,7 +27,6 @@
 
 <!-- Files definition -->
     <file_definition src="./file_def_nemo-oce.xml"/>     <!--  NEMO ocean dynamics                     -->
-    <file_definition src="./file_def_nemo-ice.xml"/>     <!--  NEMO ocean sea ice                      -->
     <file_definition src="./file_def_nemo-innerttrc.xml"/> <!--  NEMO ocean inert passive tracer           -->
 
 <!-- Axis definition -->
diff --git a/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg b/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg
index fb7cc82244b1e374e61e1f2418b52c68aa5c4b0f..6dacc95db14c3a1a6e5a3252448bd06bc82bc18d 100644
--- a/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg
+++ b/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg
@@ -85,7 +85,7 @@
                      ! Type of air-sea fluxes 
    ln_blk      = .true.    !  Bulk formulation                          (T => fill namsbc_blk )
                      ! Sea-ice :
-   nn_ice      = 2         !  =0 no ice boundary condition
+   nn_ice      = 0         !  =0 no ice boundary condition
       !                    !  =1 use observed ice-cover                 (  => fill namsbc_iif )
       !                    !  =2 or 3 for SI3 and CICE, respectively
                      ! Misc. options of sbc : 
@@ -361,6 +361,10 @@
 !-----------------------------------------------------------------------
 &namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T)
 !-----------------------------------------------------------------------
+      nn_mxlice = 0           ! type of scaling under sea-ice
+      !                       !    = 0 no scaling under sea-ice
+      !                       !    = 1 scaling with constant sea-ice thickness
+      !                       !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model )
 /
 !!======================================================================
 !!                  ***  Diagnostics namelists  ***                   !!
diff --git a/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg b/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg
index ed652bd07d82e67edda89f63269bbf7241f86ea0..cb10389e928920e65cc06244b4d40d8a147b525c 100644
--- a/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg
+++ b/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg
@@ -298,6 +298,10 @@
 !-----------------------------------------------------------------------
 &namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T)
 !-----------------------------------------------------------------------
+      nn_mxlice = 2           ! type of scaling under sea-ice
+      !                       !    = 0 no scaling under sea-ice
+      !                       !    = 1 scaling with constant sea-ice thickness
+      !                       !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model )
 /
 !!======================================================================
 !!                  ***  Diagnostics namelists  ***                   !!
diff --git a/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg b/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg
index 4b74ffed5fa3efe0875aa6974611b577e395a3b5..5a082933dd7c2d717d7b4c906ffd40aed4708699 100644
--- a/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg
+++ b/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg
@@ -299,6 +299,10 @@
 !-----------------------------------------------------------------------
 &namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T)
 !-----------------------------------------------------------------------
+      nn_mxlice = 2           ! type of scaling under sea-ice
+      !                       !    = 0 no scaling under sea-ice
+      !                       !    = 1 scaling with constant sea-ice thickness
+      !                       !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model )
 /
 !!======================================================================
 !!                  ***  Diagnostics namelists  ***                   !!
diff --git a/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg b/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg
index 63973baf27a0c7fdef34c341cbe6e2d22a3b3853..94bb5c091f998f6dd59cc4c92c6429fd76b3bed3 100644
--- a/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg
+++ b/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg
@@ -357,6 +357,10 @@
 !-----------------------------------------------------------------------
 &namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T)
 !-----------------------------------------------------------------------
+      nn_mxlice = 2           ! type of scaling under sea-ice
+      !                       !    = 0 no scaling under sea-ice
+      !                       !    = 1 scaling with constant sea-ice thickness
+      !                       !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model ) 
 /
 !!======================================================================
 !!                  ***  Diagnostics namelists  ***                   !!
diff --git a/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg b/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
index 7f302e9a0770f7aafe9596bf6d99a9449fe77169..c8353d3ad6ea8628ec7d3748dc60578b8b9c2cd3 100644
--- a/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
+++ b/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
@@ -395,6 +395,10 @@
       !                       !        = 2 add a tke source just at the base of the ML
       !                       !        = 3 as = 1 applied on HF part of the stress           (ln_cpl=T)
       ln_mxhsw    = .false.   !  surface mixing length scale = F(wave height)
+      nn_mxlice = 2           ! type of scaling under sea-ice
+      !                       !    = 0 no scaling under sea-ice
+      !                       !    = 1 scaling with constant sea-ice thickness
+      !                       !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model )
 /
 !-----------------------------------------------------------------------
 &namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T)
diff --git a/cfgs/SHARED/namelist_ref b/cfgs/SHARED/namelist_ref
index cde6a6b1673c16b17825f79c5db7426a70b7a6e6..ae9d42ee79e125df420a6142699d2a3f4de78754 100644
--- a/cfgs/SHARED/namelist_ref
+++ b/cfgs/SHARED/namelist_ref
@@ -1199,7 +1199,7 @@
    !                       !                 = 2 first vertical derivative of mixing length bounded by 1
    !                       !                 = 3 as =2 with distinct dissipative an mixing length scale
    ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
-      nn_mxlice    = 2        ! type of scaling under sea-ice
+      nn_mxlice    = 0        ! type of scaling under sea-ice
       !                       !    = 0 no scaling under sea-ice
       !                       !    = 1 scaling with constant sea-ice thickness
       !                       !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model )
@@ -1245,7 +1245,7 @@
    !                       !           = 1 roughness uses rn_hsri and is weigthed by 1-TANH(10*fr_i)
    !                       !           = 2 roughness uses rn_hsri and is weighted by 1-fr_i
    !                       !           = 3 roughness uses rn_hsri and is weighted by 1-MIN(1,4*fr_i)
-   nn_mxlice     =     1   !  mixing under sea ice
+   nn_mxlice     =     0   !  mixing under sea ice
                            !     = 0 No scaling under sea-ice
                            !     = 1 scaling with constant Ice-ocean roughness (rn_hsri)
                            !     = 2 scaling with mean sea-ice thickness
diff --git a/cfgs/WED025/EXPREF/namelist_cfg b/cfgs/WED025/EXPREF/namelist_cfg
index b4293c93155c90857513313f949ead0749694af6..71602bf7bb6bca7377554cf011a04e04a7f02cce 100644
--- a/cfgs/WED025/EXPREF/namelist_cfg
+++ b/cfgs/WED025/EXPREF/namelist_cfg
@@ -552,6 +552,10 @@
 !-----------------------------------------------------------------------
 &namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T)
 !-----------------------------------------------------------------------
+      nn_mxlice = 2           ! type of scaling under sea-ice
+      !                       !    = 0 no scaling under sea-ice
+      !                       !    = 1 scaling with constant sea-ice thickness
+      !                       !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model )
 /
 !-----------------------------------------------------------------------
 &namzdf_gls    !   GLS vertical diffusion                               (ln_zdfgls =T)
diff --git a/doc/namelists/namzdf_tke b/doc/namelists/namzdf_tke
index 85990fe702a48f44706e7551821e6dfb73671634..d5d10df2c1acffe4a5b3d63741eb18604c2e5989 100644
--- a/doc/namelists/namzdf_tke
+++ b/doc/namelists/namzdf_tke
@@ -13,7 +13,7 @@
    !                       !                 = 2 first vertical derivative of mixing length bounded by 1
    !                       !                 = 3 as =2 with distinct dissipative an mixing length scale
    ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
-      nn_mxlice    = 2        ! type of scaling under sea-ice
+      nn_mxlice    = 0        ! type of scaling under sea-ice
       !                       !    = 0 no scaling under sea-ice
       !                       !    = 1 scaling with constant sea-ice thickness
       !                       !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model )
diff --git a/src/ICE/icectl.F90 b/src/ICE/icectl.F90
index e66fe77985a8370c07564c29a8c63fb9b154e4c8..ad20e6734e4f96454372ef4f8ea548145b473d4a 100644
--- a/src/ICE/icectl.F90
+++ b/src/ICE/icectl.F90
@@ -691,49 +691,52 @@ CONTAINS
       CALL prt_ctl_info(' ========== ')
       CALL prt_ctl_info(' - Cell values : ')
       CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
-      CALL prt_ctl(tab2d_1=e1e2t      , clinfo1=' cell area   :')
-      CALL prt_ctl(tab2d_1=at_i       , clinfo1=' at_i        :')
-      CALL prt_ctl(tab2d_1=ato_i      , clinfo1=' ato_i       :')
-      CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' vt_i        :')
-      CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' vt_s        :')
-      CALL prt_ctl(tab2d_1=divu_i     , clinfo1=' divu_i      :')
-      CALL prt_ctl(tab2d_1=delta_i    , clinfo1=' delta_i     :')
-      CALL prt_ctl(tab2d_1=stress1_i  , clinfo1=' stress1_i   :')
-      CALL prt_ctl(tab2d_1=stress2_i  , clinfo1=' stress2_i   :')
-      CALL prt_ctl(tab2d_1=stress12_i , clinfo1=' stress12_i  :')
-      CALL prt_ctl(tab2d_1=strength   , clinfo1=' strength    :')
-      CALL prt_ctl(tab2d_1=delta_i    , clinfo1=' delta_i     :')
-      CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
+      CALL prt_ctl(tab2d_1=e1e2t      , clinfo1=' cell area   :', mask1=tmask)
+      CALL prt_ctl(tab2d_1=at_i       , clinfo1=' at_i        :', mask1=tmask)
+      CALL prt_ctl(tab2d_1=ato_i      , clinfo1=' ato_i       :', mask1=tmask)
+      CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' vt_i        :', mask1=tmask)
+      CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' vt_s        :', mask1=tmask)
+      CALL prt_ctl(tab2d_1=divu_i     , clinfo1=' divu_i      :', mask1=tmask)
+      CALL prt_ctl(tab2d_1=delta_i    , clinfo1=' delta_i     :', mask1=tmask)
+      CALL prt_ctl(tab2d_1=stress1_i  , clinfo1=' stress1_i   :', mask1=tmask)
+      CALL prt_ctl(tab2d_1=stress2_i  , clinfo1=' stress2_i   :', mask1=tmask)
+      CALL prt_ctl(tab2d_1=stress12_i , clinfo1=' stress12_i  :')   ! should be fmask
+      CALL prt_ctl(tab2d_1=strength   , clinfo1=' strength    :', mask1=tmask)
+      CALL prt_ctl(tab2d_1=delta_i    , clinfo1=' delta_i     :', mask1=tmask)
+      CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' u_ice       :', mask1=umask,   &
+         &         tab2d_2=v_ice      , clinfo2=' v_ice       :', mask2=vmask)
 
       DO jl = 1, jpl
          CALL prt_ctl_info(' ')
          CALL prt_ctl_info(' - Category : ', ivar=jl)
          CALL prt_ctl_info('   ~~~~~~~~~~')
-         CALL prt_ctl(tab2d_1=h_i        (:,:,jl)        , clinfo1= ' h_i         : ')
-         CALL prt_ctl(tab2d_1=h_s        (:,:,jl)        , clinfo1= ' h_s         : ')
-         CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' t_su        : ')
-         CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' t_snow      : ')
-         CALL prt_ctl(tab2d_1=s_i        (:,:,jl)        , clinfo1= ' s_i         : ')
-         CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' o_i         : ')
-         CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' a_i         : ')
-         CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' v_i         : ')
-         CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' v_s         : ')
-         CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' e_snow      : ')
-         CALL prt_ctl(tab2d_1=sv_i       (:,:,jl)        , clinfo1= ' sv_i        : ')
-         CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' oa_i        : ')
+         CALL prt_ctl(tab2d_1=h_i (:,:,jl)  , clinfo1= ' h_i         : ', mask1=tmask)
+         CALL prt_ctl(tab2d_1=h_s (:,:,jl)  , clinfo1= ' h_s         : ', mask1=tmask)
+         CALL prt_ctl(tab2d_1=t_su(:,:,jl)  , clinfo1= ' t_su        : ', mask1=tmask)
+         CALL prt_ctl(tab2d_1=t_s (:,:,1,jl), clinfo1= ' t_snow      : ', mask1=tmask)
+         CALL prt_ctl(tab2d_1=s_i (:,:,jl)  , clinfo1= ' s_i         : ', mask1=tmask)
+         CALL prt_ctl(tab2d_1=o_i (:,:,jl)  , clinfo1= ' o_i         : ', mask1=tmask)
+         CALL prt_ctl(tab2d_1=a_i (:,:,jl)  , clinfo1= ' a_i         : ', mask1=tmask)
+         CALL prt_ctl(tab2d_1=v_i (:,:,jl)  , clinfo1= ' v_i         : ', mask1=tmask)
+         CALL prt_ctl(tab2d_1=v_s (:,:,jl)  , clinfo1= ' v_s         : ', mask1=tmask)
+         CALL prt_ctl(tab2d_1=e_s (:,:,1,jl), clinfo1= ' e_snow      : ', mask1=tmask)
+         CALL prt_ctl(tab2d_1=sv_i(:,:,jl)  , clinfo1= ' sv_i        : ', mask1=tmask)
+         CALL prt_ctl(tab2d_1=oa_i(:,:,jl)  , clinfo1= ' oa_i        : ', mask1=tmask)
 
          DO jk = 1, nlay_i
             CALL prt_ctl_info(' - Layer : ', ivar=jk)
-            CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' t_i       : ')
-            CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' e_i       : ')
+            CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' t_i       : ', mask1=tmask)
+            CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' e_i       : ', mask1=tmask)
          END DO
       END DO
 
       CALL prt_ctl_info(' ')
       CALL prt_ctl_info(' - Stresses : ')
       CALL prt_ctl_info('   ~~~~~~~~~~ ')
-      CALL prt_ctl(tab2d_1=utau       , clinfo1= ' utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
-      CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
+      CALL prt_ctl(tab2d_1=utau       , clinfo1= ' utau      : ', mask1 = umask,  &
+         &         tab2d_2=vtau       , clinfo2= ' vtau      : ', mask2 = vmask)
+      CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' utau_ice  : ', mask1 = umask,  &
+         &         tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ', mask2 = vmask)
 
    END SUBROUTINE ice_prt3D
 
diff --git a/src/OCE/ASM/asminc.F90 b/src/OCE/ASM/asminc.F90
index aba92a15dc73f389de5be7fc51c1a5d1b425b783..54371a62c3b73a6d626d97937dbbf52d20a565b8 100644
--- a/src/OCE/ASM/asminc.F90
+++ b/src/OCE/ASM/asminc.F90
@@ -617,10 +617,10 @@ CONTAINS
 !!gm
 
             IF( ln_zps .AND. .NOT. ln_c1d .AND. .NOT. ln_isfcav)           &
-               &  CALL zps_hde    ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
+               &  CALL zps_hde    ( kt, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
                &                              rhd, gru , grv               )  ! of t, s, rd at the last ocean level
             IF( ln_zps .AND. .NOT. ln_c1d .AND.       ln_isfcav)                       &
-               &  CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi,    &  ! Partial steps for top cell (ISF)
+               &  CALL zps_hde_isf( nit000, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi,    &  ! Partial steps for top cell (ISF)
                &                                  rhd, gru , grv , grui, grvi          )  ! of t, s, rd at the last ocean level
 
             DEALLOCATE( t_bkginc )
diff --git a/src/OCE/DIA/diahth.F90 b/src/OCE/DIA/diahth.F90
index 2cee430d390d650c650df5e2745bb32a5cfa2ce3..0b8650acf6c31eacf0ef56eb870ebf898a7aee5f 100644
--- a/src/OCE/DIA/diahth.F90
+++ b/src/OCE/DIA/diahth.F90
@@ -106,12 +106,13 @@ CONTAINS
       IF( ln_timing )   CALL timing_start('dia_hth')
 
       IF( kt == nit000 ) THEN
-         l_hth = .FALSE.
-         IF(   iom_use( 'mlddzt'   ) .OR. iom_use( 'mldr0_3'  ) .OR. iom_use( 'mldr0_1'  )    .OR.  & 
-            &  iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' )    .OR.  &    
-            &  iom_use( '20d'      ) .OR. iom_use( '26d'      ) .OR. iom_use( '28d'      )    .OR.  &    
-            &  iom_use( 'hc300'    ) .OR. iom_use( 'hc700'    ) .OR. iom_use( 'hc2000'   )    .OR.  &    
-            &  iom_use( 'pycndep'  ) .OR. iom_use( 'tinv'     ) .OR. iom_use( 'depti'    )  ) l_hth = .TRUE.
+         !
+         l_hth = iom_use( 'mlddzt'   ) .OR. iom_use( 'mldr0_3'  ) .OR. iom_use( 'mldr0_1'  )    .OR.  & 
+            &    iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' )    .OR.  &    
+            &    iom_use( '20d'      ) .OR. iom_use( '26d'      ) .OR. iom_use( '28d'      )    .OR.  &    
+            &    iom_use( 'hc300'    ) .OR. iom_use( 'hc700'    ) .OR. iom_use( 'hc2000'   )    .OR.  &    
+            &    iom_use( 'pycndep'  ) .OR. iom_use( 'tinv'     ) .OR. iom_use( 'depti'    )
+         !
          !                                      ! allocate dia_hth array
          IF( l_hth ) THEN 
             IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' )
@@ -124,11 +125,12 @@ CONTAINS
 
       IF( l_hth ) THEN
          !
-         IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
-            ! initialization
-            ztinv  (:,:) = 0._wp  
-            zdepinv(:,:) = 0._wp  
-            zmaxdzT(:,:) = 0._wp  
+         ! initialization
+         IF( iom_use( 'tinv'   ) )   ztinv  (:,:) = 0._wp  
+         IF( iom_use( 'depti'  ) )   zdepinv(:,:) = 0._wp  
+         IF( iom_use( 'mlddzt' ) )   zmaxdzT(:,:) = 0._wp  
+         IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' )   &
+            &                    .OR. iom_use( 'mldr10_3' ) .OR. iom_use( 'pycndep'  ) ) THEN
             DO_2D( 1, 1, 1, 1 )
                zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 
                hth     (ji,jj) = zztmp
@@ -137,6 +139,8 @@ CONTAINS
                zrho10_3(ji,jj) = zztmp
                zpycn   (ji,jj) = zztmp
             END_2D
+         ENDIF
+         IF( iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
             IF( nla10 > 1 ) THEN 
                DO_2D( 1, 1, 1, 1 )
                   zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 
@@ -144,25 +148,9 @@ CONTAINS
                   zrho0_1(ji,jj) = zztmp
                END_2D
             ENDIF
-      
-            ! Preliminary computation
-            ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
-            DO_2D( 1, 1, 1, 1 )
-               IF( tmask(ji,jj,nla10) == 1. ) THEN
-                  zu  =  1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80   * ts(ji,jj,nla10,jp_sal,Kmm)  &
-                     &           - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)   &
-                     &           - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm)
-                  zv  =  5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00   * ts(ji,jj,nla10,jp_sal,Kmm)  &
-                     &           - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)
-                  zut =    11.25 -  0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01   * ts(ji,jj,nla10,jp_sal,Kmm)
-                  zvt =    38.00 -  0.750 * ts(ji,jj,nla10,jp_tem,Kmm)
-                  zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
-                  zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
-               ELSE
-                  zdelr(ji,jj) = 0._wp
-               ENDIF
-            END_2D
-
+         ENDIF
+     
+         IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
             ! ------------------------------------------------------------- !
             ! thermocline depth: strongest vertical gradient of temperature !
             ! turbocline depth (mixing layer depth): avt = zavt5            !
@@ -198,6 +186,25 @@ CONTAINS
          !
          IF(  iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' ) .OR.  &    
             &  iom_use( 'pycndep' ) .OR. iom_use( 'tinv'     ) .OR. iom_use( 'depti'    )  ) THEN
+            !
+            ! Preliminary computation
+            ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
+            DO_2D( 1, 1, 1, 1 )
+               IF( tmask(ji,jj,nla10) == 1. ) THEN
+                  zu  =  1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80   * ts(ji,jj,nla10,jp_sal,Kmm)  &
+                     &           - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)   &
+                     &           - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm)
+                  zv  =  5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00   * ts(ji,jj,nla10,jp_sal,Kmm)  &
+                     &           - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)
+                  zut =    11.25 -  0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01   * ts(ji,jj,nla10,jp_sal,Kmm)
+                  zvt =    38.00 -  0.750 * ts(ji,jj,nla10,jp_tem,Kmm)
+                  zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
+                  zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
+               ELSE
+                  zdelr(ji,jj) = 0._wp
+               ENDIF
+            END_2D
+            !
             ! ------------------------------------------------------------- !
             ! MLD: abs( tn - tn(10m) ) = ztem2                              !
             ! Top of thermocline: tn = tn(10m) - ztem2                      !
diff --git a/src/OCE/DOM/domzgr.F90 b/src/OCE/DOM/domzgr.F90
index 76fc34641055338a14ffb3fec079600941146de4..85f89ed2c9b67d30da99d6ad550cf09c8d31cfdb 100644
--- a/src/OCE/DOM/domzgr.F90
+++ b/src/OCE/DOM/domzgr.F90
@@ -324,9 +324,9 @@ CONTAINS
 #if defined key_qco && key_isf
          DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpk )        ! vertical sum at partial cell xxxx other level  
             IF( jk == k_top(ji,jj) ) THEN                               ! first ocean point : partial cell
-               gdept_0(ji,jj,jk) = gdepw_0(ji,jj,jk  ) + 0.5_wp * e3w_0(ji,jj,jk)   ! = risfdep + 1/2 e3w_0(mikt)
+               pdept(ji,jj,jk) = pdepw(ji,jj,jk  ) + 0.5_wp * pe3w(ji,jj,jk)   ! = risfdep + 1/2 e3w_0(mikt)
             ELSE                                                        !  other levels
-               gdept_0(ji,jj,jk) = gdept_0(ji,jj,jk-1) +          e3w_0(ji,jj,jk) 
+               pdept(ji,jj,jk) = pdept(ji,jj,jk-1) +          pe3w(ji,jj,jk) 
             ENDIF
          END_3D
 #endif
diff --git a/src/OCE/DOM/dtatsd.F90 b/src/OCE/DOM/dtatsd.F90
index 7ffccc4cfbc85e95a31aa1518138450e0fc5329a..5863789cc4d8e4c1f09657a00610197ef9c7c599 100644
--- a/src/OCE/DOM/dtatsd.F90
+++ b/src/OCE/DOM/dtatsd.F90
@@ -199,8 +199,7 @@ CONTAINS
          ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk)
       END_3D
       !
-! JC I think it's more convenient to consider the general sco case as the rule
-!      IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
+      IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==!
          !
          IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
             IF( kt == nit000 .AND. lwp )THEN
@@ -236,31 +235,35 @@ CONTAINS
             ptsd(ji,jj,jpk,jp_sal) = 0._wp
          END_2D
          !
-!      ELSE                                !==   z- or zps- coordinate   ==!
-!         !
-!         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
-!            ptsd(ji,jj,jk,jp_tem) = ptsd(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)    ! Mask
-!            ptsd(ji,jj,jk,jp_sal) = ptsd(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
-!         END_3D
-!         !
-!         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
-!            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-!               ik = mbkt(ji,jj)
-!               IF( ik > 1 ) THEN
-!                  zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-!                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
-!                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
-!               ENDIF
-!               ik = mikt(ji,jj)
-!               IF( ik > 1 ) THEN
-!                  zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
-!                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
-!                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
-!               END IF
-!            END_2D
-!         ENDIF
-!         !
-!      ENDIF
+      ELSE                                !==   z- or zps- coordinate   ==!
+         !
+         ! We must keep this definition in a case different from the general case of s-coordinate as we don't
+         ! want to use "underground" values (levels below ocean bottom) to be able to start the model from
+         ! masked temp and sal (read for example in a restart or in output.init)
+         !
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+            ptsd(ji,jj,jk,jp_tem) = ptsd(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)    ! Mask
+            ptsd(ji,jj,jk,jp_sal) = ptsd(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
+         END_3D
+         !
+         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+               ik = mbkt(ji,jj)
+               IF( ik > 1 ) THEN
+                  zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
+                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
+               ENDIF
+               ik = mikt(ji,jj)
+               IF( ik > 1 ) THEN
+                  zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
+                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
+                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
+               END IF
+            END_2D
+         ENDIF
+         !
+      ENDIF
       !
       IF( .NOT.ln_tsd_dmp ) THEN                   !==   deallocate T & S structure   ==!
          !                                              (data used only for initialisation)
diff --git a/src/OCE/DYN/dynhpg.F90 b/src/OCE/DYN/dynhpg.F90
index 4ca2b7a465845357eb90f4392b603527caa5aa50..145c7f9e17407167fc54774c8fe282f407e31f8b 100644
--- a/src/OCE/DYN/dynhpg.F90
+++ b/src/OCE/DYN/dynhpg.F90
@@ -328,7 +328,7 @@ CONTAINS
       ENDIF
 
       ! Partial steps: Compute NOW horizontal gradient of t, s, rd at the last ocean level
-      CALL zps_hde( kt, Kmm, jpts, ts(:,:,:,:,Kmm), zgtsu, zgtsv, rhd, zgru , zgrv )
+      CALL zps_hde( kt, jpts, ts(:,:,:,:,Kmm), zgtsu, zgtsv, rhd, zgru , zgrv )
 
       ! Local constant initialization
       zcoef0 = - grav * 0.5_wp
diff --git a/src/OCE/DYN/dynvor.F90 b/src/OCE/DYN/dynvor.F90
index 0dde705b2e3a7b94dab93160215c6c541169b732..03dda78ceb4f6ce993419062208a9240e3ae7d75 100644
--- a/src/OCE/DYN/dynvor.F90
+++ b/src/OCE/DYN/dynvor.F90
@@ -1098,7 +1098,7 @@ CONTAINS
             !
          CASE DEFAULT                        !* F-point metric term :   pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f)
             ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 0, 0, 0, 0 )
                di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj  ) - e2v(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
                dj_e1u_2e1e2f(ji,jj) = ( e1u(ji  ,jj+1) - e1u(ji,jj) )  * 0.5 * r1_e1e2f(ji,jj)
             END_2D
diff --git a/src/OCE/IOM/iom.F90 b/src/OCE/IOM/iom.F90
index c137d5737cf3db2ba5eb286e3538b171ea18c25e..d66b13d1952f6e3973963acd142394caca868e29 100644
--- a/src/OCE/IOM/iom.F90
+++ b/src/OCE/IOM/iom.F90
@@ -752,7 +752,7 @@ CONTAINS
       ! =============
       clname   = trim(cdname)
       IF ( .NOT. Agrif_Root() .AND. .NOT. lliof ) THEN
-         iln    = INDEX(clname,'/')
+         iln    = INDEX(clname,'/', BACK=.TRUE.)
          cltmpn = clname(1:iln)
          clname = clname(iln+1:LEN_TRIM(clname))
          clname=TRIM(cltmpn)//TRIM(Agrif_CFixed())//'_'//TRIM(clname)
diff --git a/src/OCE/IOM/iom_nf90.F90 b/src/OCE/IOM/iom_nf90.F90
index 5f664536bf2121f546ba07e4adddf03b4eee8f36..c49d9538e007c312f57e69b2fcb05786be34a42f 100644
--- a/src/OCE/IOM/iom_nf90.F90
+++ b/src/OCE/IOM/iom_nf90.F90
@@ -556,6 +556,7 @@ CONTAINS
       LOGICAL               :: lchunk               ! logical switch to activate chunking and compression
       !                                             ! when appropriate (currently chunking is applied to 4d fields only)
       INTEGER               :: idlv                 ! local variable
+      CHARACTER(LEN=256)    :: ccname               ! local variable
       !---------------------------------------------------------------------
       !
       clinfo = '          iom_nf90_rp0123d, file: '//TRIM(iom_file(kiomid)%name)//', var: '//TRIM(cdvar)
@@ -571,7 +572,10 @@ CONTAINS
          ! define the dimension variables if it is not already done
          DO jd = 1, 2
             CALL iom_nf90_check(NF90_INQUIRE_DIMENSION(if90id,jd,iom_file(kiomid)%cn_var(jd),iom_file(kiomid)%dimsz(jd,jd)),clinfo)
-            CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(iom_file(kiomid)%cn_var(jd)), NF90_FLOAT , (/ 1, 2 /),   &
+            ccname = TRIM(iom_file(kiomid)%cn_var(jd))
+            IF ( ccname == 'x') ccname = 'nav_lon'
+            IF ( ccname == 'y') ccname = 'nav_lat'
+            CALL iom_nf90_check(NF90_DEF_VAR( if90id, ccname, NF90_FLOAT , (/ 1, 2 /),   &
                &                              iom_file(kiomid)%nvid(jd) ), clinfo)
          END DO
          iom_file(kiomid)%dimsz(2,1) = iom_file(kiomid)%dimsz(2,2)   ! second dim of first  variable
diff --git a/src/OCE/LBC/lbc_lnk_pt2pt_generic.h90 b/src/OCE/LBC/lbc_lnk_pt2pt_generic.h90
index 8d336f7d2b613768f35af2fe6a416f6dadf745fc..395a3e93c037cf53ac82eb911726c24926168d3b 100644
--- a/src/OCE/LBC/lbc_lnk_pt2pt_generic.h90
+++ b/src/OCE/LBC/lbc_lnk_pt2pt_generic.h90
@@ -120,7 +120,7 @@
       !                   !                    ip0i ip1i        im1i im0i
       !
       iwewe(:) = (/ jpwe,jpea,jpwe,jpea /)   ;   issnn(:) = (/ jpso,jpso,jpno,jpno /)
-      !cd     sides:     west  east south north      ;   corners: so-we, so-ea, no-we, no-ea
+      !     sides:     west  east south north      ;   corners: so-we, so-ea, no-we, no-ea
       isizei(1:4) = (/ ihls, ihls,  ipi,  ipi /)   ;   isizei(5:8) = ihls              ! i- count
       isizej(1:4) = (/ Nj_0, Nj_0, ihls, ihls /)   ;   isizej(5:8) = ihls              ! j- count
       ishtSi(1:4) = (/ ip1i, im1i, ip0i, ip0i /)   ;   ishtSi(5:8) = ishtSi( iwewe )   ! i- shift send data
diff --git a/src/OCE/LBC/mppini.F90 b/src/OCE/LBC/mppini.F90
index 31efed848812f6ebef740445d854d863cee24273..e5cd21c354c4e766e2e79045f80c25ffeda116bb 100644
--- a/src/OCE/LBC/mppini.F90
+++ b/src/OCE/LBC/mppini.F90
@@ -175,6 +175,7 @@ CONTAINS
          ENDIF
             WRITE(numout,*) '      avoid use of mpi_allgather at the north fold  ln_nnogather = ', ln_nnogather
             WRITE(numout,*) '      halo width (applies to both rows and columns)       nn_hls = ', nn_hls
+            WRITE(numout,*) '      choice of communication method                     nn_comm = ', nn_comm
       ENDIF
       !
       IF(lwm)   WRITE( numond, nammpp )
@@ -854,6 +855,7 @@ CONTAINS
          IF(lwp) THEN
             WRITE(numout,*)
             WRITE(numout,*)  '  -----------------------------------------------------------'
+            CALL FLUSH(numout)
          ENDIF
          CALL mppsync
          CALL mppstop( ld_abort = .TRUE. )
@@ -1152,11 +1154,11 @@ CONTAINS
       !!
       !!----------------------------------------------------------------------
       INTEGER ::   inumsave
-      INTEGER ::   jh
+      INTEGER ::   ji,jj,jh
       INTEGER ::   ipi, ipj
       INTEGER ::   iiwe, iiea, iist, iisz 
       INTEGER ::   ijso, ijno, ijst, ijsz 
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zmsk
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zmsk0, zmsk
       LOGICAL , DIMENSION(Ni_0,Nj_0,1)      ::   lloce
       !!----------------------------------------------------------------------
       !
@@ -1175,13 +1177,21 @@ CONTAINS
          ipi = Ni_0 + 2*jh   ! local domain size
          ipj = Nj_0 + 2*jh
          !
-         ALLOCATE( zmsk(ipi,ipj) )
-         zmsk(jh+1:jh+Ni_0,jh+1:jh+Nj_0) = REAL(COUNT(lloce, dim = 3), wp)   ! define inner domain -> need REAL to use lbclnk
-         CALL lbc_lnk('mppini', zmsk, 'T', 1._wp, khls = jh)                 ! fill halos
-         ! Beware, coastal F points can be used in the code -> we may need communications for these points F points even if tmask = 0
-         ! -> the mask we must use here is equal to 1 as soon as one of the 4 neighbours is oce (sum of the mask, not multiplication)
-         zmsk(jh+1:jh+Ni_0,jh+1:jh+Nj_0) = zmsk(jh+1:jh+Ni_0,jh+1  :jh+Nj_0  ) + zmsk(jh+1+1:jh+Ni_0+1,jh+1  :jh+Nj_0  )   &
-            &                            + zmsk(jh+1:jh+Ni_0,jh+1+1:jh+Nj_0+1) + zmsk(jh+1+1:jh+Ni_0+1,jh+1+1:jh+Nj_0+1)
+         ALLOCATE( zmsk0(ipi,ipj), zmsk(ipi,ipj) )
+         zmsk0(jh+1:jh+Ni_0,jh+1:jh+Nj_0) = REAL(COUNT(lloce, dim = 3), wp)   ! define inner domain -> need REAL to use lbclnk
+         CALL lbc_lnk('mppini', zmsk0, 'T', 1._wp, khls = jh)                 ! fill halos
+         ! Beware about the mask we must use here :
+         DO jj = jh+1, jh+Nj_0
+            DO ji = jh+1, jh+Ni_0
+               zmsk(ji,jj) = zmsk0(ji,jj)   &
+                  !  1) dynvor may use scale factors on i+1 (e2v for di_e2v_2e1e2f) and j+1 (e1u for dj_e1u_2e1e2f) even if land
+                  ! -> the mask must be > 1 if south/west neighbours is oce as we may need to send these arrays to these neighbours
+                  &        + zmsk0(ji-1,jj) + zmsk0(ji,jj-1)   &
+                  !  2) coastal F points can be used, so we may need communications for these points F points even IF tmask = 0
+                  ! -> the mask must be > 1 as soon as one of the 3 neighbours is oce: (i,j+1) (i+1,j) (i+1,j+1)
+                  &        + zmsk0(ji+1,jj) + zmsk0(ji,jj+1) + zmsk0(ji+1,jj+1)
+            END DO
+         END DO
          CALL lbc_lnk('mppini', zmsk, 'T', 1._wp, khls = jh)                 ! fill halos again!
          !        
          iiwe = jh   ;   iiea = Ni_0   ! bottom-left corner - 1 of the sent data
@@ -1206,7 +1216,7 @@ IF( nn_comm == 1 ) THEN       ! SM: NOT WORKING FOR NEIGHBOURHOOD COLLECTIVE COM
          !
          iiwe = iiwe-jh   ;   iiea = iiea+jh   ! bottom-left corner - 1 of the received data
          ijso = ijso-jh   ;   ijno = ijno+jh
-         ! do not send if we send only land points
+         ! do not recv if we recv only land points
          IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijst+1:ijst+ijsz) )) == 0 )   mpiRnei(jh,jpwe) = -1
          IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijst+1:ijst+ijsz) )) == 0 )   mpiRnei(jh,jpea) = -1
          IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijso+1:ijso+jh  ) )) == 0 )   mpiRnei(jh,jpso) = -1
@@ -1225,7 +1235,7 @@ ENDIF
             IF( mpiRnei(jh,jpno) > -1 )   mpiRnei(jh, (/jpnw,jpne/) ) = -1   ! NW and NE corners will be received through North nei
         ENDIF
          !
-         DEALLOCATE( zmsk )
+         DEALLOCATE( zmsk0, zmsk )
          !
          CALL mpp_ini_nc(jh)    ! Initialize/Update communicator for neighbourhood collective communications
          !
diff --git a/src/OCE/SBC/geo2ocean.F90 b/src/OCE/SBC/geo2ocean.F90
index 21e80952ca7db917f0c3d6267fa915c1eae40b5f..182e5e03073cc458a418e34b745b0fef8f355368 100644
--- a/src/OCE/SBC/geo2ocean.F90
+++ b/src/OCE/SBC/geo2ocean.F90
@@ -163,71 +163,71 @@ CONTAINS
          !                  
          zlam = plamt(ji,jj)     ! north pole direction & modulous (at t-point)
          zphi = pphit(ji,jj)
-         zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-         zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
+         zxnpt = 0._wp- 2._wp* COS( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)
+         zynpt = 0._wp- 2._wp* SIN( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)
          znnpt = zxnpt*zxnpt + zynpt*zynpt
          !
          zlam = plamu(ji,jj)     ! north pole direction & modulous (at u-point)
          zphi = pphiu(ji,jj)
-         zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-         zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
+         zxnpu = 0._wp- 2._wp* COS( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)
+         zynpu = 0._wp- 2._wp* SIN( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)
          znnpu = zxnpu*zxnpu + zynpu*zynpu
          !
          zlam = plamv(ji,jj)     ! north pole direction & modulous (at v-point)
          zphi = pphiv(ji,jj)
-         zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-         zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
+         zxnpv = 0._wp- 2._wp* COS( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)
+         zynpv = 0._wp- 2._wp* SIN( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)
          znnpv = zxnpv*zxnpv + zynpv*zynpv
          !
          zlam = plamf(ji,jj)     ! north pole direction & modulous (at f-point)
          zphi = pphif(ji,jj)
-         zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-         zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
+         zxnpf = 0._wp- 2._wp* COS( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)
+         zynpf = 0._wp- 2._wp* SIN( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)
          znnpf = zxnpf*zxnpf + zynpf*zynpf
          !
          zlam = plamv(ji,jj  )   ! j-direction: v-point segment direction (around t-point)
          zphi = pphiv(ji,jj  )
          zlan = plamv(ji,jj-1)
          zphh = pphiv(ji,jj-1)
-         zxvvt =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-            &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-         zyvvt =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-            &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
+         zxvvt =  2._wp* COS( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)   &
+            &  -  2._wp* COS( rad*zlan ) * TAN( rpi/4._wp- rad*zphh/2._wp)
+         zyvvt =  2._wp* SIN( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)   &
+            &  -  2._wp* SIN( rad*zlan ) * TAN( rpi/4._wp- rad*zphh/2._wp)
          znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  )
-         znvvt = MAX( znvvt, 1.e-14 )
+         znvvt = MAX( znvvt, 1.e-14_wp )
          !
          zlam = plamf(ji,jj  )   ! j-direction: f-point segment direction (around u-point)
          zphi = pphif(ji,jj  )
          zlan = plamf(ji,jj-1)
          zphh = pphif(ji,jj-1)
-         zxffu =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-            &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-         zyffu =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-            &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
+         zxffu =  2._wp* COS( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)   &
+            &  -  2._wp* COS( rad*zlan ) * TAN( rpi/4._wp- rad*zphh/2._wp)
+         zyffu =  2._wp* SIN( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)   &
+            &  -  2._wp* SIN( rad*zlan ) * TAN( rpi/4._wp- rad*zphh/2._wp)
          znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu )  )
-         znffu = MAX( znffu, 1.e-14 )
+         znffu = MAX( znffu, 1.e-14_wp )
          !
          zlam = plamf(ji  ,jj)   ! i-direction: f-point segment direction (around v-point)
          zphi = pphif(ji  ,jj)
          zlan = plamf(ji-1,jj)
          zphh = pphif(ji-1,jj)
-         zxffv =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-            &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-         zyffv =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-            &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
+         zxffv =  2._wp* COS( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)   &
+            &  -  2._wp* COS( rad*zlan ) * TAN( rpi/4._wp- rad*zphh/2._wp)
+         zyffv =  2._wp* SIN( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)   &
+            &  -  2._wp* SIN( rad*zlan ) * TAN( rpi/4._wp- rad*zphh/2._wp)
          znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv )  )
-         znffv = MAX( znffv, 1.e-14 )
+         znffv = MAX( znffv, 1.e-14_wp )
          !
          zlam = plamu(ji,jj+1)   ! j-direction: u-point segment direction (around f-point)
          zphi = pphiu(ji,jj+1)
          zlan = plamu(ji,jj  )
          zphh = pphiu(ji,jj  )
-         zxuuf =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-            &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-         zyuuf =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-            &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
+         zxuuf =  2._wp* COS( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)   &
+            &  -  2._wp* COS( rad*zlan ) * TAN( rpi/4._wp- rad*zphh/2._wp)
+         zyuuf =  2._wp* SIN( rad*zlam ) * TAN( rpi/4._wp- rad*zphi/2._wp)   &
+            &  -  2._wp* SIN( rad*zlan ) * TAN( rpi/4._wp- rad*zphh/2._wp)
          znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  )
-         znuuf = MAX( znuuf, 1.e-14 )
+         znuuf = MAX( znuuf, 1.e-14_wp )
          !
          !                       ! cosinus and sinus using dot and cross products
          gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt
@@ -249,21 +249,21 @@ CONTAINS
       ! =============== !
 
       DO_2D( 0, 1, 0, 0 )
-         IF( MOD( ABS( plamv(ji,jj) - plamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
-            gsint(ji,jj) = 0.
-            gcost(ji,jj) = 1.
+         IF( MOD( ABS( plamv(ji,jj) - plamv(ji,jj-1) ), 360._wp) < 1.e-8_wp ) THEN
+            gsint(ji,jj) = 0._wp
+            gcost(ji,jj) = 1._wp
          ENDIF
-         IF( MOD( ABS( plamf(ji,jj) - plamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
-            gsinu(ji,jj) = 0.
-            gcosu(ji,jj) = 1.
+         IF( MOD( ABS( plamf(ji,jj) - plamf(ji,jj-1) ), 360._wp) < 1.e-8_wp ) THEN
+            gsinu(ji,jj) = 0._wp
+            gcosu(ji,jj) = 1._wp
          ENDIF
-         IF(      ABS( pphif(ji,jj) - pphif(ji-1,jj) )         < 1.e-8 ) THEN
-            gsinv(ji,jj) = 0.
+         IF(      ABS( pphif(ji,jj) - pphif(ji-1,jj) )           < 1.e-8_wp ) THEN
+            gsinv(ji,jj) = 0._wp
             gcosv(ji,jj) = 1.
          ENDIF
-         IF( MOD( ABS( plamu(ji,jj) - plamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN
-            gsinf(ji,jj) = 0.
-            gcosf(ji,jj) = 1.
+         IF( MOD( ABS( plamu(ji,jj) - plamu(ji,jj+1) ), 360._wp) < 1.e-8_wp ) THEN
+            gsinf(ji,jj) = 0._wp
+            gcosf(ji,jj) = 1._wp
          ENDIF
       END_2D
 
@@ -367,8 +367,8 @@ CONTAINS
       CHARACTER(len=1)            , INTENT( IN    ) ::  cgrid
       REAL(wp), DIMENSION(jpi,jpj), INTENT(   OUT ) ::  pxx , pyy , pzz
       !!
-      REAL(wp), PARAMETER :: rpi = 3.141592653E0
-      REAL(wp), PARAMETER :: rad = rpi / 180.e0
+      REAL(wp), PARAMETER :: rpi = 3.141592653e0_wp
+      REAL(wp), PARAMETER :: rad = rpi / 180.e0_wp
       INTEGER ::   ig     !
       INTEGER ::   ierr   ! local integer
       !!----------------------------------------------------------------------
diff --git a/src/OCE/SBC/sbcblk.F90 b/src/OCE/SBC/sbcblk.F90
index a686596427053376cce721eb1c0af58acc9ebb99..aa3668c3300bb8d3a46998a76f2f56014332d540 100644
--- a/src/OCE/SBC/sbcblk.F90
+++ b/src/OCE/SBC/sbcblk.F90
@@ -369,6 +369,8 @@ CONTAINS
                sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp   ! use standard pressure in Pa
             ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN
                sf(jfpr)%fnow(:,:,1:ipka) = 0._wp        ! no precip or no snow or no surface currents
+            ELSEIF( jfpr == jp_wndi .OR. jfpr == jp_wndj ) THEN
+               sf(jfpr)%fnow(:,:,1:ipka) = 0._wp
             ELSEIF( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) THEN
                IF( .NOT. ln_abl ) THEN
                   DEALLOCATE( sf(jfpr)%fnow )   ! deallocate as not used in this case
@@ -866,11 +868,11 @@ CONTAINS
          CALL iom_put( "vtau_oce", ztau_j(:,:)*tmask(:,:,1) )  ! vtau at T-points!
 
          IF(sn_cfctl%l_prtctl) THEN
-            CALL prt_ctl( tab2d_1=pssq   , clinfo1=' blk_oce_1: pssq   : ')
-            CALL prt_ctl( tab2d_1=wndm   , clinfo1=' blk_oce_1: wndm   : ')
+            CALL prt_ctl( tab2d_1=pssq   , clinfo1=' blk_oce_1: pssq   : ', mask1=tmask )
+            CALL prt_ctl( tab2d_1=wndm   , clinfo1=' blk_oce_1: wndm   : ', mask1=tmask )
             CALL prt_ctl( tab2d_1=utau   , clinfo1=' blk_oce_1: utau   : ', mask1=umask,   &
                &          tab2d_2=vtau   , clinfo2='            vtau   : ', mask2=vmask )
-            CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd     : ')
+            CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd     : ', mask1=tmask )
          ENDIF
          !
       ENDIF ! ln_blk / ln_abl
@@ -943,8 +945,10 @@ CONTAINS
       qns(:,:) = qns(:,:) * tmask(:,:,1)
       !
 #if defined key_si3
-      qns_oce(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:)                             ! non solar without emp (only needed by SI3)
-      qsr_oce(:,:) = qsr(:,:)
+      IF ( nn_ice == 2 ) THEN
+         qns_oce(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:)                  ! non solar without emp (only needed by SI3)
+         qsr_oce(:,:) = qsr(:,:)
+      ENDIF
 #endif
       !
       CALL iom_put( "rho_air"  , rhoa*tmask(:,:,1) )       ! output air density [kg/m^3]
@@ -965,11 +969,11 @@ CONTAINS
       ENDIF
       !
       IF(sn_cfctl%l_prtctl) THEN
-         CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw  : ')
-         CALL prt_ctl(tab2d_1=psen , clinfo1=' blk_oce_2: psen  : ' )
-         CALL prt_ctl(tab2d_1=plat , clinfo1=' blk_oce_2: plat  : ' )
-         CALL prt_ctl(tab2d_1=qns  , clinfo1=' blk_oce_2: qns   : ' )
-         CALL prt_ctl(tab2d_1=emp  , clinfo1=' blk_oce_2: emp   : ')
+         CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw  : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=psen , clinfo1=' blk_oce_2: psen  : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=plat , clinfo1=' blk_oce_2: plat  : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=qns  , clinfo1=' blk_oce_2: qns   : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=emp  , clinfo1=' blk_oce_2: emp   : ', mask1=tmask )
       ENDIF
       !
    END SUBROUTINE blk_oce_2
@@ -1092,8 +1096,8 @@ CONTAINS
          END_2D
          CALL lbc_lnk( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )
          !
-         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   &
-            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' )
+         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : ', mask1=umask   &
+            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ', mask2=vmask )
       ELSE ! ln_abl
 
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
@@ -1105,7 +1109,7 @@ CONTAINS
 
       ENDIF ! ln_blk  / ln_abl
       !
-      IF(sn_cfctl%l_prtctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
+      IF(sn_cfctl%l_prtctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ', mask1=tmask )
       !
    END SUBROUTINE blk_ice_1
 
@@ -1137,6 +1141,7 @@ CONTAINS
       REAL(wp) ::   zst, zst3, zsq, zsipt    ! local variable
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
       REAL(wp) ::   zztmp, zzblk, zztmp1, z1_rLsub   !   -      -
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zmsk   ! temporary mask for prt_ctl
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqlw        ! long wave heat sensitivity over ice
@@ -1301,12 +1306,23 @@ CONTAINS
       ENDIF
       !
       IF(sn_cfctl%l_prtctl) THEN
-         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
-         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
-         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl)
-         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl)
-         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
-         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
+         ALLOCATE(zmsk(jpi,jpj,jpl))
+         DO jl = 1, jpl
+            zmsk(:,:,jl) = tmask(:,:,1)
+         END DO
+         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', mask1=zmsk,   &
+            &         tab3d_2=z_qsb   , clinfo2=' z_qsb    : '         , mask2=zmsk, kdim=jpl)
+         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', mask1=zmsk,   &
+            &         tab3d_2=dqla_ice, clinfo2=' dqla_ice : '         , mask2=zmsk, kdim=jpl)
+         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice: z_dqsb   : ', mask1=zmsk,   &
+            &         tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : '         , mask2=zmsk, kdim=jpl)
+         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', mask1=zmsk,   &
+            &         tab3d_2=qsr_ice , clinfo2=' qsr_ice  : '         , mask2=zmsk, kdim=jpl)
+         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', mask1=zmsk,   &
+            &         tab3d_2=qns_ice , clinfo2=' qns_ice  : '         , mask2=zmsk, kdim=jpl)
+         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', mask1=tmask,   &
+            &         tab2d_2=sprecip , clinfo2=' sprecip  : '         , mask2=tmask         )
+         DEALLOCATE(zmsk)
       ENDIF
 
       !#LB:
diff --git a/src/OCE/SBC/sbcfwb.F90 b/src/OCE/SBC/sbcfwb.F90
index cc5ca7b32f4a2d0ee887b954e8d9d2978f5af4be..58a17f028f4a7cc404859b8c7c4a66f98cd04406 100644
--- a/src/OCE/SBC/sbcfwb.F90
+++ b/src/OCE/SBC/sbcfwb.F90
@@ -104,11 +104,11 @@ CONTAINS
          ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes
          ! and in case of no melt, it can generate HSSW.
          !
-#if ! defined key_si3 && ! defined key_cice
-         snwice_mass_b(:,:) = 0.e0               ! no sea-ice model is being used : no snow+ice mass
-         snwice_mass  (:,:) = 0.e0
-         snwice_fmass (:,:) = 0.e0
-#endif
+         IF( nn_ice == 0 ) THEN
+            snwice_mass_b(:,:) = 0.e0               ! no sea-ice model is being used : no snow+ice mass
+            snwice_mass  (:,:) = 0.e0
+            snwice_fmass (:,:) = 0.e0
+         ENDIF
          !
       ENDIF
 
diff --git a/src/OCE/SBC/sbcmod.F90 b/src/OCE/SBC/sbcmod.F90
index 220693366245e7a91363b61cbec9ba4547d15560..30ee55c7c0dbb66d89024718ff621af8b2c4f73d 100644
--- a/src/OCE/SBC/sbcmod.F90
+++ b/src/OCE/SBC/sbcmod.F90
@@ -39,6 +39,7 @@ MODULE sbcmod
    USE sbcice_if      ! surface boundary condition: ice-if sea-ice model
 #if defined key_si3
    USE icestp         ! surface boundary condition: SI3 sea-ice model
+   USE ice
 #endif
    USE sbcice_cice    ! surface boundary condition: CICE sea-ice model
    USE sbccpl         ! surface boundary condition: coupled formulation
@@ -325,8 +326,14 @@ CONTAINS
       IF( ln_apr_dyn )    CALL sbc_apr_init              ! Atmo Pressure Forcing initialization
       !
 #if defined key_si3
-      IF( lk_agrif .AND. nn_ice == 0 ) THEN            ! allocate ice arrays in case agrif + ice-model + no-ice in child grid
-                          IF( sbc_ice_alloc() /= 0 )   CALL ctl_stop('STOP', 'sbc_ice_alloc : unable to allocate arrays' )
+      IF( nn_ice == 0 ) THEN
+         ! allocate ice arrays in case agrif + ice-model + no-ice in child grid
+         jpl = 1 ; nlay_i = 1 ; nlay_s = 1
+         IF( sbc_ice_alloc() /= 0 )   CALL ctl_stop('STOP', 'sbc_ice_alloc : unable to allocate arrays' )
+#if defined key_agrif
+         CALL Agrif_Declare_Var_ice  !  "      "   "   "      "  Sea ice
+#endif
+
       ELSEIF( nn_ice == 2 ) THEN
                           CALL ice_init( Kbb, Kmm, Kaa )         ! ICE initialization
       ENDIF
diff --git a/src/OCE/SBC/sbcrnf.F90 b/src/OCE/SBC/sbcrnf.F90
index be3f221b8fbaf4c8d779cdcade213ae84b890cf6..758d4f4ba0c3038309ca5c2a251381b8e3ff438a 100644
--- a/src/OCE/SBC/sbcrnf.F90
+++ b/src/OCE/SBC/sbcrnf.F90
@@ -374,9 +374,9 @@ CONTAINS
          IF( .NOT. sn_dep_rnf%ln_clim ) THEN   ;   WRITE(rn_dep_file, '(a,"_y",i4)' ) TRIM( rn_dep_file ), nyear    ! add year
             IF( sn_dep_rnf%clftyp == 'monthly' )   WRITE(rn_dep_file, '(a,"m",i2)'  ) TRIM( rn_dep_file ), nmonth   ! add month
          ENDIF
-         CALL iom_open ( rn_dep_file, inum )                             ! open file
-         CALL iom_get  ( inum, jpdom_global, sn_dep_rnf%clvar, h_rnf )   ! read the river mouth array
-         CALL iom_close( inum )                                          ! close file
+         CALL iom_open ( rn_dep_file, inum )                                                 ! open file
+         CALL iom_get  ( inum, jpdom_global, sn_dep_rnf%clvar, h_rnf, kfill = jpfillcopy )   ! read the river mouth. no 0 on halos!
+         CALL iom_close( inum )                                                              ! close file
          !
          nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
diff --git a/src/OCE/SBC/sbcssm.F90 b/src/OCE/SBC/sbcssm.F90
index 7ff88f827e6dade6f157e99f8a0ba51006673415..1db7fb00a76649dbac077eb4d7e9f56a6c78a676 100644
--- a/src/OCE/SBC/sbcssm.F90
+++ b/src/OCE/SBC/sbcssm.F90
@@ -239,7 +239,7 @@ CONTAINS
             !
             IF( zf_sbc /= REAL( nn_fsbc, wp ) ) THEN      ! nn_fsbc has changed between 2 runs
                IF(lwp) WRITE(numout,*) '   restart with a change in the frequency of mean from ', zf_sbc, ' to ', nn_fsbc
-               zcoef = REAL( nn_fsbc - 1, wp ) / zf_sbc
+               zcoef = REAL( nn_fsbc - 1, wp ) / ( zf_sbc - 1._wp )   ! zf_sbc /= 1 as it was written in the restart
                ssu_m(:,:) = zcoef * ssu_m(:,:)
                ssv_m(:,:) = zcoef * ssv_m(:,:)
                sst_m(:,:) = zcoef * sst_m(:,:)
diff --git a/src/OCE/TRA/eosbn2.F90 b/src/OCE/TRA/eosbn2.F90
index fcefee7a761644f5a28e11502683c358cfe32d44..81fda25e0249bd92da7d324c7e133ee8beb125d1 100644
--- a/src/OCE/TRA/eosbn2.F90
+++ b/src/OCE/TRA/eosbn2.F90
@@ -44,10 +44,8 @@ MODULE eosbn2
    USE stopts         ! Stochastic T/S fluctuations
    !
    USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+   USE lib_mpp        ! for ctl_stop
    USE prtctl         ! Print control
-   USE lbclnk         ! ocean lateral boundary conditions
    USE timing         ! Timing
 
    IMPLICIT NONE
diff --git a/src/OCE/TRA/traldf_lap_blp.F90 b/src/OCE/TRA/traldf_lap_blp.F90
index 3113c0fcb1f6049b0ea397044249d54029e18418..16e5df16cdfb7ca8c542f8fc40d6061a84979c13 100644
--- a/src/OCE/TRA/traldf_lap_blp.F90
+++ b/src/OCE/TRA/traldf_lap_blp.F90
@@ -240,8 +240,8 @@ CONTAINS
       !
       IF (nn_hls==1) CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp )     ! Lateral boundary conditions (unchanged sign)
       !                                               ! Partial top/bottom cell: GRADh( zlap )
-      IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom
-      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, Kmm, kjpt, zlap, zglu, zglv )              ! only bottom
+      IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom
+      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, kjpt, zlap, zglu, zglv )              ! only bottom
       ENDIF
       !
       SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pt_rhs)  ==!
diff --git a/src/OCE/TRA/traldf_triad.F90 b/src/OCE/TRA/traldf_triad.F90
index 84ccde88b8cb2cc5348a5ef451a7ca004c4a270b..70eddec54f8399bc0e4531d1198843206aa64952 100644
--- a/src/OCE/TRA/traldf_triad.F90
+++ b/src/OCE/TRA/traldf_triad.F90
@@ -21,7 +21,6 @@ MODULE traldf_triad
    USE traldf_iso     ! lateral diffusion (Madec operator)         (tra_ldf_iso routine)
    USE diaptr         ! poleward transport diagnostics
    USE diaar5         ! AR5 diagnostics
-   USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
    !
    USE in_out_manager ! I/O manager
    USE iom            ! I/O library
diff --git a/src/OCE/TRA/zpshde.F90 b/src/OCE/TRA/zpshde.F90
index 67d8dca6f7b590d2f990a7d31c4897f865ac734b..2b786d65e30aee97c385ab35e5f84951a51eb130 100644
--- a/src/OCE/TRA/zpshde.F90
+++ b/src/OCE/TRA/zpshde.F90
@@ -40,11 +40,9 @@ MODULE zpshde
    !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE zps_hde( kt, Kmm, kjpt, pta, pgtu, pgtv,  &
-      &                               prd, pgru, pgrv )
+   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv, prd, pgru, pgrv )
       !!
       INTEGER                     , INTENT(in   )           ::  kt          ! ocean time-step index
-      INTEGER                     , INTENT(in   )           ::  Kmm         ! ocean time level index
       INTEGER                     , INTENT(in   )           ::  kjpt        ! number of tracers
       REAL(wp), DIMENSION(:,:,:,:), INTENT(in   )           ::  pta         ! 4D tracers fields
       REAL(wp), DIMENSION(:,:,:)  , INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts
@@ -56,13 +54,13 @@ CONTAINS
       IF( PRESENT(prd)  ) THEN ; itrd = is_tile(prd)  ; ELSE ; itrd = 0 ; ENDIF
       IF( PRESENT(pgru) ) THEN ; itgr = is_tile(pgru) ; ELSE ; itgr = 0 ; ENDIF
 
-      CALL zps_hde_t( kt, Kmm, kjpt, pta, is_tile(pta), pgtu, pgtv, is_tile(pgtu), &
-         &                           prd, itrd,         pgru, pgrv, itgr )
+      CALL zps_hde_t( kt, kjpt, pta, is_tile(pta), pgtu, pgtv, is_tile(pgtu), &
+         &                      prd, itrd,         pgru, pgrv, itgr )
    END SUBROUTINE zps_hde
 
 
-   SUBROUTINE zps_hde_t( kt, Kmm, kjpt, pta, ktta, pgtu, pgtv, ktgt,   &
-      &                                 prd, ktrd, pgru, pgrv, ktgr )
+   SUBROUTINE zps_hde_t( kt, kjpt, pta, ktta, pgtu, pgtv, ktgt,   &
+      &                            prd, ktrd, pgru, pgrv, ktgr )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE zps_hde  ***
       !!
@@ -87,13 +85,13 @@ CONTAINS
       !!                  |   |____                ____|   |
       !!              ___ |   |   |           ___  |   |   |
       !!
-      !!      case 1->   e3w(i+1,:,:,Kmm) >= e3w(i,:,:,Kmm) ( and e3w(:,j+1,:,Kmm) >= e3w(:,j,:,Kmm) ) then
-      !!          t~ = t(i+1,j  ,k) + (e3w(i+1,j,k,Kmm) - e3w(i,j,k,Kmm)) * dk(Ti+1)/e3w(i+1,j,k,Kmm)
-      !!        ( t~ = t(i  ,j+1,k) + (e3w(i,j+1,k,Kmm) - e3w(i,j,k,Kmm)) * dk(Tj+1)/e3w(i,j+1,k,Kmm)  )
+      !!      case 1->   e3w_0(i+1,:,:) >= e3w_0(i,:,:) ( and e3w_0(:,j+1,:) >= e3w_0(:,j,:) ) then
+      !!          t~ = t(i+1,j  ,k) + (e3w_0(i+1,j,k) - e3w_0(i,j,k)) * dk(Ti+1)/e3w_0(i+1,j,k)
+      !!        ( t~ = t(i  ,j+1,k) + (e3w_0(i,j+1,k) - e3w_0(i,j,k)) * dk(Tj+1)/e3w_0(i,j+1,k)  )
       !!          or
-      !!      case 2->   e3w(i+1,:,:,Kmm) <= e3w(i,:,:,Kmm) ( and e3w(:,j+1,:,Kmm) <= e3w(:,j,:,Kmm) ) then
-      !!          t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i+1,j,k,Kmm)) * dk(Ti)/e3w(i,j,k,Kmm)
-      !!        ( t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i,j+1,k,Kmm)) * dk(Tj)/e3w(i,j,k,Kmm) )
+      !!      case 2->   e3w_0(i+1,:,:) <= e3w_0(i,:,:) ( and e3w_0(:,j+1,:) <= e3w_0(:,j,:) ) then
+      !!          t~ = t(i,j,k) + (e3w_0(i,j,k) - e3w_0(i+1,j,k)) * dk(Ti)/e3w_0(i,j,k)
+      !!        ( t~ = t(i,j,k) + (e3w_0(i,j,k) - e3w_0(i,j+1,k)) * dk(Tj)/e3w_0(i,j,k) )
       !!          Idem for di(s) and dj(s)
       !!
       !!      For rho, we call eos which will compute rd~(t~,s~) at the right
@@ -107,7 +105,6 @@ CONTAINS
       !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points
       !!----------------------------------------------------------------------
       INTEGER                                , INTENT(in   )           ::  kt          ! ocean time-step index
-      INTEGER                                , INTENT(in   )           ::  Kmm         ! ocean time level index
       INTEGER                                , INTENT(in   )           ::  kjpt        ! number of tracers
       INTEGER                                , INTENT(in   )           ::  ktta, ktgt, ktrd, ktgr
       REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(in   )           ::  pta         ! 4D tracers fields
@@ -132,19 +129,18 @@ CONTAINS
          DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )              ! Gradient of density at the last level
             iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
             ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
-!!gm BUG ? when applied to before fields, e3w(:,:,k,Kbb) should be used....
-            ze3wu = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
-            ze3wv = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
+            ze3wu = e3w_0(ji+1,jj  ,iku) - e3w_0(ji,jj,iku)
+            ze3wv = e3w_0(ji  ,jj+1,ikv) - e3w_0(ji,jj,ikv)
             !
             ! i- direction
             IF( ze3wu >= 0._wp ) THEN      ! case 1
-               zmaxu =  ze3wu / e3w(ji+1,jj,iku,Kmm)
+               zmaxu =  ze3wu / e3w_0(ji+1,jj,iku)
                ! interpolated values of tracers
                zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
                ! gradient of  tracers
                pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
             ELSE                           ! case 2
-               zmaxu = -ze3wu / e3w(ji,jj,iku,Kmm)
+               zmaxu = -ze3wu / e3w_0(ji,jj,iku)
                ! interpolated values of tracers
                zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
                ! gradient of tracers
@@ -153,13 +149,13 @@ CONTAINS
             !
             ! j- direction
             IF( ze3wv >= 0._wp ) THEN      ! case 1
-               zmaxv =  ze3wv / e3w(ji,jj+1,ikv,Kmm)
+               zmaxv =  ze3wv / e3w_0(ji,jj+1,ikv)
                ! interpolated values of tracers
                ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
                ! gradient of tracers
                pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
             ELSE                           ! case 2
-               zmaxv =  -ze3wv / e3w(ji,jj,ikv,Kmm)
+               zmaxv =  -ze3wv / e3w_0(ji,jj,ikv)
                ! interpolated values of tracers
                ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
                ! gradient of tracers
@@ -176,13 +172,13 @@ CONTAINS
          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
-            ze3wu  = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
-            ze3wv  = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
-            IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)     ! i-direction: case 1
-            ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)     ! -     -      case 2
+            ze3wu  = e3w_0(ji+1,jj  ,iku) - e3w_0(ji,jj,iku)
+            ze3wv  = e3w_0(ji  ,jj+1,ikv) - e3w_0(ji,jj,ikv)
+            IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_0(ji  ,jj,iku)     ! i-direction: case 1
+            ELSE                        ;   zhi(ji,jj) = gdept_0(ji+1,jj,iku)     ! -     -      case 2
             ENDIF
-            IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)     ! j-direction: case 1
-            ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)     ! -     -      case 2
+            IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_0(ji,jj  ,ikv)     ! j-direction: case 1
+            ELSE                        ;   zhj(ji,jj) = gdept_0(ji,jj+1,ikv)     ! -     -      case 2
             ENDIF
          END_2D
          !
@@ -192,8 +188,8 @@ CONTAINS
          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )              ! Gradient of density at the last level
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
-            ze3wu  = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
-            ze3wv  = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
+            ze3wu  = e3w_0(ji+1,jj  ,iku) - e3w_0(ji,jj,iku)
+            ze3wv  = e3w_0(ji  ,jj+1,ikv) - e3w_0(ji,jj,ikv)
             IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
             ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
             ENDIF
@@ -210,11 +206,10 @@ CONTAINS
    END SUBROUTINE zps_hde_t
 
 
-   SUBROUTINE zps_hde_isf( kt, Kmm, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  &
-      &                                   prd, pgru, pgrv, pgrui, pgrvi )
+   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  &
+      &                              prd, pgru, pgrv, pgrui, pgrvi )
       !!
       INTEGER                     , INTENT(in   )           ::  kt           ! ocean time-step index
-      INTEGER                     , INTENT(in   )           ::  Kmm          ! ocean time level index
       INTEGER                     , INTENT(in   )           ::  kjpt         ! number of tracers
       REAL(wp), DIMENSION(:,:,:,:), INTENT(in   )           ::  pta          ! 4D tracers fields
       REAL(wp), DIMENSION(:,:,:)  , INTENT(  out)           ::  pgtu, pgtv   ! hor. grad. of ptra at u- & v-pts
@@ -229,13 +224,13 @@ CONTAINS
       IF( PRESENT(pgru)  ) THEN ; itgr  = is_tile(pgru)  ; ELSE ; itgr  = 0 ; ENDIF
       IF( PRESENT(pgrui) ) THEN ; itgri = is_tile(pgrui) ; ELSE ; itgri = 0 ; ENDIF
 
-      CALL zps_hde_isf_t( kt, Kmm, kjpt, pta, is_tile(pta), pgtu, pgtv, is_tile(pgtu), pgtui, pgtvi, is_tile(pgtui),  &
-      &                                  prd, itrd,         pgru, pgrv, itgr,          pgrui, pgrvi, itgri )
+      CALL zps_hde_isf_t( kt, kjpt, pta, is_tile(pta), pgtu, pgtv, is_tile(pgtu), pgtui, pgtvi, is_tile(pgtui),  &
+      &                             prd, itrd,         pgru, pgrv, itgr,          pgrui, pgrvi, itgri )
    END SUBROUTINE zps_hde_isf
 
 
-   SUBROUTINE zps_hde_isf_t( kt, Kmm, kjpt, pta, ktta, pgtu, pgtv, ktgt, pgtui, pgtvi, ktgti,  &
-      &                                     prd, ktrd, pgru, pgrv, ktgr, pgrui, pgrvi, ktgri )
+   SUBROUTINE zps_hde_isf_t( kt, kjpt, pta, ktta, pgtu, pgtv, ktgt, pgtui, pgtvi, ktgti,  &
+      &                                prd, ktrd, pgru, pgrv, ktgr, pgrui, pgrvi, ktgri )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE zps_hde_isf  ***
       !!
@@ -261,13 +256,13 @@ CONTAINS
       !!                  |   |____                ____|   |
       !!              ___ |   |   |           ___  |   |   |
       !!
-      !!      case 1->   e3w(i+1,j,k,Kmm) >= e3w(i,j,k,Kmm) ( and e3w(i,j+1,k,Kmm) >= e3w(i,j,k,Kmm) ) then
-      !!          t~ = t(i+1,j  ,k) + (e3w(i+1,j  ,k,Kmm) - e3w(i,j,k,Kmm)) * dk(Ti+1)/e3w(i+1,j  ,k,Kmm)
-      !!        ( t~ = t(i  ,j+1,k) + (e3w(i  ,j+1,k,Kmm) - e3w(i,j,k,Kmm)) * dk(Tj+1)/e3w(i  ,j+1,k,Kmm)  )
+      !!      case 1->   e3w_0(i+1,j,k) >= e3w_0(i,j,k) ( and e3w_0(i,j+1,k) >= e3w_0(i,j,k) ) then
+      !!          t~ = t(i+1,j  ,k) + (e3w_0(i+1,j  ,k) - e3w_0(i,j,k)) * dk(Ti+1)/e3w_0(i+1,j  ,k)
+      !!        ( t~ = t(i  ,j+1,k) + (e3w_0(i  ,j+1,k) - e3w_0(i,j,k)) * dk(Tj+1)/e3w_0(i  ,j+1,k)  )
       !!          or
-      !!      case 2->   e3w(i+1,j,k,Kmm) <= e3w(i,j,k,Kmm) ( and e3w(i,j+1,k,Kmm) <= e3w(i,j,k,Kmm) ) then
-      !!          t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i+1,j  ,k,Kmm)) * dk(Ti)/e3w(i,j,k,Kmm)
-      !!        ( t~ = t(i,j,k) + (e3w(i,j,k,Kmm) - e3w(i  ,j+1,k,Kmm)) * dk(Tj)/e3w(i,j,k,Kmm) )
+      !!      case 2->   e3w_0(i+1,j,k) <= e3w_0(i,j,k) ( and e3w_0(i,j+1,k) <= e3w_0(i,j,k) ) then
+      !!          t~ = t(i,j,k) + (e3w_0(i,j,k) - e3w_0(i+1,j  ,k)) * dk(Ti)/e3w_0(i,j,k)
+      !!        ( t~ = t(i,j,k) + (e3w_0(i,j,k) - e3w_0(i  ,j+1,k)) * dk(Tj)/e3w_0(i,j,k) )
       !!          Idem for di(s) and dj(s)
       !!
       !!      For rho, we call eos which will compute rd~(t~,s~) at the right
@@ -283,7 +278,6 @@ CONTAINS
       !!              - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points
       !!----------------------------------------------------------------------
       INTEGER                                , INTENT(in   )           ::  kt           ! ocean time-step index
-      INTEGER                                , INTENT(in   )           ::  Kmm          ! ocean time level index
       INTEGER                                , INTENT(in   )           ::  kjpt         ! number of tracers
       INTEGER                                , INTENT(in   )           ::  ktta, ktgt, ktgti, ktrd, ktgr, ktgri
       REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(in   )           ::  pta          ! 4D tracers fields
@@ -313,18 +307,18 @@ CONTAINS
 
             iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
             ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
-            ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
-            ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
+            ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku)
+            ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv)
             !
             ! i- direction
             IF( ze3wu >= 0._wp ) THEN      ! case 1
-               zmaxu =  ze3wu / e3w(ji+1,jj,iku,Kmm)
+               zmaxu =  ze3wu / e3w_0(ji+1,jj,iku)
                ! interpolated values of tracers
                zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
                ! gradient of  tracers
                pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
             ELSE                           ! case 2
-               zmaxu = -ze3wu / e3w(ji,jj,iku,Kmm)
+               zmaxu = -ze3wu / e3w_0(ji,jj,iku)
                ! interpolated values of tracers
                zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
                ! gradient of tracers
@@ -333,13 +327,13 @@ CONTAINS
             !
             ! j- direction
             IF( ze3wv >= 0._wp ) THEN      ! case 1
-               zmaxv =  ze3wv / e3w(ji,jj+1,ikv,Kmm)
+               zmaxv =  ze3wv / e3w_0(ji,jj+1,ikv)
                ! interpolated values of tracers
                ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
                ! gradient of tracers
                pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
             ELSE                           ! case 2
-               zmaxv =  -ze3wv / e3w(ji,jj,ikv,Kmm)
+               zmaxv =  -ze3wv / e3w_0(ji,jj,ikv)
                ! interpolated values of tracers
                ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
                ! gradient of tracers
@@ -359,14 +353,14 @@ CONTAINS
 
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
-            ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
-            ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
+            ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku)
+            ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv)
             !
-            IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)    ! i-direction: case 1
-            ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)    ! -     -      case 2
+            IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_0(ji  ,jj,iku)    ! i-direction: case 1
+            ELSE                        ;   zhi(ji,jj) = gdept_0(ji+1,jj,iku)    ! -     -      case 2
             ENDIF
-            IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)    ! j-direction: case 1
-            ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)    ! -     -      case 2
+            IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_0(ji,jj  ,ikv)    ! j-direction: case 1
+            ELSE                        ;   zhj(ji,jj) = gdept_0(ji,jj+1,ikv)    ! -     -      case 2
             ENDIF
 
          END_2D
@@ -379,8 +373,8 @@ CONTAINS
          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
-            ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
-            ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
+            ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku)
+            ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv)
 
             IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
             ELSE                        ;   pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
@@ -404,20 +398,20 @@ CONTAINS
             !
             ! (ISF) case partial step top and bottom in adjacent cell in vertical
             ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
-            ! in this case e3w(i,j,k,Kmm) - e3w(i,j+1,k,Kmm) is not the distance between Tj~ and Tj
+            ! in this case e3w_0(i,j,k) - e3w_0(i,j+1,k) is not the distance between Tj~ and Tj
             ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
-            ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
-            ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm)
+            ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku)
+            ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)
 
             ! i- direction
             IF( ze3wu >= 0._wp ) THEN      ! case 1
-               zmaxu = ze3wu / e3w(ji+1,jj,ikup1,Kmm)
+               zmaxu = ze3wu / e3w_0(ji+1,jj,ikup1)
                ! interpolated values of tracers
                zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) )
                ! gradient of tracers
                pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
             ELSE                           ! case 2
-               zmaxu = - ze3wu / e3w(ji,jj,ikup1,Kmm)
+               zmaxu = - ze3wu / e3w_0(ji,jj,ikup1)
                ! interpolated values of tracers
                zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) )
                ! gradient of  tracers
@@ -426,13 +420,13 @@ CONTAINS
             !
             ! j- direction
             IF( ze3wv >= 0._wp ) THEN      ! case 1
-               zmaxv =  ze3wv / e3w(ji,jj+1,ikvp1,Kmm)
+               zmaxv =  ze3wv / e3w_0(ji,jj+1,ikvp1)
                ! interpolated values of tracers
                ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) )
                ! gradient of tracers
                pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
             ELSE                           ! case 2
-               zmaxv =  - ze3wv / e3w(ji,jj,ikvp1,Kmm)
+               zmaxv =  - ze3wv / e3w_0(ji,jj,ikvp1)
                ! interpolated values of tracers
                ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) )
                ! gradient of tracers
@@ -451,15 +445,15 @@ CONTAINS
 
             iku = miku(ji,jj)
             ikv = mikv(ji,jj)
-            ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
-            ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm)
+            ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku)
+            ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)
             !
-            IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)    ! i-direction: case 1
-            ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)    ! -     -      case 2
+            IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_0(ji  ,jj,iku)    ! i-direction: case 1
+            ELSE                        ;   zhi(ji,jj) = gdept_0(ji+1,jj,iku)    ! -     -      case 2
             ENDIF
 
-            IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)    ! j-direction: case 1
-            ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)    ! -     -      case 2
+            IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_0(ji,jj  ,ikv)    ! j-direction: case 1
+            ELSE                        ;   zhj(ji,jj) = gdept_0(ji,jj+1,ikv)    ! -     -      case 2
             ENDIF
 
          END_2D
@@ -470,8 +464,8 @@ CONTAINS
          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             iku = miku(ji,jj)
             ikv = mikv(ji,jj)
-            ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
-            ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm)
+            ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku)
+            ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)
 
             IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1
             ELSE                      ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj  ,iku) - zri(ji,jj    ) ) ! i: 2
diff --git a/src/OCE/ZDF/zdfgls.F90 b/src/OCE/ZDF/zdfgls.F90
index e3f7c96e82ac51ce504c9cab9d6788fb9692368c..b25670799dae097c9853dfdda687a3743aeea10a 100644
--- a/src/OCE/ZDF/zdfgls.F90
+++ b/src/OCE/ZDF/zdfgls.F90
@@ -967,8 +967,13 @@ CONTAINS
             CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean     sea-ice thickness'
             CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max      sea-ice thickness'
             CASE DEFAULT
-               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 ')
+               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 ') 
+         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
          WRITE(numout,*)
       ENDIF
 
diff --git a/src/OCE/ZDF/zdftke.F90 b/src/OCE/ZDF/zdftke.F90
index 4f4ae8ac309bf8eacdb54115f2ce771c7912f55c..68dd9f6ed89537ea357a57cb963ec863cfb94e1d 100644
--- a/src/OCE/ZDF/zdftke.F90
+++ b/src/OCE/ZDF/zdftke.F90
@@ -761,6 +761,11 @@ CONTAINS
             CASE DEFAULT
                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 ') 
+            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
          ENDIF
          WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc
          WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc
diff --git a/src/OCE/step.F90 b/src/OCE/step.F90
index 99e8a952d32f080c71df9f7c6d1149c6331dfa47..ff2759f8c9abfc5bd28d1eb1b36d3e3a00d526f1 100644
--- a/src/OCE/step.F90
+++ b/src/OCE/step.F90
@@ -189,11 +189,11 @@ CONTAINS
       IF( ln_zps .OR. l_ldfslp ) CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
 
       IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
-            &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
+            &            CALL zps_hde    ( kstp, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
             &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level
 
       IF( ln_zps .AND.       ln_isfcav)                                                &
-            &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+            &            CALL zps_hde_isf( kstp, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
             &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
 
       IF( l_ldfslp ) THEN                             ! slope of lateral mixing
diff --git a/src/OCE/stpmlf.F90 b/src/OCE/stpmlf.F90
index 407e38686c29bdcd88ea282ef622f2945d8bb03e..bae696704a881cbd31d2ec554ec1d9d6d2fe869c 100644
--- a/src/OCE/stpmlf.F90
+++ b/src/OCE/stpmlf.F90
@@ -193,11 +193,11 @@ CONTAINS
       IF( ln_zps .OR. l_ldfslp ) CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
 
       IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
-            &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
+            &            CALL zps_hde    ( kstp, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
             &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level
 
       IF( ln_zps .AND.       ln_isfcav)                                                &
-            &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+            &            CALL zps_hde_isf( kstp, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
             &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
 
       IF( l_ldfslp ) THEN                             ! slope of lateral mixing
diff --git a/src/OCE/stprk3.F90 b/src/OCE/stprk3.F90
index 5b1ca4f34801596683cde6fd69891b7276d0d332..337b55f62beb10f4c27aa80a4aa551ced1d860f7 100644
--- a/src/OCE/stprk3.F90
+++ b/src/OCE/stprk3.F90
@@ -154,20 +154,20 @@ CONTAINS
       !  VERTICAL PHYSICS
 !!st                         CALL zdf_phy( kstp, Nbb, Nnn, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
                          CALL zdf_phy( kstp, Nbb, Nbb, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
-!!gm
+!!gm gdep
       !  LATERAL  PHYSICS
       !
-      IF( l_ldfslp ) THEN                             ! slope of lateral mixing
-!!gm gdep
-                         CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
+      IF( ln_zps .OR. l_ldfslp ) CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
 
-         IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
-            &            CALL zps_hde    ( kstp, Nbb, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
+      IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
+            &            CALL zps_hde    ( kstp, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
             &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level
 
-         IF( ln_zps .AND.       ln_isfcav)                                                &
-            &            CALL zps_hde_isf( kstp, Nbb, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+      IF( ln_zps .AND.       ln_isfcav)                                                &
+            &            CALL zps_hde_isf( kstp, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
             &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
+
+      IF( l_ldfslp ) THEN                             ! slope of lateral mixing
          IF( ln_traldf_triad ) THEN
                          CALL ldf_slp_triad( kstp, Nbb, Nbb )             ! before slope for triad operator
          ELSE
diff --git a/src/OFF/dtadyn.F90 b/src/OFF/dtadyn.F90
index c55f569943ca129e6139f79a06c9b45a1f1695de..253d699db5d2b13d0e708af11cfb31d1515b4c69 100644
--- a/src/OFF/dtadyn.F90
+++ b/src/OFF/dtadyn.F90
@@ -675,10 +675,10 @@ CONTAINS
 
       ! Partial steps: before Horizontal DErivative
       IF( ln_zps  .AND. .NOT. ln_isfcav)                            &
-         &            CALL zps_hde    ( kt, Kmm, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
+         &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
          &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level
       IF( ln_zps .AND.        ln_isfcav)                            &
-         &            CALL zps_hde_isf( kt, Kmm, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF)
+         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF)
          &                                        rhd, gru , grv , grui, grvi )  ! of t, s, rd at the first ocean level
 
          rn2b(:,:,:) = rn2(:,:,:)                ! needed for zdfmxl
diff --git a/src/TOP/TRP/trctrp.F90 b/src/TOP/TRP/trctrp.F90
index 03b07bdc86c1562e9707bcbc5800c4717f47e9ae..f164444ed42f6aefe1428a0e3ce8413d35aaea49 100644
--- a/src/TOP/TRP/trctrp.F90
+++ b/src/TOP/TRP/trctrp.F90
@@ -70,8 +70,8 @@ CONTAINS
          !
          !                                                         ! Partial top/bottom cell: GRADh( trb )  
          IF( ln_zps ) THEN
-            IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kt, Kmm, jptra, tr(:,:,:,:,Kbb), pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom
-            ELSE                 ; CALL zps_hde    ( kt, Kmm, jptra, tr(:,:,:,:,Kbb), gtru, gtrv )                                      !  only bottom
+            IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kt, jptra, tr(:,:,:,:,Kbb), pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom
+            ELSE                 ; CALL zps_hde    ( kt, jptra, tr(:,:,:,:,Kbb), gtru, gtrv )                                      !  only bottom
             ENDIF
          ENDIF
          !
diff --git a/src/TOP/trcini.F90 b/src/TOP/trcini.F90
index 486f5208c7b5f69f06a94adc89a751f3c7e20852..394f1cd16b77cf421edafa76afb8f4bc25316cb7 100644
--- a/src/TOP/trcini.F90
+++ b/src/TOP/trcini.F90
@@ -234,7 +234,6 @@ CONTAINS
       !!                     ***  ROUTINE trc_ini_state ***
       !! ** Purpose :          Initialisation of passive tracer concentration 
       !!----------------------------------------------------------------------
-      USE zpshde          ! partial step: hor. derivative   (zps_hde routine)
       USE trcrst          ! passive tracers restart
       USE trcdta          ! initialisation from files
       !