Index: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn3d.F90
===================================================================
--- NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn3d.F90	(revision 11224)
+++ NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn3d.F90	(revision 11234)
@@ -186,5 +186,5 @@
       igrd = 3                      ! Copying tangential velocity into bdy points
       IF( llrim0 ) THEN   ;   ibeg = 1                       ;   iend = idx%nblenrim0(igrd)
-      ELSE               ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)
+      ELSE                ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)
       ENDIF
       DO jb = ibeg, iend
Index: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyvol.F90
===================================================================
--- NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyvol.F90	(revision 11224)
+++ NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyvol.F90	(revision 11234)
@@ -99,5 +99,5 @@
             ii = idx%nbi(jb,jgrd)
             ij = idx%nbj(jb,jgrd)
-            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove
+            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice
             zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij)
          END DO
@@ -106,5 +106,5 @@
             ii = idx%nbi(jb,jgrd)
             ij = idx%nbj(jb,jgrd)
-            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove
+            IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice
             zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1)
          END DO
@@ -128,5 +128,5 @@
                ii = idx%nbi(jb,jgrd)
                ij = idx%nbj(jb,jgrd)
-               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice
+               !IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove ?
                pua2d(ii,ij) = pua2d(ii,ij) - idx%flagu(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii+1,ij)
          END DO
@@ -135,5 +135,5 @@
                ii = idx%nbi(jb,jgrd)
                ij = idx%nbj(jb,jgrd)
-               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! sum : else halo couted twice
+               !IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove ?
                pva2d(ii,ij) = pva2d(ii,ij) - idx%flagv(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii,ij+1)
          END DO
@@ -154,5 +154,5 @@
                   ii = idx%nbi(jb,jgrd)
                   ij = idx%nbj(jb,jgrd)
-                  IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
+                  IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE
                   ztranst = ztranst + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij)
             END DO
@@ -161,5 +161,5 @@
                   ii = idx%nbi(jb,jgrd)
                   ij = idx%nbj(jb,jgrd)
-                  IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
+                  IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE
                   ztranst = ztranst + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1)
             END DO
Index: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90	(revision 11224)
+++ NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90	(revision 11234)
@@ -64,4 +64,6 @@
    USE diatmb          ! Top,middle,bottom output
 
+   USE iom   ! to remove
+
    IMPLICIT NONE
    PRIVATE
@@ -104,8 +106,6 @@
       !
       ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
-      !
       IF( ln_dynvor_een .OR. ln_dynvor_eeT )   &
-         &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
-         &               ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) )
+         &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2)   )
          !
       ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj)                    , STAT=ierr(3) )
@@ -153,14 +153,17 @@
       INTEGER  ::   ikbv, iktv            !   -      -
       REAL(wp) ::   r1_2dt_b, z2dt_bf               ! local scalars
-      REAL(wp) ::   zx1, zx2, zu_spg, zhura, z1_hu  !   -      -
-      REAL(wp) ::   zy1, zy2, zv_spg, zhvra, z1_hv  !   -      -
+      REAL(wp) ::   zx1, zx2, zhura        , z1_hu  !   -      -
+      REAL(wp) ::   zy1, zy2, zhvra        , z1_hv  !   -      -
       REAL(wp) ::   za0, za1, za2, za3              !   -      -
-      REAL(wp) ::   zmdi, zztmp            , z1_ht  !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf
-      REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc
-      REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv
+      REAL(wp) ::   zmdi, zztmp, zldg      , z1_ht  !   -      -
+      REAL(wp) ::   zhu_bck, zhv_bck                !   -      -
+      REAL(wp) ::   zun_save, zvn_save              !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e
+      REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zu_spg, zssh_frc
+      REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zv_spg, zhdiv
       REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e
       REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e
       REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
+      REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV         ! fluxes
       !
       REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_wd_dl = .True. 
@@ -182,9 +185,8 @@
       zwdramp = r_rn_wdmin1               ! simplest ramp 
 !     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
-      !                                         ! reciprocal of baroclinic time step 
-      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt
-      ELSE                                        ;   z2dt_bf = 2.0_wp * rdt
-      ENDIF
-      r1_2dt_b = 1.0_wp / z2dt_bf 
+      !                                         ! inverse of baroclinic time step 
+      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   r1_2dt_b = 1._wp / (         rdt )
+      ELSE                                        ;   r1_2dt_b = 1._wp / ( 2._wp * rdt )
+      ENDIF
       !
       ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart 
@@ -210,134 +212,7 @@
             ll_fw_start =.FALSE.
          ENDIF
-         !
-         ! Set averaging weights and cycle length:
+         !                    ! Set averaging weights and cycle length:
          CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
          !
-      ENDIF
-      !
-      IF( ln_isfcav ) THEN    ! top+bottom friction (ocean cavities)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
-               zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
-            END DO
-         END DO
-      ELSE                    ! bottom friction only
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
-               zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
-            END DO
-         END DO
-      ENDIF
-      !
-      ! Set arrays to remove/compute coriolis trend.
-      ! Do it once at kt=nit000 if volume is fixed, else at each long time step.
-      ! Note that these arrays are also used during barotropic loop. These are however frozen
-      ! although they should be updated in the variable volume case. Not a big approximation.
-      ! To remove this approximation, copy lines below inside barotropic loop
-      ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
-      !
-      IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
-         !
-         SELECT CASE( nvor_scheme )
-         CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme)
-            SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
-            CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    &
-                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp  
-                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
-                  END DO
-               END DO
-            CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      &
-                        &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   &
-                        &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      &
-                        &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   )
-                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
-                  END DO
-               END DO
-            END SELECT
-            CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
-            !
-            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
-            DO jj = 2, jpj
-               DO ji = 2, jpi
-                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
-                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
-                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
-                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
-               END DO
-            END DO
-            !
-         CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme)
-            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
-            DO jj = 2, jpj
-               DO ji = 2, jpi
-                  z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
-                  ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
-                  ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht
-                  ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
-                  ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht
-               END DO
-            END DO
-            !
-         CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT !
-            !
-            zwz(:,:) = 0._wp
-            zhf(:,:) = 0._wp
-            
-!!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed 
-!!gm    A priori a better value should be something like :
-!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 
-!!gm                     divided by the sum of the corresponding mask 
-!!gm 
-!!            
-            IF( .NOT.ln_sco ) THEN
-  
-   !!gm  agree the JC comment  : this should be done in a much clear way
-  
-   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
-   !     Set it to zero for the time being 
-   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
-   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
-   !              ENDIF
-   !              zhf(:,:) = gdepw_0(:,:,jk+1)
-               !
-            ELSE
-               !
-               !zhf(:,:) = hbatf(:,:)
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
-                        &              + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
-                        &       / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          &
-                        &              + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  )
-                  END DO
-               END DO
-            ENDIF
-            !
-            DO jj = 1, jpjm1
-               zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
-            END DO
-            !
-            DO jk = 1, jpkm1
-               DO jj = 1, jpjm1
-                  zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
-               END DO
-            END DO
-            CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
-            ! JC: TBC. hf should be greater than 0 
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array
-               END DO
-            END DO
-            zwz(:,:) = ff_f(:,:) * zwz(:,:)
-         END SELECT
       ENDIF
       !
@@ -348,4 +223,5 @@
          CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
       ENDIF
+      !
                           
       ! -----------------------------------------------------------------------------
@@ -354,26 +230,14 @@
       !      
       !
-      !                                   !* e3*d/dt(Ua) (Vertically integrated)
-      !                                   ! --------------------------------------------------
-      zu_frc(:,:) = 0._wp
-      zv_frc(:,:) = 0._wp
-      !
-      DO jk = 1, jpkm1
-         zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-         zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)         
-      END DO
-      !
-      zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:)
-      zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:)
-      !
-      !
-      !                                   !* baroclinic momentum trend (remove the vertical mean trend)
-      DO jk = 1, jpkm1                    ! -----------------------------------------------------------
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk)
-            END DO
-         END DO
+      !                                   !=  zu_frc =  1/H e3*d/dt(Ua)  =!  (Vertical mean of Ua, the 3D trends)
+      !                                   !  ---------------------------  !
+      zu_frc(:,:) = SUM( e3u_n(:,:,:) * ua(:,:,:) * umask(:,:,:) , DIM=3 ) * r1_hu_n(:,:)
+      zv_frc(:,:) = SUM( e3v_n(:,:,:) * va(:,:,:) * vmask(:,:,:) , DIM=3 ) * r1_hv_n(:,:)
+      !
+      !
+      !                                   !=  Ua => baroclinic trend  =!   (remove its vertical mean)
+      DO jk = 1, jpkm1                    !  ------------------------  !
+         ua(:,:,jk) = ( ua(:,:,jk) - zu_frc(:,:) ) * umask(:,:,jk)
+         va(:,:,jk) = ( va(:,:,jk) - zv_frc(:,:) ) * vmask(:,:,jk)
       END DO
       
@@ -381,113 +245,23 @@
 !!gm  Is it correct to do so ?   I think so...
       
-      
-      !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
-      !                                   ! --------------------------------------------------------
-      !
-      zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes 
-      zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)
-      !
-      SELECT CASE( nvor_scheme )
-      CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1   ! vector opt.
-               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj)                    &
-                  &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   &
-                  &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   )
-                  !
-               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj)                    &
-                  &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   & 
-                  &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   ) 
-            END DO  
-         END DO  
-         !         
-      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
-               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
-               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
-               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
-               ! energy conserving formulation for planetary vorticity term
-               zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-               zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
-            END DO
-         END DO
-         !
-      CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
-                 &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
-               zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
-                 &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
-               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
-               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
-            END DO
-         END DO
-         !
-      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
-                &                                         + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
-                &                                         + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
-                &                                         + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
-               zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
-                &                                         + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
-                &                                         + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
-                &                                         + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
-            END DO
-         END DO
-         !
-      END SELECT
-      !
-      !                                   !* Right-Hand-Side of the barotropic momentum equation
-      !                                   ! ----------------------------------------------------
-      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient
-         IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters
+      !                                   !=  remove 2D Coriolis and pressure gradient trends  =!
+      !                                   !  -------------------------------------------------  !
+      !
+      IF( kt == nit000 .OR. .NOT. ln_linssh )   CALL dyn_cor_2D_init   ! Set zwz, the barotropic Coriolis force coefficient
+      !       ! recompute zwz = f/depth  at every time step for (.NOT.ln_linssh) as the water colomn height changes
+      !
+      !                                         !* 2D Coriolis trends
+      zhU(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes 
+      zhV(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)        ! NB: FULL domain : put a value in last row and column
+      !
+      CALL dyn_cor_2d( ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV,  &   ! <<== in
+         &                                     zu_trd, zv_trd   )   ! ==>> out
+      !
+      IF( .NOT.ln_linssh ) THEN                 !* surface pressure gradient   (variable volume only)
+         !
+         IF( ln_wd_il ) THEN                       ! W/D : limiter applied to spgspg
+            CALL wad_spg( sshn, zcpx, zcpy )          ! Calculating W/D gravity filters, zcpx and zcpy
             DO jj = 2, jpjm1
-               DO ji = 2, jpim1 
-                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
-                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
-                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
-                     &                                                         > rn_wdmin1 + rn_wdmin2
-                  ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( &
-                     &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
-                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
-                  IF(ll_tmp1) THEN
-                     zcpx(ji,jj) = 1.0_wp
-                  ELSEIF(ll_tmp2) THEN
-                     ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-                     zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
-                                 &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
-                     zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
-                  ELSE
-                     zcpx(ji,jj) = 0._wp
-                  ENDIF
-                  !
-                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                &
-                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
-                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
-                     &                                                       > rn_wdmin1 + rn_wdmin2
-                  ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( &
-                     &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
-                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
-  
-                  IF(ll_tmp1) THEN
-                     zcpy(ji,jj) = 1.0_wp
-                  ELSE IF(ll_tmp2) THEN
-                     ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-                     zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
-                        &             / (sshn(ji,jj+1) - sshn(ji,jj  )) )
-                     zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
-                  ELSE
-                     zcpy(ji,jj) = 0._wp
-                  ENDIF
-               END DO
-            END DO
-            !
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
+               DO ji = 2, jpim1                ! SPG with the application of W/D gravity filters
                   zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
                      &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth
@@ -496,7 +270,5 @@
                END DO
             END DO
-            !
-         ELSE
-            !
+         ELSE                                      ! now suface pressure gradient
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -516,75 +288,33 @@
       END DO 
       !
-      !                                         ! Add bottom stress contribution from baroclinic velocities:      
-      IF (ln_bt_fw) THEN
-         DO jj = 2, jpjm1                          
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ikbu = mbku(ji,jj)       
-               ikbv = mbkv(ji,jj)    
-               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
-               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
-            END DO
-         END DO
-      ELSE
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ikbu = mbku(ji,jj)       
-               ikbv = mbkv(ji,jj)    
-               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
-               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
-            END DO
-         END DO
-      ENDIF
-      !
-      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag
-      IF( ln_wd_il ) THEN
-         zztmp = -1._wp / rdtbt
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu_frc(ji,jj) = zu_frc(ji,jj) + & 
-               & MAX(r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ), zztmp ) * zwx(ji,jj) *  wdrampu(ji,jj)
-               zv_frc(ji,jj) = zv_frc(ji,jj) + & 
-               & MAX(r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ), zztmp ) * zwy(ji,jj) *  wdrampv(ji,jj)
-            END DO
-         END DO
-      ELSE
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj)
-               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj)
-            END DO
-         END DO
-      END IF
-      !
-      IF( ln_isfcav ) THEN       ! Add TOP stress contribution from baroclinic velocities:      
-         IF( ln_bt_fw ) THEN
-            DO jj = 2, jpjm1
+      !                                   !=  Add bottom stress contribution from baroclinic velocities  =!
+      !                                   !  -----------------------------------------------------------  !
+      CALL drg_init( zu_frc, zv_frc,  zCdU_u, zCdU_v )      ! also provide the barotropic drag coefficients
+      !
+      !                                   !=  Add atmospheric pressure forcing  =!
+      !                                   !  ----------------------------------  !
+      IF( ln_apr_dyn ) THEN
+         IF( ln_bt_fw ) THEN                          ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
+            DO jj = 2, jpjm1              
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  iktu = miku(ji,jj)
-                  iktv = mikv(ji,jj)
-                  zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
-                  zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
-               END DO
-            END DO
-         ELSE
-            DO jj = 2, jpjm1
+                  zu_frc(ji,jj) = zu_frc(ji,jj) + grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
+                  zv_frc(ji,jj) = zv_frc(ji,jj) + grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
+               END DO
+            END DO
+         ELSE                                         ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
+            zztmp = grav * r1_2
+            DO jj = 2, jpjm1              
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  iktu = miku(ji,jj)
-                  iktv = mikv(ji,jj)
-                  zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
-                  zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
-               END DO
-            END DO
-         ENDIF
-         !
-         ! Note that the "unclipped" top friction parameter is used even with explicit drag
-         DO jj = 2, jpjm1              
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj)
-               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj)
-            END DO
-         END DO
-      ENDIF
-      !       
+                  zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)  &
+                       &                                   + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
+                  zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)  &
+                       &                                   + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
+               END DO
+            END DO
+         ENDIF 
+      ENDIF
+      !
+      !                                   !=  Add atmospheric pressure forcing  =!
+      !                                   !  ----------------------------------  !
       IF( ln_bt_fw ) THEN                        ! Add wind forcing
          DO jj = 2, jpjm1
@@ -604,50 +334,28 @@
       ENDIF  
       !
-      IF( ln_apr_dyn ) THEN                     ! Add atm pressure forcing
-         IF( ln_bt_fw ) THEN
-            DO jj = 2, jpjm1              
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
-                  zv_spg =  grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
-                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
-                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
-               END DO
-            END DO
-         ELSE
-            zztmp = grav * r1_2
-            DO jj = 2, jpjm1              
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zu_spg = zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
-                      &             + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
-                  zv_spg = zztmp * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
-                      &             + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
-                  zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
-                  zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
-               END DO
-            END DO
-         ENDIF 
-      ENDIF
-      !                                   !* Right-Hand-Side of the barotropic ssh equation
-      !                                   ! -----------------------------------------------
-      !                                         ! Surface net water flux and rivers
-      IF (ln_bt_fw) THEN
-         zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
-      ELSE
+      !              !----------------!
+      !              !==  sssh_frc  ==!   Right-Hand-Side of the barotropic ssh equation   (over the FULL domain)
+      !              !----------------!
+      !                                   !=  Net water flux forcing applied to a water column  =!
+      !                                   ! ---------------------------------------------------  !
+      IF (ln_bt_fw) THEN                          ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
+         zssh_frc(:,:) = r1_rau0 * ( emp(:,:)             - rnf(:,:)              + fwfisf(:,:)                  )
+      ELSE                                        ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
          zztmp = r1_rau0 * r1_2
-         zssh_frc(:,:) = zztmp * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
-                &                 + fwfisf(:,:) + fwfisf_b(:,:)                     )
-      ENDIF
-      !
-      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary
+         zssh_frc(:,:) = zztmp * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) + fwfisf(:,:) + fwfisf_b(:,:)  )
+      ENDIF
+      !                                   !=  Add Stokes drift divergence  =!   (if exist)
+      IF( ln_sdw ) THEN                   !  -----------------------------  !
          zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)
       ENDIF
       !
 #if defined key_asminc
-      !                                         ! Include the IAU weighted SSH increment
+      !                                   !=  Add the IAU weighted SSH increment  =!
+      !                                   !  ------------------------------------  !
       IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
          zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:)
       ENDIF
 #endif
-      !                                   !* Fill boundary data arrays for AGRIF
+      !                                   != Fill boundary data arrays for AGRIF
       !                                   ! ------------------------------------
 #if defined key_agrif
@@ -671,5 +379,10 @@
          vb_e   (:,:) = 0._wp
       ENDIF
-
+      !
+      IF( ln_linssh ) THEN    ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0)
+         zhup2_e(:,:) = hu_n(:,:)
+         zhvp2_e(:,:) = hv_n(:,:)
+         zhtp2_e(:,:) = ht_n(:,:)
+      ENDIF
       !
       IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                    
@@ -693,6 +406,4 @@
       ENDIF
       !
-      !
-      !
       ! Initialize sums:
       ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)          
@@ -714,12 +425,13 @@
          !
          l_full_nf_update = jn == icycle   ! false: disable full North fold update (performances) for jn = 1 to icycle-1
-         !                                                !  ------------------
-         !                                                !* Update the forcing (BDY and tides)
-         !                                                !  ------------------
-         ! Update only tidal forcing at open boundaries
+         !
+         !                    !==  Update the forcing ==! (BDY and tides)
+         !
          IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 )
          IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, kt_offset= noffset   )
          !
-         ! Set extrapolation coefficients for predictor step:
+         !                    !==  extrapolation at mid-step  ==!   (jn+1/2)
+         !
+         !                       !* Set extrapolation coefficients for predictor step:
          IF ((jn<3).AND.ll_init) THEN      ! Forward           
            za1 = 1._wp                                          
@@ -731,6 +443,9 @@
            za3 =  0.281105_wp              ! za3 = bet
          ENDIF
-
-         ! Extrapolate barotropic velocities at step jit+0.5:
+         !
+         !                       !* Extrapolate barotropic velocities at mid-step (jn+1/2)
+         !--        m+1/2               m                m-1           m-2       --!
+         !--       u      = (3/2+beta) u   -(1/2+2beta) u      + beta u          --!
+         !-------------------------------------------------------------------------!
          ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)
          va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)
@@ -739,35 +454,11 @@
             !                                             !  ------------------
             ! Extrapolate Sea Level at step jit+0.5:
+            !--         m+1/2                 m                  m-1             m-2       --!
+            !--      ssh      = (3/2+beta) ssh   -(1/2+2beta) ssh      + beta ssh          --!
+            !--------------------------------------------------------------------------------!
             zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
             
-            ! set wetting & drying mask at tracer points for this barotropic sub-step 
-            IF ( ln_wd_dl ) THEN 
-               !
-               IF ( ln_wd_dl_rmp ) THEN 
-                  DO jj = 1, jpj                                 
-                     DO ji = 1, jpi   ! vector opt.  
-                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
-!                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin2 ) THEN 
-                           ztwdmask(ji,jj) = 1._wp
-                        ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
-                           ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1)) ) 
-                        ELSE 
-                           ztwdmask(ji,jj) = 0._wp
-                        END IF
-                     END DO
-                  END DO
-               ELSE
-                  DO jj = 1, jpj                                 
-                     DO ji = 1, jpi   ! vector opt.  
-                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN 
-                           ztwdmask(ji,jj) = 1._wp
-                        ELSE 
-                           ztwdmask(ji,jj) = 0._wp
-                        ENDIF
-                     END DO
-                  END DO
-               ENDIF 
-               !
-            ENDIF 
+            ! set wetting & drying mask at tracer points for this barotropic mid-step
+            IF( ln_wd_dl )   CALL wad_tmsk( zsshp2_e, ztwdmask )
             !
             DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
@@ -786,17 +477,15 @@
             zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:)
             zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:)
-         ELSE
-            zhup2_e(:,:) = hu_n(:,:)
-            zhvp2_e(:,:) = hv_n(:,:)
-            zhtp2_e(:,:) = ht_n(:,:)
-         ENDIF
-         !                                                !* after ssh
-         !                                                !  -----------
-         !
-         ! Enforce volume conservation at open boundaries:
+            !
+         ENDIF
+         !
+         !                    !==  after SSH  ==!   (jn+1)
+         !
+         !                             ! update (ua_e,va_e) to enforce volume conservation at open boundaries
+         !                             ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d
          IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e )
          !
-         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
-         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
+         zhU(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5
+         zhV(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
          !
 #if defined key_agrif
@@ -805,74 +494,58 @@
             IF((nbondi == -1).OR.(nbondi == 2)) THEN
                DO jj = 1, jpj
-                  zwx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)
-                  zwy(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)
+                  zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)
+                  zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)
                END DO
             ENDIF
             IF((nbondi ==  1).OR.(nbondi == 2)) THEN
                DO jj=1,jpj
-                  zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
-                  zwy(nlci-nbghostcells  :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells  :nlci-1,jj)
+                  zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)
+                  zhV(nlci-nbghostcells  :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells  :nlci-1,jj)
                END DO
             ENDIF
             IF((nbondj == -1).OR.(nbondj == 2)) THEN
                DO ji=1,jpi
-                  zwy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)
-                  zwx(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)
+                  zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)
+                  zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)
                END DO
             ENDIF
             IF((nbondj ==  1).OR.(nbondj == 2)) THEN
                DO ji=1,jpi
-                  zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
-                  zwx(ji,nlcj-nbghostcells  :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells  :nlcj-1)
+                  zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)
+                  zhU(ji,nlcj-nbghostcells  :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells  :nlcj-1)
                END DO
             ENDIF
          ENDIF
 #endif
-         IF( ln_wd_il )   CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
-
-         IF ( ln_wd_dl ) THEN 
-            !
-            ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 
-            !
-            DO jj = 1, jpjm1                                 
-               DO ji = 1, jpim1   
-                  IF ( zwx(ji,jj) > 0.0 ) THEN 
-                     zuwdmask(ji, jj) = ztwdmask(ji  ,jj) 
-                  ELSE 
-                     zuwdmask(ji, jj) = ztwdmask(ji+1,jj)  
-                  END IF 
-                  zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj)
-                  un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj)
-
-                  IF ( zwy(ji,jj) > 0.0 ) THEN
-                     zvwdmask(ji, jj) = ztwdmask(ji, jj  )
-                  ELSE 
-                     zvwdmask(ji, jj) = ztwdmask(ji, jj+1)  
-                  END IF 
-                  zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 
-                  vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj)
-               END DO
-            END DO
+         IF( ln_wd_il )   CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rdtbt)
+
+         IF( ln_wd_dl ) THEN           ! un_e and vn_e are set to zero at faces where 
+            !                          ! the direction of the flow is from dry cells
+            CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask )
             !
          ENDIF    
-         
-         ! Sum over sub-time-steps to compute advective velocities
-         za2 = wgtbtp2(jn)
-         un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)
-         vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)
-         
-         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True) 
+         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 
          IF ( ln_wd_dl_bc ) THEN
             zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:)
             zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:)
          END IF 
+         
+         ! Sum over sub-time-steps to compute advective velocities
+         ! 
+         za2 = wgtbtp2(jn)
+         un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:)
+         vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:)
 
          ! Set next sea level:
          DO jj = 2, jpjm1                                 
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
-                  &             + zwy(ji,jj) - zwy(ji,jj-1)   ) * r1_e1e2t(ji,jj)
-            END DO
-         END DO
+               zhdiv(ji,jj) = (   zhU(ji,jj) - zhU(ji-1,jj)   &
+                  &             + zhV(ji,jj) - zhV(ji,jj-1)   ) * r1_e1e2t(ji,jj)
+            END DO
+         END DO
+         !     Compute Sea Level at step jit+1
+         !--           m+1        m                               m+1/2          --!
+         !--        ssh    =  ssh   - delta_t' * [ frc + div( flux      ) ]      --!
+         !-------------------------------------------------------------------------!
          ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:)
          
@@ -899,94 +572,25 @@
             CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )
          ENDIF   
-         !                                 
-         ! Half-step back interpolation of SSH for surface pressure computation:
-         !----------------------------------------------------------------------
-         IF ((jn==1).AND.ll_init) THEN
-           za0=1._wp                        ! Forward-backward
-           za1=0._wp                           
-           za2=0._wp
-           za3=0._wp
-         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
-           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps
-           za1=-0.1666666666666_wp          ! za1 = gam
-           za2= 0.0833333333333_wp          ! za2 = eps
-           za3= 0._wp              
-         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 
-            IF (rn_bt_alpha==0._wp) THEN
-               za0=0.614_wp                 ! za0 = 1/2 +   gam + 2*eps
-               za1=0.285_wp                 ! za1 = 1/2 - 2*gam - 3*eps
-               za2=0.088_wp                 ! za2 = gam
-               za3=0.013_wp                 ! za3 = eps
-            ELSE
-               zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
-               zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
-               za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
-               za1 = 1._wp - za0 - zgamma - zepsilon
-               za2 = zgamma
-               za3 = zepsilon
-            ENDIF 
-         ENDIF
-         !
+         !         
+         ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2
+         !--            m+1/2           m+1              m               m-1              m-2     --!
+         !--        ssh'    =  za0 * ssh     +  za1 * ssh   +  za2 * ssh      +  za3 * ssh        --!
+         !------------------------------------------------------------------------------------------!
+         CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 )   ! coeficients of the interpolation
          zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   &
             &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
-         
-         IF( ln_wd_il ) THEN                   ! Calculating and applying W/D gravity filters
-           DO jj = 2, jpjm1
-              DO ji = 2, jpim1 
-                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
-                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) .AND.            &
-                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) &
-                     &                                                             > rn_wdmin1 + rn_wdmin2
-                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji+1,jj))  > 1.E-12 ).AND.( &
-                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                &
-                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
-   
-                IF(ll_tmp1) THEN
-                  zcpx(ji,jj) = 1.0_wp
-                ELSE IF(ll_tmp2) THEN
-                  ! no worries about  zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj) = 0, it won't happen ! here
-                  zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +     ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
-                              &    / (zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj)) )
-                ELSE
-                  zcpx(ji,jj) = 0._wp
-                ENDIF
-                !
-                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
-                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) .AND.            &
-                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) &
-                     &                                                             > rn_wdmin1 + rn_wdmin2
-                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji,jj+1))  > 1.E-12 ).AND.( &
-                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
-                     &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
-   
-                IF(ll_tmp1) THEN
-                  zcpy(ji,jj) = 1.0_wp
-                ELSEIF(ll_tmp2) THEN
-                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here
-                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +     ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
-                              &    / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  )) )
-                ELSE
-                  zcpy(ji,jj) = 0._wp
-                ENDIF
-              END DO
-           END DO
-         ENDIF
-         !
-         ! Compute associated depths at U and V points:
-         IF( .NOT.ln_linssh  .AND. .NOT.ln_dynadv_vec ) THEN     !* Vector form
-            !                                        
-            DO jj = 2, jpjm1                            
-               DO ji = 2, jpim1
-                  zx1 = r1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    &
-                     &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
-                     &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
-                  zy1 = r1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  &
-                     &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
-                     &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
-                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
-                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
-               END DO
-            END DO
-            !
+         !
+         !                             ! Surface pressure gradient
+         zldg = ( 1._wp - rn_scal_load ) * grav    ! local factor
+         DO jj = 2, jpjm1                            
+            DO ji = 2, jpim1
+               zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
+               zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
+            END DO
+         END DO
+         IF( ln_wd_il ) THEN        ! W/D : gravity filters applied on pressure gradient
+            CALL wad_spg( zsshp2_e, zcpx, zcpy )   ! Calculating W/D gravity filters
+            zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
+            zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
          ENDIF
          !
@@ -994,65 +598,6 @@
          ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
          ! at each time step. We however keep them constant here for optimization.
-         ! Recall that zwx and zwy arrays hold fluxes at this stage:
-         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5
-         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
-         !
-         SELECT CASE( nvor_scheme )
-         CASE( np_ENT )             ! energy conserving scheme (t-point)
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1   ! vector opt.
-
-               z1_hu = ssumask(ji,jj) / ( zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) )
-               z1_hv = ssvmask(ji,jj) / ( zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) )
-            
-               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                   &
-                  &               * (  e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) )   &
-                  &                  + e1e2t(ji  ,jj)*zhtp2_e(ji  ,jj)*ff_t(ji  ,jj) * ( va_e(ji  ,jj) + va_e(ji  ,jj-1) )   )
-                  !
-               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
-                  &               * (  e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) )   & 
-                  &                  + e1e2t(ji,jj  )*zhtp2_e(ji,jj  )*ff_t(ji,jj  ) * ( ua_e(ji,jj  ) + ua_e(ji-1,jj  ) )   ) 
-            END DO  
-         END DO  
-         !         
-         CASE( np_ENE, np_MIX )     ! energy conserving scheme (f-point)
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj)
-                  zy2 = ( zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
-                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
-                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
-                  zu_trd(ji,jj) = r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-                  zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
-               END DO
-            END DO
-            !
-         CASE( np_ENS )             ! enstrophy conserving scheme (f-point)
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
-                   &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
-                  zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
-                   &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
-                  zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
-                  zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
-               END DO
-            END DO
-            !
-         CASE( np_EET , np_EEN )   ! energy & enstrophy scheme (using e3t or e3f)
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  )  &
-                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  )  &
-                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1)  & 
-                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1)  )
-                  zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1)  & 
-                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1)  &
-                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  )  & 
-                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  )  )
-               END DO
-            END DO
-            ! 
-         END SELECT
+         ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
+         CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV,    zu_trd, zv_trd   )
          !
          ! Add tidal astronomical forcing if defined
@@ -1060,8 +605,6 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
-                  zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
-                  zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg
-                  zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg
+                  zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
+                  zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
                END DO
             END DO
@@ -1077,36 +620,22 @@
                END DO
             END DO
-         ENDIF 
-         !
-         ! Surface pressure trend:
-         IF( ln_wd_il ) THEN
-           DO jj = 2, jpjm1
-              DO ji = 2, jpim1 
-                 ! Add surface pressure gradient
-                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
-                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
-                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj) 
-                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj)
-              END DO
-           END DO
-         ELSE
-           DO jj = 2, jpjm1
-              DO ji = fs_2, fs_jpim1   ! vector opt.
-                 ! Add surface pressure gradient
-                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
-                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
-                 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg
-                 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg
-              END DO
-           END DO
-         END IF
-
+         ENDIF
          !
          ! Set next velocities:
+         !     Compute barotropic speeds at step jit+1    (h : total height of the water colomn)
+         !--                              VECTOR FORM
+         !--     m+1                    m               /                                     m+1/2           \     --!
+         !--    u     =                u   + delta_t' * \ grad_x( ssh') -         f * k vect u      +     frc /     --!
+         !--                                                                                                        --!
+         !--                             FLUX FORM                                                                  --!
+         !--   m+1    ___1___   /  m    m               /                  m+1/2              m+1/2    n      \  \  --!
+         !--  u    =     m+1   |  h  * u   + delta_t' * \ grad_x( ssh') - h     * f * k vect u      + h * frc /   | --!
+         !--            h       \                                                                                /  --!
+         !------------------------------------------------------------------------------------------------------------!
          IF( ln_dynadv_vec .OR. ln_linssh ) THEN      !* Vector form
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ua_e(ji,jj) = (                                 un_e(ji,jj)   & 
-                            &     + rdtbt * (                      zwx(ji,jj)   &
+                            &     + rdtbt * (                   zu_spg(ji,jj)   &
                             &                                 + zu_trd(ji,jj)   &
                             &                                 + zu_frc(ji,jj) ) & 
@@ -1114,9 +643,8 @@
 
                   va_e(ji,jj) = (                                 vn_e(ji,jj)   &
-                            &     + rdtbt * (                      zwy(ji,jj)   &
+                            &     + rdtbt * (                   zv_spg(ji,jj)   &
                             &                                 + zv_trd(ji,jj)   &
                             &                                 + zv_frc(ji,jj) ) &
                             &   ) * ssvmask(ji,jj)
- 
                END DO
             END DO
@@ -1125,22 +653,22 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-
-                  zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
-                  zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
-
-                  zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))
-                  zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))
-
-                  ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   & 
-                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
-                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
-                            &               +    hu_n(ji,jj)  * zu_frc(ji,jj) ) &
-                            &   ) * zhura
-
-                  va_e(ji,jj) = (                hv_e(ji,jj)  *   vn_e(ji,jj)   &
-                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
-                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
-                            &               +    hv_n(ji,jj)  * zv_frc(ji,jj) ) &
-                            &   ) * zhvra
+                  !                    ! backward extrapolated depth used in spg terms at jn+1/2
+                  zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    &
+                       &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
+                  zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    &
+                       &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj)
+                  !                    ! inverse depth at jn+1
+                  z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) )
+                  z1_hv = ssvmask(ji,jj) / ( hu_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) )
+                  !
+                  ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      & 
+                       &            + rdtbt * (  zhu_bck        * zu_spg (ji,jj)  &   !
+                       &                       + zhup2_e(ji,jj) * zu_trd (ji,jj)  &   !
+                       &                       +  hu_n  (ji,jj) * zu_frc (ji,jj)  )   ) * z1_hu
+                  !
+                  va_e(ji,jj) = (               hv_e  (ji,jj) *   vn_e (ji,jj)      &
+                       &            + rdtbt * (  zhv_bck        * zv_spg (ji,jj)  &   !
+                       &                       + zhvp2_e(ji,jj) * zv_trd (ji,jj)  &   !
+                       &                       +  hv_n  (ji,jj) * zv_frc (ji,jj)  )   ) * z1_hu
                END DO
             END DO
@@ -1213,20 +741,26 @@
       ! Set advection velocity correction:
       IF (ln_bt_fw) THEN
-         zwx(:,:) = un_adv(:,:)
-         zwy(:,:) = vn_adv(:,:)
          IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
-            un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
-            vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
-            !
-            ! Update corrective fluxes for next time step:
-            un_bf(:,:)  = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))
-            vn_bf(:,:)  = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zun_save = un_adv(ji,jj)
+                  zvn_save = vn_adv(ji,jj)
+                  !                          ! apply the previously computed correction 
+                  un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) )
+                  vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) )
+                  !                          ! Update corrective fluxes for next time step
+                  un_bf(ji,jj)  = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )
+                  vn_bf(ji,jj)  = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )
+                  !                          ! Save integrated transport for next computation
+                  ub2_b(ji,jj) = zun_save
+                  vb2_b(ji,jj) = zvn_save
+               END DO
+            END DO
          ELSE
-            un_bf(:,:) = 0._wp
-            vn_bf(:,:) = 0._wp 
-         END IF         
-         ! Save integrated transport for next computation
-         ub2_b(:,:) = zwx(:,:)
-         vb2_b(:,:) = zwy(:,:)
+            un_bf(:,:) = 0._wp            ! corrective fluxes for next time step set to zero
+            vn_bf(:,:) = 0._wp
+            ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation
+            vb2_b(:,:) = vn_adv(:,:)
+         END IF
       ENDIF
 
@@ -1473,4 +1007,5 @@
       REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
       REAL(wp), DIMENSION(jpi,jpj) ::   zcu
+      INTEGER  :: inum
       !!----------------------------------------------------------------------
       !
@@ -1579,4 +1114,497 @@
    END SUBROUTINE dyn_spg_ts_init
 
+   
+   SUBROUTINE dyn_cor_2d_init
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE dyn_cor_2d_init  ***
+      !!
+      !! ** Purpose : Set time splitting options
+      !! Set arrays to remove/compute coriolis trend.
+      !! Do it once during initialization if volume is fixed, else at each long time step.
+      !! Note that these arrays are also used during barotropic loop. These are however frozen
+      !! although they should be updated in the variable volume case. Not a big approximation.
+      !! To remove this approximation, copy lines below inside barotropic loop
+      !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step
+      !!
+      !! Compute zwz = f / ( height of the water colomn )
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji ,jj, jk              ! dummy loop indices
+      REAL(wp) ::   z1_ht
+      REAL(wp), DIMENSION(jpi,jpj) :: zhf
+      !!----------------------------------------------------------------------
+      !
+      !
+      SELECT CASE( nvor_scheme )
+      CASE( np_EEN )                != EEN scheme using e3f (energy & enstrophy scheme)
+         SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point
+         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    &
+                       &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp  
+                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
+               END DO
+            END DO
+         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      &
+                       &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   &
+                       &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      &
+                       &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   )
+                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
+               END DO
+            END DO
+         END SELECT
+         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
+         !
+         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
+         DO jj = 2, jpj
+            DO ji = 2, jpi
+               ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
+               ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
+               ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
+               ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
+            END DO
+         END DO
+         !
+      CASE( np_EET )                  != EEN scheme using e3t (energy conserving scheme)
+         ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
+         DO jj = 2, jpj
+            DO ji = 2, jpi
+               z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
+               ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
+               ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht
+               ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
+               ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht
+            END DO
+         END DO
+         !
+      CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT !
+         !
+         zwz(:,:) = 0._wp
+         zhf(:,:) = 0._wp
+         
+         !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed 
+!!gm    A priori a better value should be something like :
+!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 
+!!gm                     divided by the sum of the corresponding mask 
+!!gm 
+!!            
+         IF( .NOT.ln_sco ) THEN
+  
+   !!gm  agree the JC comment  : this should be done in a much clear way
+  
+   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
+   !     Set it to zero for the time being 
+   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
+   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
+   !              ENDIF
+   !              zhf(:,:) = gdepw_0(:,:,jk+1)
+            !
+         ELSE
+            !
+            !zhf(:,:) = hbatf(:,:)
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          &
+                       &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      &
+                       &     / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          &
+                       &            + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  )
+               END DO
+            END DO
+         ENDIF
+         !
+         DO jj = 1, jpjm1
+            zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
+         END DO
+         !
+         DO jk = 1, jpkm1
+            DO jj = 1, jpjm1
+               zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)
+            END DO
+         END DO
+         CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )
+         ! JC: TBC. hf should be greater than 0 
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj)
+            END DO
+         END DO
+         zwz(:,:) = ff_f(:,:) * zwz(:,:)
+      END SELECT
+      
+   END SUBROUTINE dyn_cor_2d_init
+
+
+
+   SUBROUTINE dyn_cor_2d( ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV,    zu_trd, zv_trd   )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE dyn_cor_2d  ***
+      !!
+      !! ** Purpose : Compute u and v coriolis trends
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji ,jj               ! dummy loop indices
+      REAL(wp) ::   zx1, zx2, zy1, zy2   !   -      -
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd
+      !!----------------------------------------------------------------------
+      SELECT CASE( nvor_scheme )
+      CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj)                    &
+                  &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   &
+                  &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   )
+                  !
+               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj)                    &
+                  &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   & 
+                  &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   ) 
+            END DO  
+         END DO  
+         !         
+      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
+               zy2 = ( zhV(ji,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
+               zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
+               zx2 = ( zhU(ji  ,jj) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
+               ! energy conserving formulation for planetary vorticity term
+               zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
+               zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
+            END DO
+         END DO
+         !
+      CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zy1 =   r1_8 * ( zhV(ji  ,jj-1) + zhV(ji+1,jj-1) &
+                 &            + zhV(ji  ,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
+               zx1 = - r1_8 * ( zhU(ji-1,jj  ) + zhU(ji-1,jj+1) &
+                 &            + zhU(ji  ,jj  ) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
+               zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
+               zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
+            END DO
+         END DO
+         !
+      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
+                &                                         + ftnw(ji+1,jj) * zhV(ji+1,jj  ) &
+                &                                         + ftse(ji,jj  ) * zhV(ji  ,jj-1) &
+                &                                         + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
+               zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
+                &                                         + ftse(ji,jj+1) * zhU(ji  ,jj+1) &
+                &                                         + ftnw(ji,jj  ) * zhU(ji-1,jj  ) &
+                &                                         + ftne(ji,jj  ) * zhU(ji  ,jj  ) )
+            END DO
+         END DO
+         !
+      END SELECT
+      !
+   END SUBROUTINE dyn_cor_2D
+
+
+   SUBROUTINE wad_tmsk( pssh, ptmsk )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE wad_lmt  ***
+      !!                    
+      !! ** Purpose :   set wetting & drying mask at tracer points 
+      !!              for the current barotropic sub-step 
+      !!
+      !! ** Method  :   ??? 
+      !!
+      !! ** Action  :  ptmsk : wetting & drying t-mask
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pssh    !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   ptmsk   !
+      !
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_wd_dl_rmp ) THEN     
+         DO jj = 1, jpj
+            DO ji = 1, jpi                    
+               IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
+                  !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN 
+                  ptmsk(ji,jj) = 1._wp
+               ELSEIF( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin1 ) THEN
+                  ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1) )
+               ELSE 
+                  ptmsk(ji,jj) = 0._wp
+               ENDIF
+            END DO
+         END DO
+      ELSE  
+         DO jj = 1, jpj
+            DO ji = 1, jpi                              
+               IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
+               ELSE                                                 ;   ptmsk(ji,jj) = 0._wp
+               ENDIF
+            END DO
+         END DO
+      ENDIF
+      !
+   END SUBROUTINE wad_tmsk
+
+
+   SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE wad_lmt  ***
+      !!                    
+      !! ** Purpose :   set wetting & drying mask at tracer points 
+      !!              for the current barotropic sub-step 
+      !!
+      !! ** Method  :   ??? 
+      !!
+      !! ** Action  :  ptmsk : wetting & drying t-mask
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pTmsk              ! W & D t-mask
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   phU, phV, pu, pv   ! ocean velocities and transports
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pUmsk, pVmsk       ! W & D u- and v-mask
+      !
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpim1   ! not jpi-column
+            IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
+            ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)  
+            ENDIF
+            phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
+            pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
+         END DO
+      END DO
+      !
+      DO jj = 1, jpjm1   ! not jpj-row
+         DO ji = 1, jpi
+            IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
+            ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)  
+            ENDIF
+            phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 
+            pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
+         END DO
+      END DO
+      !
+   END SUBROUTINE wad_Umsk
+
+
+   SUBROUTINE wad_spg( sshn, zcpx, zcpy )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE  wad_sp  ***
+      !!
+      !! ** Purpose : 
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji ,jj               ! dummy loop indices
+      LOGICAL  ::   ll_tmp1, ll_tmp2
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: sshn
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
+      !!----------------------------------------------------------------------
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1 
+            ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
+                 &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
+                 &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
+                 &                                                         > rn_wdmin1 + rn_wdmin2
+            ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( &
+                 &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
+                 &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
+            IF(ll_tmp1) THEN
+               zcpx(ji,jj) = 1.0_wp
+            ELSEIF(ll_tmp2) THEN
+               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
+                    &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
+               zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
+            ELSE
+               zcpx(ji,jj) = 0._wp
+            ENDIF
+            !
+            ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                &
+                 &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
+                 &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
+                 &                                                       > rn_wdmin1 + rn_wdmin2
+            ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( &
+                 &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
+                 &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
+            
+            IF(ll_tmp1) THEN
+               zcpy(ji,jj) = 1.0_wp
+            ELSE IF(ll_tmp2) THEN
+               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
+                    &             / (sshn(ji,jj+1) - sshn(ji,jj  )) )
+               zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
+            ELSE
+               zcpy(ji,jj) = 0._wp
+            ENDIF
+         END DO
+      END DO
+            
+   END SUBROUTINE wad_spg
+     
+
+
+   SUBROUTINE drg_init( pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE drg_init  ***
+      !!                    
+      !! ** Purpose : - add the baroclinic top/bottom drag contribution to 
+      !!              the baroclinic part of the barotropic RHS
+      !!              - compute the barotropic drag coefficients
+      !!
+      !! ** Method  :   computation done over the INNER domain only 
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pCdU_u , pCdU_v    ! barotropic drag coefficients
+      !
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      INTEGER  ::   ikbu, ikbv, iktu, iktv
+      REAL(wp) ::   zztmp
+      REAL(wp), DIMENSION(jpi,jpj) ::   zu_i, zv_i
+      !!----------------------------------------------------------------------
+      !
+      !                    !==  Set the barotropic drag coef.  ==!
+      !
+      IF( ln_isfcav ) THEN          ! top+bottom friction (ocean cavities)
+         
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1     ! INNER domain
+               pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
+               pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
+            END DO
+         END DO
+      ELSE                          ! bottom friction only
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1  ! INNER domain
+               pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
+               pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
+            END DO
+         END DO
+      ENDIF
+      !
+      !                    !==  BOTTOM stress contribution from baroclinic velocities  ==!
+      !
+      IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
+         
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1  ! INNER domain
+               ikbu = mbku(ji,jj)       
+               ikbv = mbkv(ji,jj)    
+               zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj)
+               zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
+            END DO
+         END DO
+      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
+         
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1   ! INNER domain
+               ikbu = mbku(ji,jj)       
+               ikbv = mbkv(ji,jj)    
+               zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj)
+               zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      !
+      IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
+         zztmp = -1._wp / rdtbt
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1    ! INNER domain
+               pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
+                    &                              r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
+               pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 & 
+                    &                              r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  )
+            END DO
+         END DO
+      ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
+         
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1    ! INNER domain
+               pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj)
+               pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj)
+            END DO
+         END DO
+      END IF
+      !
+      !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case)
+      !
+      IF( ln_isfcav ) THEN
+         !
+         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
+            
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1   ! INNER domain
+                  iktu = miku(ji,jj)
+                  iktv = mikv(ji,jj)
+                  zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj)
+                  zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
+               END DO
+            END DO
+         ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
+            
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1      ! INNER domain
+                  iktu = miku(ji,jj)
+                  iktv = mikv(ji,jj)
+                  zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj)
+                  zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
+               END DO
+            END DO
+         ENDIF
+         !
+         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
+         
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1    ! INNER domain
+               pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj)
+               pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj)
+            END DO
+         END DO
+         !
+      ENDIF
+      !
+   END SUBROUTINE drg_init
+
+   SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in
+      &                      za0, za1, za2, za3 )   ! ==> out
+      !!----------------------------------------------------------------------
+      INTEGER ,INTENT(in   ) ::   jn                   ! index of sub time step
+      LOGICAL ,INTENT(in   ) ::   ll_init              !
+      REAL(wp),INTENT(  out) ::   za0, za1, za2, za3   ! Half-step back interpolation coefficient
+      !
+      REAL(wp) ::   zepsilon, zgamma                   !   -      -
+      !!----------------------------------------------------------------------
+      !                             ! set Half-step back interpolation coefficient
+      IF    ( jn==1 .AND. ll_init ) THEN   !* Forward-backward
+         za0 = 1._wp                        
+         za1 = 0._wp                           
+         za2 = 0._wp
+         za3 = 0._wp
+      ELSEIF( jn==2 .AND. ll_init ) THEN   !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
+         za0 = 1.0833333333333_wp                 ! za0 = 1-gam-eps
+         za1 =-0.1666666666666_wp                 ! za1 = gam
+         za2 = 0.0833333333333_wp                 ! za2 = eps
+         za3 = 0._wp              
+      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 
+         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion  
+            za0 = 0.614_wp                        ! za0 = 1/2 +   gam + 2*eps
+            za1 = 0.285_wp                        ! za1 = 1/2 - 2*gam - 3*eps
+            za2 = 0.088_wp                        ! za2 = gam
+            za3 = 0.013_wp                        ! za3 = eps
+         ELSE                                 ! no time diffusion
+            zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
+            zgamma   = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
+            za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
+            za1 = 1._wp - za0 - zgamma - zepsilon
+            za2 = zgamma
+            za3 = zepsilon
+         ENDIF 
+      ENDIF
+   END SUBROUTINE ts_bck_interp
+
+
    !!======================================================================
 END MODULE dynspg_ts
