Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceistate.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceistate.F90	(revision 12582)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceistate.F90	(revision 12583)
@@ -18,5 +18,5 @@
    USE oce            ! dynamics and tracers variables
    USE dom_oce        ! ocean domain
-   USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd 
+   USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd
    USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b
    USE eosbn2         ! equation of state
@@ -60,5 +60,5 @@
    INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m)
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
-   !   
+   !
    !! * Substitutions
 #  include "do_loop_substitute.h90"
@@ -77,7 +77,7 @@
       !!
       !! ** Method  :   This routine will put some ice where ocean
-      !!                is at the freezing point, then fill in ice 
-      !!                state variables using prescribed initial 
-      !!                values in the namelist            
+      !!                is at the freezing point, then fill in ice
+      !!                state variables using prescribed initial
+      !!                values in the namelist
       !!
       !! ** Steps   :   1) Set initial surface and basal temperatures
@@ -91,5 +91,5 @@
       !!              where there is no ice
       !!--------------------------------------------------------------------
-      INTEGER, INTENT(in) :: kt            ! time step 
+      INTEGER, INTENT(in) :: kt            ! time step
       INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
       !
@@ -117,5 +117,5 @@
       ! basal temperature (considered at freezing point)   [Kelvin]
       CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
-      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
+      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)
       !
       ! surface temperature and conductivity
@@ -142,5 +142,5 @@
       e_i (:,:,:,:) = 0._wp
       e_s (:,:,:,:) = 0._wp
-      
+
       ! general fields
       a_i (:,:,:) = 0._wp
@@ -215,5 +215,5 @@
             IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
                &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc.
-               &                              * si(jp_ati)%fnow(:,:,1) 
+               &                              * si(jp_ati)%fnow(:,:,1)
             zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1)
             !
@@ -224,5 +224,5 @@
             !
             ! change the switch for the following
-            WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1) 
+            WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1)
             ELSEWHERE                         ;   zswitch(:,:) = 0._wp
             END WHERE
@@ -231,5 +231,5 @@
             !                          !---------------!
             ! no ice if (sst - Tfreez) >= thresold
-            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp 
+            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp
             ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1)
             END WHERE
@@ -244,5 +244,5 @@
                zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
                ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
-               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
+               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
                zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
             ELSEWHERE
@@ -265,5 +265,5 @@
             zhpnd_ini(:,:) = 0._wp
          ENDIF
-         
+
          !-------------!
          ! fill fields !
@@ -292,5 +292,5 @@
          ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
             &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
-         
+
          ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl)
          CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   &
@@ -338,5 +338,5 @@
          DO jl = 1, jpl
             DO_3D_11_11( 1, nlay_i )
-               t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
+               t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)
                ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
                e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
@@ -354,5 +354,5 @@
          END WHERE
          v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
-          
+
          ! specific temperatures for coupled runs
          tn_ice(:,:,:) = t_su(:,:,:)
@@ -374,45 +374,47 @@
          ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rau0
          !
-         IF( .NOT.ln_linssh ) THEN
-            !
-            WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:)
-            ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
-            !
-            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors                
-               e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:)
-               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
-               e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm)
-            END DO
-            !
-            ! Reconstruction of all vertical scale factors at now and before time-steps
-            ! =========================================================================
-            ! Horizontal scale factor interpolations
-            ! --------------------------------------
-            CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
-            CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
-            CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
-            CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
-            CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
-            ! Vertical scale factor interpolations
-            ! ------------------------------------
-            CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )
-            CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
-            CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
-            CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
-            CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
-            ! t- and w- points depth
-            ! ----------------------
-            !!gm not sure of that....
-            gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
-            gdepw(:,:,1,Kmm) = 0.0_wp
-            gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
-            DO jk = 2, jpk
-               gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk  ,Kmm)
-               gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm)
-               gde3w(:,:,jk) = gdept(:,:,jk  ,Kmm) - ssh (:,:,Kmm)
-            END DO
-         ENDIF
+         IF( .NOT.ln_linssh )   CALL dom_vvl_zgr( Kbb, Kmm, Kaa )   ! interpolation scale factor, depth and water column
+! !!st
+!          IF( .NOT.ln_linssh ) THEN
+!             !
+!             WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:)
+!             ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
+!             !
+!             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
+!                e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:)
+!                e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
+!                e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm)
+!             END DO
+!             !
+!             ! Reconstruction of all vertical scale factors at now and before time-steps
+!             ! =========================================================================
+!             ! Horizontal scale factor interpolations
+!             ! --------------------------------------
+!             CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
+!             CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
+!             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
+!             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
+!             CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
+!             ! Vertical scale factor interpolations
+!             ! ------------------------------------
+!             CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )
+!             CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
+!             CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
+!             CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+!             CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
+!             ! t- and w- points depth
+!             ! ----------------------
+!             !!gm not sure of that....
+!             gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
+!             gdepw(:,:,1,Kmm) = 0.0_wp
+!             gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
+!             DO jk = 2, jpk
+!                gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk  ,Kmm)
+!                gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm)
+!                gde3w(:,:,jk) = gdept(:,:,jk  ,Kmm) - ssh (:,:,Kmm)
+!             END DO
+!          ENDIF
       ENDIF
-      
+
       !------------------------------------
       ! 4) store fields at before time-step
@@ -429,5 +431,5 @@
       v_ice_b(:,:)     = v_ice(:,:)
       ! total concentration is needed for Lupkes parameterizations
-      at_i_b (:,:)     = at_i (:,:) 
+      at_i_b (:,:)     = at_i (:,:)
 
 !!clem: output of initial state should be written here but it is impossible because
@@ -441,8 +443,8 @@
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_istate_init  ***
-      !!        
-      !! ** Purpose :   Definition of initial state of the ice 
-      !!
-      !! ** Method  :   Read the namini namelist and check the parameter 
+      !!
+      !! ** Purpose :   Definition of initial state of the ice
+      !!
+      !! ** Method  :   Read the namini namelist and check the parameter
       !!              values called at the first timestep (nit000)
       !!
@@ -485,5 +487,5 @@
          WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst
          IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN
-            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s 
+            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s
             WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s
             WRITE(numout,*) '      initial ice concentr  in the north-south         rn_ati_ini     = ', rn_ati_ini_n,rn_ati_ini_s
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12582)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12583)
@@ -13,11 +13,11 @@
 
    !!----------------------------------------------------------------------
-   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
-   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
-   !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid
-   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
-   !!   dom_vvl_r3c      : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points
-   !!   dom_vvl_rst      : read/write restart file
-   !!   dom_vvl_ctl      : Check the vvl options
+   !!   dom_qe_init     : define initial vertical scale factors, depths and column thickness
+   !!   dom_qe_sf_nxt   : Compute next vertical scale factors
+   !!   dom_qe_sf_update   : Swap vertical scale factors and update the vertical grid
+   !!   dom_qe_interpol : Interpolate vertical scale factors from one grid point to another
+   !!   dom_qe_r3c      : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points
+   !!   dom_qe_rst      : read/write restart file
+   !!   dom_qe_ctl      : Check the vvl options
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -374,4 +374,8 @@
                               + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
             gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
+            gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
+                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
+            gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
+                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
          END DO
          !
@@ -383,4 +387,6 @@
             gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
             gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
+            gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
+            gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
          END DO
          !
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 12583)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 12583)
@@ -0,0 +1,22 @@
+!!----------------------------------------------------------------------
+!!                    ***  domzgr_substitute.h90   ***
+!!----------------------------------------------------------------------
+!! ** purpose :   substitute fsdep. and fse.., the vert. depth and scale
+!!      factors depending on the vertical coord. used, using CPP macro.
+!!----------------------------------------------------------------------
+!! History :  4.2  !  2020-02  (S. Techene, G. Madec)  star coordinate
+!!----------------------------------------------------------------------
+!! NEMO/OCE 4.2 , NEMO Consortium (2020)
+!! $Id$
+!! Software governed by the CeCILL license (see ./LICENSE)
+!!----------------------------------------------------------------------
+!#if defined key_lf
+#   define  e3t(i,j,k,t)   (e3t_0(i,j,k)*(1.+r3t(i,j,t)*tmask(i,j,k)))
+#   define  e3u(i,j,k,t)   (e3u_0(i,j,k)*(1.+r3u(i,j,t)*umask(i,j,k)))
+#   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1.+r3v(i,j,t)*vmask(i,j,k)))
+#   define  e3f(i,j,k,t)   (e3f_0(i,j,k)*(1.+r3f(i,j,t)*fmask(i,j,k)))
+#   define  e3w(i,j,k,t)   (e3w_0(i,j,k)*(1.+r3t(i,j,t)))
+#   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1.+r3u(i,j,t)))
+#   define  e3vw(i,j,k,t)  (e3vw_0(i,j,k)*(1.+r3v(i,j,t)))
+!#endif
+!!----------------------------------------------------------------------
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfLF.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfLF.F90	(revision 12582)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfLF.F90	(revision 12583)
@@ -116,41 +116,41 @@
       ENDIF
 
-      IF ( ln_dynspg_ts ) THEN
-         ! Ensure below that barotropic velocities match time splitting estimate
-         ! Compute actual transport and replace it with ts estimate at "after" time step
-         zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
-         zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
-         DO jk = 2, jpkm1
-            zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
-            zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
-         END DO
-         DO jk = 1, jpkm1
-            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
-            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
-         END DO
-         !
-         IF( .NOT.ln_bt_fw ) THEN
-            ! Remove advective velocity from "now velocities"
-            ! prior to asselin filtering
-            ! In the forward case, this is done below after asselin filtering
-            ! so that asselin contribution is removed at the same time
-            DO jk = 1, jpkm1
-               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)
-               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
-            END DO
-         ENDIF
-      ENDIF
-
-      ! Update after velocity on domain lateral boundaries
-      ! --------------------------------------------------
-# if defined key_agrif
-      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
-# endif
-      !
-      CALL lbc_lnk_multi( 'dynatfLF', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. )     !* local domain boundaries
-      !
-      !                                !* BDY open boundaries
-      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
-      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
+!       IF ( ln_dynspg_ts ) THEN
+!          ! Ensure below that barotropic velocities match time splitting estimate
+!          ! Compute actual transport and replace it with ts estimate at "after" time step
+!          zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
+!          zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
+!          DO jk = 2, jpkm1
+!             zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
+!             zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
+!          END DO
+!          DO jk = 1, jpkm1
+!             puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
+!             pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
+!          END DO
+!          !
+!          IF( .NOT.ln_bt_fw ) THEN
+!             ! Remove advective velocity from "now velocities"
+!             ! prior to asselin filtering
+!             ! In the forward case, this is done below after asselin filtering
+!             ! so that asselin contribution is removed at the same time
+!             DO jk = 1, jpkm1
+!                puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)
+!                pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
+!             END DO
+!          ENDIF
+!       ENDIF
+!
+!       ! Update after velocity on domain lateral boundaries
+!       ! --------------------------------------------------
+! # if defined key_agrif
+!       CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
+! # endif
+!       !
+!       CALL lbc_lnk_multi( 'dynatfLF', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. )     !* local domain boundaries
+!       !
+!       !                                !* BDY open boundaries
+!       IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
+!       IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
 
 !!$   Do we need a call to bdy_vol here??
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90	(revision 12582)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90	(revision 12583)
@@ -36,5 +36,5 @@
 # if defined key_cice4
    USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  &
-                strocnxT,strocnyT,                               & 
+                strocnxT,strocnyT,                               &
                 sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm,     &
                 fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt,          &
@@ -45,5 +45,5 @@
 #else
    USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  &
-                strocnxT,strocnyT,                               & 
+                strocnxT,strocnyT,                               &
                 sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai,     &
                 fresh_ai,fhocn_ai,fswthru_ai,frzmlt,          &
@@ -70,5 +70,5 @@
    INTEGER             ::   jj_off
 
-   INTEGER , PARAMETER ::   jpfld   = 13   ! maximum number of files to read 
+   INTEGER , PARAMETER ::   jpfld   = 13   ! maximum number of files to read
    INTEGER , PARAMETER ::   jp_snow = 1    ! index of snow file
    INTEGER , PARAMETER ::   jp_rain = 2    ! index of rain file
@@ -109,12 +109,12 @@
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_ice_cice  ***
-      !!                   
-      !! ** Purpose :   update the ocean surface boundary condition via the 
-      !!                CICE Sea Ice Model time stepping 
-      !!
-      !! ** Method  : - Get any extra forcing fields for CICE  
+      !!
+      !! ** Purpose :   update the ocean surface boundary condition via the
+      !!                CICE Sea Ice Model time stepping
+      !!
+      !! ** Method  : - Get any extra forcing fields for CICE
       !!              - Prepare forcing fields
       !!              - CICE model time stepping
-      !!              - call the routine that computes mass and 
+      !!              - call the routine that computes mass and
       !!                heat fluxes at the ice/ocean interface
       !!
@@ -171,7 +171,7 @@
       ! there is no restart file.
       ! Values from a CICE restart file would overwrite this
-      IF( .NOT. ln_rstart ) THEN    
-         CALL nemo2cice( ts(:,:,1,jp_tem,Kmm) , sst , 'T' , 1.) 
-      ENDIF  
+      IF( .NOT. ln_rstart ) THEN
+         CALL nemo2cice( ts(:,:,1,jp_tem,Kmm) , sst , 'T' , 1.)
+      ENDIF
 #endif
 
@@ -233,38 +233,39 @@
 !!gm This should be put elsewhere....   (same remark for limsbc)
 !!gm especially here it is assumed zstar coordinate, but it can be ztilde....
-            IF( .NOT.ln_linssh ) THEN
-               !
-               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-               ENDDO
-               e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb)
-               ! Reconstruction of all vertical scale factors at now and before time-steps
-               ! =============================================================================
-               ! Horizontal scale factor interpolations
-               ! --------------------------------------
-               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
-               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
-               CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
-               CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
-               CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
-               ! Vertical scale factor interpolations
-               ! ------------------------------------
-               CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )
-               CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
-               CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
-               CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
-               CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
-               ! t- and w- points depth
-               ! ----------------------
-               gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
-               gdepw(:,:,1,Kmm) = 0.0_wp
-               gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
-               DO jk = 2, jpk
-                  gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk,Kmm)
-                  gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm)
-                  gde3w(:,:,jk)     = gdept(:,:,jk  ,Kmm) - sshn   (:,:)
-               END DO
-            ENDIF
+            IF( .NOT.ln_linssh )   CALL dom_vvl_zgr( Kbb, Kmm, Kaa )   ! interpolation scale factor, depth and water column
+            ! IF( .NOT.ln_linssh ) THEN
+            !    !
+            !    DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
+            !       e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+            !       e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+            !    ENDDO
+            !    e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb)
+            !    ! Reconstruction of all vertical scale factors at now and before time-steps
+            !    ! =============================================================================
+            !    ! Horizontal scale factor interpolations
+            !    ! --------------------------------------
+            !    CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
+            !    CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
+            !    CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
+            !    CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
+            !    CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
+            !    ! Vertical scale factor interpolations
+            !    ! ------------------------------------
+            !    CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )
+            !    CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
+            !    CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
+            !    CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+            !    CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
+            !    ! t- and w- points depth
+            !    ! ----------------------
+            !    gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
+            !    gdepw(:,:,1,Kmm) = 0.0_wp
+            !    gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
+            !    DO jk = 2, jpk
+            !       gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk,Kmm)
+            !       gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm)
+            !       gde3w(:,:,jk)     = gdept(:,:,jk  ,Kmm) - sshn   (:,:)
+            !    END DO
+            ! ENDIF
          ENDIF
       ENDIF
@@ -272,5 +273,5 @@
    END SUBROUTINE cice_sbc_init
 
-   
+
    SUBROUTINE cice_sbc_in( kt, ksbc )
       !!---------------------------------------------------------------------
@@ -281,5 +282,5 @@
       INTEGER, INTENT(in   ) ::   ksbc ! surface forcing type
       !
-      INTEGER  ::   ji, jj, jl                   ! dummy loop indices      
+      INTEGER  ::   ji, jj, jl                   ! dummy loop indices
       REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zpice
       REAL(wp), DIMENSION(jpi,jpj,ncat) :: ztmpn
@@ -293,8 +294,8 @@
       ztmp(:,:)=0.0
 
-! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on 
+! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
 ! the first time-step)
 
-! forced and coupled case 
+! forced and coupled case
 
       IF( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
@@ -356,5 +357,5 @@
 !  Convert to GBM
             IF(ksbc == jp_flx) THEN
-               ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 
+               ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl)
             ELSE
                ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
@@ -380,6 +381,6 @@
          CALL nemo2cice(ztmp,Tair,'T', 1. )    ! Air temperature (K)
          CALL nemo2cice(ztmp,potT,'T', 1. )    ! Potential temp (K)
-! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows  
-         ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )    
+! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows
+         ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )
                                                ! Constant (101000.) atm pressure assumed
          CALL nemo2cice(ztmp,rhoa,'T', 1. )    ! Air density (kg/m^3)
@@ -389,12 +390,12 @@
          CALL nemo2cice(ztmp,zlvl,'T', 1. )    ! Atmos level height (m)
 
-! May want to check all values are physically realistic (as in CICE routine 
+! May want to check all values are physically realistic (as in CICE routine
 ! prepare_forcing)?
 
 ! Divide shortwave into spectral bands (as in prepare_forcing)
          ztmp(:,:)=qsr_ice(:,:,1)*frcvdr       ! visible direct
-         CALL nemo2cice(ztmp,swvdr,'T', 1. )             
+         CALL nemo2cice(ztmp,swvdr,'T', 1. )
          ztmp(:,:)=qsr_ice(:,:,1)*frcvdf       ! visible diffuse
-         CALL nemo2cice(ztmp,swvdf,'T', 1. )              
+         CALL nemo2cice(ztmp,swvdf,'T', 1. )
          ztmp(:,:)=qsr_ice(:,:,1)*frcidr       ! near IR direct
          CALL nemo2cice(ztmp,swidr,'T', 1. )
@@ -406,12 +407,12 @@
 ! Snowfall
 ! Ensure fsnow is positive (as in CICE routine prepare_forcing)
-      IF( iom_use('snowpre') )   CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit  
-      ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)  
-      CALL nemo2cice(ztmp,fsnow,'T', 1. ) 
+      IF( iom_use('snowpre') )   CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit
+      ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)
+      CALL nemo2cice(ztmp,fsnow,'T', 1. )
 
 ! Rainfall
       IF( iom_use('precip') )   CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
       ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
-      CALL nemo2cice(ztmp,frain,'T', 1. ) 
+      CALL nemo2cice(ztmp,frain,'T', 1. )
 
 ! Freezing/melting potential
@@ -482,5 +483,5 @@
       INTEGER, INTENT( in  ) ::   kt   ! ocean time step
       INTEGER, INTENT( in  ) ::   ksbc ! surface forcing type
-      
+
       INTEGER  ::   ji, jj, jl                 ! dummy loop indices
       REAL(wp), DIMENSION(jpi,jpj) :: ztmp1, ztmp2
@@ -490,6 +491,6 @@
          IF(lwp) WRITE(numout,*)'cice_sbc_out'
       ENDIF
-      
-! x comp of ocean-ice stress 
+
+! x comp of ocean-ice stress
       CALL cice2nemo(strocnx,ztmp1,'F', -1. )
       ss_iou(:,:)=0.0
@@ -500,5 +501,5 @@
       CALL lbc_lnk( 'sbcice_cice', ss_iou , 'U', -1. )
 
-! y comp of ocean-ice stress 
+! y comp of ocean-ice stress
       CALL cice2nemo(strocny,ztmp1,'F', -1. )
       ss_iov(:,:)=0.0
@@ -513,32 +514,32 @@
 ! Combine wind stress and ocean-ice stress
 ! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
-! strocnx and strocny already weighted by ice fraction in CICE so not done here 
+! strocnx and strocny already weighted by ice fraction in CICE so not done here
 
       utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
-      vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)     
- 
-! Also need ice/ocean stress on T points so that taum can be updated 
-! This interpolation is already done in CICE so best to use those values 
-      CALL cice2nemo(strocnxT,ztmp1,'T',-1.) 
-      CALL cice2nemo(strocnyT,ztmp2,'T',-1.) 
- 
-! Update taum with modulus of ice-ocean stress 
-! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here 
-taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2) 
-
-! Freshwater fluxes 
+      vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)
+
+! Also need ice/ocean stress on T points so that taum can be updated
+! This interpolation is already done in CICE so best to use those values
+      CALL cice2nemo(strocnxT,ztmp1,'T',-1.)
+      CALL cice2nemo(strocnyT,ztmp2,'T',-1.)
+
+! Update taum with modulus of ice-ocean stress
+! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here
+taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
+
+! Freshwater fluxes
 
       IF(ksbc == jp_flx) THEN
 ! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip)
 ! What we want here is evap*(1-aice)-tprecip*(1-aice) hence manipulation below
-! Not ideal since aice won't be the same as in the atmosphere.  
+! Not ideal since aice won't be the same as in the atmosphere.
 ! Better to use evap and tprecip? (but for now don't read in evap in this case)
          emp(:,:)  = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
       ELSE IF(ksbc == jp_blk) THEN
-         emp(:,:)  = (1.0-fr_i(:,:))*emp(:,:)        
+         emp(:,:)  = (1.0-fr_i(:,:))*emp(:,:)
       ELSE IF(ksbc == jp_purecpl) THEN
-! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above) 
+! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)
 ! This is currently as required with the coupling fields from the UM atmosphere
-         emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:) 
+         emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:)
       ENDIF
 
@@ -560,5 +561,5 @@
       emp(:,:)=emp(:,:)-ztmp1(:,:)
       fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
-      
+
       CALL lbc_lnk_multi( 'sbcice_cice', emp , 'T', 1., sfx , 'T', 1. )
 
@@ -634,5 +635,5 @@
       !!                    ***  ROUTINE cice_sbc_hadgam  ***
       !! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
-      !! 
+      !!
       !!
       !!---------------------------------------------------------------------
@@ -657,5 +658,5 @@
       CALL cice2nemo(vvel,v_ice,'F', -1. )
       !
-      ! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out  
+      ! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out
       !
       ! Snow and ice thicknesses (CO_2 and CO_3)
@@ -689,8 +690,8 @@
       !!---------------------------------------------------------------------
       !! ** Method  :   READ monthly flux file in NetCDF files
-      !!      
-      !!  snowfall    
-      !!  rainfall    
-      !!  sublimation rate    
+      !!
+      !!  snowfall
+      !!  rainfall
+      !!  sublimation rate
       !!  topmelt (category)
       !!  botmelt (category)
@@ -709,5 +710,5 @@
       TYPE(FLD_N) ::   sn_snow, sn_rain, sn_sblm               ! informations about the fields to be read
       TYPE(FLD_N) ::   sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
-      TYPE(FLD_N) ::   sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5 
+      TYPE(FLD_N) ::   sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
       !!
       NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm,   &
@@ -727,6 +728,6 @@
          !            !    file          ! frequency !  variable    ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   ! landmask
          !            !    name          !  (hours)  !   name       !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      ! file
-         sn_snow = FLD_N( 'snowfall_1m'  ,    -1.    ,  'snowfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    ) 
-         sn_rain = FLD_N( 'rainfall_1m'  ,    -1.    ,  'rainfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    ) 
+         sn_snow = FLD_N( 'snowfall_1m'  ,    -1.    ,  'snowfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
+         sn_rain = FLD_N( 'rainfall_1m'  ,    -1.    ,  'rainfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
          sn_sblm = FLD_N( 'sublim_1m'    ,    -1.    ,  'sublim'    ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
          sn_top1 = FLD_N( 'topmeltn1_1m' ,    -1.    ,  'topmeltn1' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
@@ -754,5 +755,5 @@
          slf_i(jp_bot2) = sn_bot2   ;   slf_i(jp_bot3) = sn_bot3   ;   slf_i(jp_bot4) = sn_bot4
          slf_i(jp_bot5) = sn_bot5
-         
+
          ! set sf structure
          ALLOCATE( sf(jpfld), STAT=ierror )
@@ -792,5 +793,5 @@
       ! control print (if less than 100 time-step asked)
       IF( nitend-nit000 <= 100 .AND. lwp ) THEN
-         WRITE(numout,*) 
+         WRITE(numout,*)
          WRITE(numout,*) '        read forcing fluxes for CICE OK'
          CALL FLUSH(numout)
@@ -802,17 +803,17 @@
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nemo2cice  ***
-      !! ** Purpose :   Transfer field in NEMO array to field in CICE array.  
+      !! ** Purpose :   Transfer field in NEMO array to field in CICE array.
 #if defined key_nemocice_decomp
-      !!             
+      !!
       !!                NEMO and CICE PE sub domains are identical, hence
-      !!                there is no need to gather or scatter data from 
+      !!                there is no need to gather or scatter data from
       !!                one PE configuration to another.
 #else
-      !!                Automatically gather/scatter between 
+      !!                Automatically gather/scatter between
       !!                different processors and blocks
       !! ** Method :    A. Ensure all haloes are filled in NEMO field (pn)
       !!                B. Gather pn into global array (png)
       !!                C. Map png into CICE global array (pcg)
-      !!                D. Scatter pcg to CICE blocks (pc) + update haloes  
+      !!                D. Scatter pcg to CICE blocks (pc) + update haloes
 #endif
       !!---------------------------------------------------------------------
@@ -858,5 +859,5 @@
       IF( jpnij > 1) THEN
          CALL mppsync
-         CALL mppgather (pn,0,png) 
+         CALL mppgather (pn,0,png)
          CALL mppsync
       ELSE
@@ -869,5 +870,5 @@
 ! (may be OK but not 100% sure)
 
-      IF(nproc==0) THEN     
+      IF(nproc==0) THEN
 !        pcg(:,:)=0.0
          DO jn=1,jpnij
@@ -890,5 +891,5 @@
          CASE ( 'T' )
             grid_loc=field_loc_center
-         CASE ( 'F' )                              
+         CASE ( 'F' )
             grid_loc=field_loc_NEcorner
       END SELECT
@@ -897,5 +898,5 @@
          CASE ( -1 )
             field_type=field_type_vector
-         CASE ( 1 )                              
+         CASE ( 1 )
             field_type=field_type_scalar
       END SELECT
@@ -916,9 +917,9 @@
       !! ** Purpose :   Transfer field in CICE array to field in NEMO array.
 #if defined key_nemocice_decomp
-      !!             
+      !!
       !!                NEMO and CICE PE sub domains are identical, hence
-      !!                there is no need to gather or scatter data from 
+      !!                there is no need to gather or scatter data from
       !!                one PE configuration to another.
-#else 
+#else
       !!                Automatically deal with scatter/gather between
       !!                different processors and blocks
@@ -926,5 +927,5 @@
       !!                B. Map pcg into NEMO global array (png)
       !!                C. Scatter png into NEMO field (pn) for each processor
-      !!                D. Ensure all haloes are filled in pn 
+      !!                D. Ensure all haloes are filled in pn
 #endif
       !!---------------------------------------------------------------------
@@ -958,5 +959,5 @@
          CASE ( 'T' )
             grid_loc=field_loc_center
-         CASE ( 'F' )                              
+         CASE ( 'F' )
             grid_loc=field_loc_NEcorner
       END SELECT
@@ -965,5 +966,5 @@
          CASE ( -1 )
             field_type=field_type_vector
-         CASE ( 1 )                              
+         CASE ( 1 )
             field_type=field_type_scalar
       END SELECT
@@ -979,5 +980,5 @@
 #else
 
-!      A. Gather CICE blocks (pc) into global array (pcg) 
+!      A. Gather CICE blocks (pc) into global array (pcg)
 
       CALL gather_global(pcg, pc, 0, distrb_info)
@@ -1005,5 +1006,5 @@
       IF( jpnij > 1) THEN
          CALL mppsync
-         CALL mppscatter (png,0,pn) 
+         CALL mppscatter (png,0,pn)
          CALL mppsync
       ELSE
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfLF.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfLF.F90	(revision 12582)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traatfLF.F90	(revision 12583)
@@ -103,13 +103,13 @@
       ENDIF
 
-      ! Update after tracer on domain lateral boundaries
-      !
-#if defined key_agrif
-      CALL Agrif_tra                     ! AGRIF zoom boundaries
-#endif
-      !                                              ! local domain boundaries  (T-point, unchanged sign)
-      CALL lbc_lnk_multi( 'traatfLF', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )
-      !
-      IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries
+!       ! Update after tracer on domain lateral boundaries
+!       !
+! #if defined key_agrif
+!       CALL Agrif_tra                     ! AGRIF zoom boundaries
+! #endif
+!       !                                              ! local domain boundaries  (T-point, unchanged sign)
+!       CALL lbc_lnk_multi( 'traatfLF', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )
+!       !
+!       IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries
 
       ! set time step size (Euler/Leapfrog)
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90	(revision 12582)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90	(revision 12583)
@@ -43,4 +43,6 @@
    USE traatfLF          ! time filtering                   (tra_atf_lf routine)
    USE dynatfLF          ! time filtering                   (dyn_atf_lf routine)
+   USE dynspg_ts         ! surface pressure gradient: split-explicit scheme (define un_adv)
+   USE bdydyn            ! ocean open boundary conditions (define bdy_dyn)
 
    IMPLICIT NONE
@@ -288,4 +290,6 @@
 !!
 !!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine
+                         CALL zdyn_ts       ( Nnn, Naa, e3u, e3v, uu, vv )                   ! barotrope ajustment
+                         CALL finalize_sbc  ( kstp, Nbb, Naa, uu, vv, ts )                   ! boundary condifions
                          CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height
                          CALL dom_qe_r3c    ( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f )            ! "now" ssh/h_0 ratio from filtrered ssh
@@ -352,4 +356,96 @@
       !
    END SUBROUTINE stplf
+
+
+   SUBROUTINE zdyn_ts (Kmm, Kaa, pe3u, pe3v, puu, pvv)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdyn_ts  ***
+      !!
+      !! ** Purpose :   Finalize after horizontal velocity.
+      !!
+      !! ** Method  : * Ensure after velocities transport matches time splitting
+      !!             estimate (ln_dynspg_ts=T)
+      !!
+      !! ** Action :   puu(Kmm),pvv(Kmm),puu(Kaa),pvv(Kaa)   now and after horizontal velocity
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT(in   ) :: Kmm, Kaa    ! before and after time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) :: pe3u, pe3v   ! scale factors
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
+      !
+      INTEGER  ::   jk   ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF ( ln_dynspg_ts ) THEN
+         ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
+         ! Ensure below that barotropic velocities match time splitting estimate
+         ! Compute actual transport and replace it with ts estimate at "after" time step
+         zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
+         zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
+         DO jk = 2, jpkm1
+            zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
+            zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
+         END DO
+         DO jk = 1, jpkm1
+            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
+            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
+         END DO
+         !
+         IF( .NOT.ln_bt_fw ) THEN
+            ! Remove advective velocity from "now velocities"
+            ! prior to asselin filtering
+            ! In the forward case, this is done below after asselin filtering
+            ! so that asselin contribution is removed at the same time
+            DO jk = 1, jpkm1
+               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)
+               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
+            END DO
+         ENDIF
+         !
+         DEALLOCATE( zue, zve )
+         !
+      ENDIF
+      !
+   END SUBROUTINE zdyn_ts
+
+
+   SUBROUTINE finalize_sbc (kt, Kbb, Kaa, puu, pvv, pts)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE finalize_sbc  ***
+      !!
+      !! ** Purpose :   Apply the boundary condition on the after velocity
+      !!
+      !! ** Method  : * Apply lateral boundary conditions on after velocity
+      !!             at the local domain boundaries through lbc_lnk call,
+      !!             at the one-way open boundaries (ln_bdy=T),
+      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
+      !!
+      !! ** Action :   puu(Kaa),pvv(Kaa)   after horizontal velocity and tracers
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
+      INTEGER                             , INTENT(in   ) :: Kbb, Kaa    ! before and after time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities to be time filtered
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers
+      !
+      ! Update after tracer and velocity on domain lateral boundaries
+      !
+#if defined key_agrif
+            CALL Agrif_tra                     ! AGRIF zoom boundaries
+            CALL Agrif_dyn( kt )               !* AGRIF zoom boundaries
+#endif
+      !                                        ! local domain boundaries  (T-point, unchanged sign)
+      CALL lbc_lnk_multi( 'finalize_sbc', puu(:,:,:,       Kaa), 'U', -1., pvv(:,:,:       ,Kaa), 'V', -1.   &
+                       &                , pts(:,:,:,jp_tem,Kaa), 'T',  1., pts(:,:,:,jp_sal,Kaa), 'T',  1. )     !* local domain boundaries
+      !
+      !                                        !* BDY open boundaries
+      IF( ln_bdy )   THEN
+                               CALL bdy_tra( kt, Kbb, pts,      Kaa )
+         IF( ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
+         IF( ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
+      ENDIF
+      !
+   END SUBROUTINE finalize_sbc
+
    !
    !!======================================================================
