Index: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dom_oce.F90	(revision 10023)
+++ NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dom_oce.F90	(revision 10030)
@@ -135,6 +135,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0                                !: uw-vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0                                !: vw-vert. scale factor [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
 
@@ -270,7 +270,5 @@
          &      e3t(jpi,jpj,jpk,Nt) , e3u(jpi,jpj,jpk,Nt) , e3v(jpi,jpj,jpk,Nt) ,                                          &
          !                                                          !
-         &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) ,         &
-         &      e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,         &               
-         &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) )                       
+         &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) ,     STAT=ierr(5) )                       
          !
       ALLOCATE( ht_0(jpi,jpj)    , hu_0(jpi,jpj)   , hv_0(jpi,jpj)    , hf_0(jpi,jpj)    ,                     &
Index: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90	(revision 10023)
+++ NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90	(revision 10030)
@@ -170,11 +170,9 @@
                                  gde3w_n = gde3w_0   !        ---          ! grid-points
          !                   !                       !                     !                                              
-          e3t_b  =   e3t_0   ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    !  scale 
-          e3u_b  =   e3u_0   ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    ! factors
-          e3v_b  =   e3v_0   ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
+           e3t_b =   e3t_0   ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    !  scale 
+           e3u_b =   e3u_0   ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    ! factors
+           e3v_b =   e3v_0   ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
                                    e3f_n =   e3f_0   !        ---          !
-          e3w_b  =   e3w_0   ;     e3w_n =   e3w_0   !        ---          !
-         e3uw_b  =  e3uw_0   ;    e3uw_n =  e3uw_0   !        ---          !
-         e3vw_b  =  e3vw_0   ;    e3vw_n =  e3vw_0   !        ---          !
+           e3w_b =   e3w_0   ;     e3w_n =   e3w_0   !        ---          !
          !                   !
                                     ht_n =    ht_0   !                     ! water column
Index: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90
===================================================================
--- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90	(revision 10023)
+++ NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90	(revision 10030)
@@ -144,6 +144,5 @@
       CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv 
       !                                !      e3t, e3u , e3v              (from 1 to jpkm1)
-      !                                !      e3w, e3uw, e3vw             (from 1 to jpk  )
-      !                                !      gdept, gdepw                (from 1 to jpk  )
+      !                                !      e3w, gdept, gdepw           (from 1 to jpk  )
       !
       !                                ! set jpk level one to the e3._0 values
@@ -153,10 +152,9 @@
       CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
       !                                !      e3t, e3u , e3v, e3f         (from 1 to jpkm1)
-      !                                !      e3w, e3uw, e3vw             (from 1 to jpk  )
-      !                                !      gdept, gdepw, gde3w         (from 1 to jpk  )
+      !                                !      e3w, gdept, gdepw, gde3w    (from 1 to jpk  )
       !
       !                                ! set one for all last level to the e3._0 value
-      e3t_n(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_n(:,:,jpk) =  e3u_0(:,:,jpk)  ;   e3v_n(:,:,jpk) =  e3v_0(:,:,jpk)
-      e3f_n(:,:,jpk) = e3f_0(:,:,jpk)
+      e3t_n(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_n(:,:,jpk) =  e3u_0(:,:,jpk)
+      e3f_n(:,:,jpk) = e3f_0(:,:,jpk)  ;   e3v_n(:,:,jpk) =  e3v_0(:,:,jpk)
       !
       !                          !* AFTER fields : (last level for OPA, 3D required for AGRIF initialisation)
@@ -215,6 +213,4 @@
       !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
       !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
-      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
-      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
       !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
       !
@@ -332,6 +328,5 @@
       ! - JC - hu_b, hv_b, hur_b, hvr_b also
       !
-      ! - GM - to be updated :   e3f_n,  e3w_n , e3uw_n , e3vw_n 
-      !                                  e3w_b , e3uw_b , e3vw_b
+      ! - GM - to be updated :   e3f_n   ,  e3w_n  , e3w_b
       !                          gdept_n , gdepw_n , gde3w_n
       !                          ht_n
@@ -356,19 +351,10 @@
       !
       !                    !==  before  ==!
-      !                                            !* ssh at u- and v-points)
-      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
-         zsshu_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
-         zsshv_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
-      END DO             ;   END DO      
-      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp )
-      !
-      !                                            !*  e3w_b , e3uw_b , e3vw_b
+      !                                            !*  e3w_b
       zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)     ! w-point
-      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! uw-point
-      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! vw-point
-      DO jk = 1, jpkm1
-          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
-         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
-         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
+      !
+      e3w_b(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
+      DO jk = 2, jpk
+         e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
       END DO
       ! 
@@ -388,25 +374,22 @@
       !      
       !                    !==   now    ==!
-      !                                            !* ssh at u- and v-points)
-      DO jj = 1, jpjm1   ;   DO ji = 1, jpim1            ! start from 1 for f-point
-         zsshu_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
-         zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
-         zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
-            &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
-      END DO             ;   END DO      
-      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp )      
-      !
-      !                                            !* e3w_n , e3uw_n , e3vw_n, e3f_n 
-      zssht_h(:,:) = ssh    (:,:,Nnn) * r1_ht_0(:,:)     ! t- & w-point
-      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! uw-point
-      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! vw-point
-      zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)     ! f-point
+      !                                            !* ssh at f-points
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1            ! start from 1 for f-point
+            zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
+               &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
+         END DO
+      END DO      
+      CALL lbc_lnk( zsshf_h(:,:),'F', 1._wp )      
+      !
+      !                                            !* e3f_n 
+      zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:)     ! f-point
+      !
       DO jk = 1, jpkm1
-         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) )
-         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) )
-          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) *  fmask(:,:,jk) )
+          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
       END DO      
-      ! 
-      zssht_h(:,:) = 1._wp + zssht_h(:,:)          !* gdept , gdepw , gde3w
+      !
+      !                                            !* gdept_n , gdepw_n , gde3w_n
+      zssht_h(:,:) = 1._wp + ssh(:,:,Nnn) * r1_ht_0(:,:) 
       !
       IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
@@ -633,6 +616,4 @@
       !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
       !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
-      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
-      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
       !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
       !
@@ -658,7 +639,6 @@
       !
       !                             !==  ssh at u- and v-points  ==!
-      !
-      DO jj = 1, jpjm1                    ! start from 1 due to f-point
-         DO ji = 1, jpim1
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1                 ! start from 1 due to f-point
             zsshu_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
             zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
@@ -670,5 +650,4 @@
       !
       !                             !==  ht, hu and hv  == !   (and their inverse)
-      !
       ht_n   (:,:) = ht_0(:,:) +  ssh   (:,:,Nnn)
       hu_n   (:,:) = hu_0(:,:) + zsshu_h(:,:)
@@ -678,5 +657,4 @@
       !      
       !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==!
-      !
       zssht_h(:,:) =  ssh   (:,:,Nnn) * r1_ht_0(:,:)
       zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
@@ -684,6 +662,5 @@
       zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)
       !
-      !                             !==  e3t  ,  e3u  ,  e3v  ,  e3f  ==!
-      !      
+      !                             !==  e3t  ,  e3u  ,  e3v  ,  e3f  ,  e3w  ==!
       DO jk = 1, jpkm1
           e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
@@ -692,14 +669,7 @@
           e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
       END DO
-      !      
-      !                             !==  e3w  ,  e3uw  ,  e3vw  ==!
-      !
-       e3w_n(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
-      e3uw_n(:,:,1) = e3uw_0(:,:,1) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,1) )
-      e3vw_n(:,:,1) = e3vw_0(:,:,1) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,1) )
+      e3w_n(:,:,1) = e3w_0(:,:,1) * (  1._wp + zssht_h(:,:) *  tmask(:,:,1)  )
       DO jk = 2, jpk
-          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX(  tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
-         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * MAX( wumask(:,:,jk-1) , wumask(:,:,jk) ) )
-         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * MAX( wvmask(:,:,jk-1) , wvmask(:,:,jk) ) )
+         e3w_n(:,:,jk) = e3w_0(:,:,jk) * (  1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) , tmask(:,:,jk) )  )
       END DO
       !
@@ -733,6 +703,4 @@
       !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
       !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
-      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
-      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
       !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
       !
@@ -778,6 +746,5 @@
       zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
       !
-      !                             !==  e3t  ,  e3u  ,  e3v  ==!
-      !      
+      !                             !==  e3t  ,  e3u  ,  e3v ,  e3w  ==!
       DO jk = 1, jpkm1
           e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
@@ -785,14 +752,7 @@
           e3v_b(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
       END DO
-      !      
-      !                             !==  e3w  ,  e3uw  ,  e3vw  ==!
-      !
-       e3w_b(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
-      e3uw_b(:,:,1) = e3uw_0(:,:,1) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,1) )
-      e3vw_b(:,:,1) = e3vw_0(:,:,1) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,1) )
+      e3w_b(:,:,1) = e3w_0(:,:,1) * (  1._wp + zssht_h(:,:) *  tmask(:,:,1)  )
       DO jk = 2, jpk
-          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX(  tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
-         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * MAX( wumask(:,:,jk-1) , wumask(:,:,jk) ) )
-         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * MAX( wvmask(:,:,jk-1) , wvmask(:,:,jk) ) )
+         e3w_b(:,:,jk) = e3w_0(:,:,jk) * (  1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) , tmask(:,:,jk) )  )
       END DO
       !   
Index: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl_RK3.F90
===================================================================
--- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl_RK3.F90	(revision 10023)
+++ NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl_RK3.F90	(revision 10030)
@@ -144,6 +144,5 @@
       CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv 
       !                                !      e3t, e3u , e3v              (from 1 to jpkm1)
-      !                                !      e3w, e3uw, e3vw             (from 1 to jpk  )
-      !                                !      gdept, gdepw                (from 1 to jpk  )
+      !                                !      e3w, gdept, gdepw           (from 1 to jpk  )
       !
       !                                ! set jpk level one to the e3._0 values
@@ -153,6 +152,5 @@
       CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
       !                                !      e3t, e3u , e3v, e3f         (from 1 to jpkm1)
-      !                                !      e3w, e3uw, e3vw             (from 1 to jpk  )
-      !                                !      gdept, gdepw, gde3w         (from 1 to jpk  )
+      !                                !      e3w, gdept, gdepw, gde3w    (from 1 to jpk  )
       !
       !                                ! set one for all last level to the e3._0 value
@@ -215,6 +213,6 @@
       !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
       !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
-      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
-      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+      !    e3uw_0                                !: uw-vert. scale factor [m]
+      !    e3vw_0                                !: vw-vert. scale factor [m]
       !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
       !
@@ -332,6 +330,5 @@
       ! - JC - hu_b, hv_b, hur_b, hvr_b also
       !
-      ! - GM - to be updated :   e3f_n,  e3w_n , e3uw_n , e3vw_n 
-      !                                  e3w_b , e3uw_b , e3vw_b
+      ! - GM - to be updated :   e3f_n   ,  e3w_n  , e3w_b
       !                          gdept_n , gdepw_n , gde3w_n
       !                          ht_n
@@ -356,19 +353,10 @@
       !
       !                    !==  before  ==!
-      !                                            !* ssh at u- and v-points)
-      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
-         zsshu_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
-         zsshv_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
-      END DO             ;   END DO      
-      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp )
-      !
-      !                                            !*  e3w_b , e3uw_b , e3vw_b
+      !                                            !*  e3w_b
       zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)     ! w-point
-      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! uw-point
-      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! vw-point
-      DO jk = 1, jpkm1
-          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
-         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
-         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
+      !
+      e3w_b(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
+      DO jk = 2, jpk
+         e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
       END DO
       ! 
@@ -388,25 +376,22 @@
       !      
       !                    !==   now    ==!
-      !                                            !* ssh at u- and v-points)
-      DO jj = 1, jpjm1   ;   DO ji = 1, jpim1            ! start from 1 for f-point
-         zsshu_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
-         zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
-         zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
-            &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
-      END DO             ;   END DO      
-      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp )      
-      !
-      !                                            !* e3w_n , e3uw_n , e3vw_n, e3f_n 
-      zssht_h(:,:) = ssh    (:,:,Nnn) * r1_ht_0(:,:)     ! t- & w-point
-      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! uw-point
-      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! vw-point
-      zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)     ! f-point
+      !                                            !* ssh at f-points
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1            ! start from 1 for f-point
+            zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
+               &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
+         END DO
+      END DO      
+      CALL lbc_lnk( zsshf_h(:,:),'F', 1._wp )      
+      !
+      !                                            !* e3f_n 
+      zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:)     ! f-point
+      !
       DO jk = 1, jpkm1
-         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) )
-         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) )
-          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) *  fmask(:,:,jk) )
+          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
       END DO      
-      ! 
-      zssht_h(:,:) = 1._wp + zssht_h(:,:)          !* gdept , gdepw , gde3w
+      !
+      !                                            !* gdept_n , gdepw_n , gde3w_n
+      zssht_h(:,:) = 1._wp + ssh(:,:,Nnn) * r1_ht_0(:,:) 
       !
       IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness 
@@ -633,6 +618,6 @@
       !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
       !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
-      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
-      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+      !    e3uw_0                                !: uw-vert. scale factor [m]
+      !    e3vw_0                                !: vw-vert. scale factor [m]
       !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
       !
@@ -658,7 +643,6 @@
       !
       !                             !==  ssh at u- and v-points  ==!
-      !
-      DO jj = 1, jpjm1                    ! start from 1 due to f-point
-         DO ji = 1, jpim1
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1                 ! start from 1 due to f-point
             zsshu_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
             zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
@@ -670,5 +654,4 @@
       !
       !                             !==  ht, hu and hv  == !   (and their inverse)
-      !
       ht_n   (:,:) = ht_0(:,:) +  ssh   (:,:,Nnn)
       hu_n   (:,:) = hu_0(:,:) + zsshu_h(:,:)
@@ -678,5 +661,4 @@
       !      
       !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==!
-      !
       zssht_h(:,:) =  ssh   (:,:,Nnn) * r1_ht_0(:,:)
       zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
@@ -684,6 +666,5 @@
       zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)
       !
-      !                             !==  e3t  ,  e3u  ,  e3v  ,  e3f  ==!
-      !      
+      !                             !==  e3t  ,  e3u  ,  e3v  ,  e3f  ,  e3w  ==!
       DO jk = 1, jpkm1
           e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
@@ -692,14 +673,7 @@
           e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
       END DO
-      !      
-      !                             !==  e3w  ,  e3uw  ,  e3vw  ==!
-      !
-       e3w_n(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
-      e3uw_n(:,:,1) = e3uw_0(:,:,1) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,1) )
-      e3vw_n(:,:,1) = e3vw_0(:,:,1) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,1) )
+      e3w_n(:,:,1) = e3w_0(:,:,1) * (  1._wp + zssht_h(:,:) *  tmask(:,:,1)  )
       DO jk = 2, jpk
-          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX(  tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
-         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * MAX( wumask(:,:,jk-1) , wumask(:,:,jk) ) )
-         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * MAX( wvmask(:,:,jk-1) , wvmask(:,:,jk) ) )
+         e3w_n(:,:,jk) = e3w_0(:,:,jk) * (  1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) , tmask(:,:,jk) )  )
       END DO
       !
@@ -733,6 +707,6 @@
       !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
       !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
-      !    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
-      !    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+      !    e3uw_0                                !: uw-vert. scale factor [m]
+      !    e3vw_0                                !: vw-vert. scale factor [m]
       !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
       !
@@ -778,6 +752,5 @@
       zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
       !
-      !                             !==  e3t  ,  e3u  ,  e3v  ==!
-      !      
+      !                             !==  e3t  ,  e3u  ,  e3v ,  e3w  ==!
       DO jk = 1, jpkm1
           e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
@@ -785,14 +758,7 @@
           e3v_b(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
       END DO
-      !      
-      !                             !==  e3w  ,  e3uw  ,  e3vw  ==!
-      !
-       e3w_b(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
-      e3uw_b(:,:,1) = e3uw_0(:,:,1) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,1) )
-      e3vw_b(:,:,1) = e3vw_0(:,:,1) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,1) )
+      e3w_b(:,:,1) = e3w_0(:,:,1) * (  1._wp + zssht_h(:,:) *  tmask(:,:,1)  )
       DO jk = 2, jpk
-          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX(  tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
-         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * MAX( wumask(:,:,jk-1) , wumask(:,:,jk) ) )
-         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * MAX( wvmask(:,:,jk-1) , wvmask(:,:,jk) ) )
+         e3w_b(:,:,jk) = e3w_0(:,:,jk) * (  1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) , tmask(:,:,jk) )  )
       END DO
       !   
Index: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90
===================================================================
--- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90	(revision 10023)
+++ NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90	(revision 10030)
@@ -102,6 +102,4 @@
       e3v_b (:,:,:) = e3v_n (:,:,:)
       !
-      e3uw_b (:,:,:) = e3uw_n (:,:,:)
-      e3vw_b (:,:,:) = e3vw_n (:,:,:)
       gdept_b(:,:,:) = gdept_n(:,:,:)
       gdepw_b(:,:,:) = gdepw_n(:,:,:)
@@ -206,6 +204,5 @@
          CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
          !                                !      e3t, e3u , e3v, e3f         (from 1 to jpkm1)
-         !                                !      e3w, e3uw, e3vw             (from 1 to jpk  )
-         !                                !      gdept, gdepw, gde3w         (from 1 to jpk  )
+         !                                !      e3w, gdept, gdepw, gde3w    (from 1 to jpk  )
          !
 
Index: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynzdf.F90	(revision 10023)
+++ NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynzdf.F90	(revision 10030)
@@ -36,6 +36,4 @@
 
    PUBLIC   dyn_zdf   !  routine called by step.F90
-
-   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise 
 
    !! * Substitutions
@@ -70,10 +68,13 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER  ::   ji, jj, jk            ! dummy loop indices
-      INTEGER  ::   iku, ikv              ! local integers
-      REAL(wp) ::   zzwi, ze3ua, z2dt_2   ! local scalars
-      REAL(wp) ::   zzws, ze3va           !   -      -
+      INTEGER  ::   ji, jj, jk    ! dummy loop indices
+      INTEGER  ::   iku, ikv      ! local integers
+      REAL(wp) ::   zzwi, zDt_2   ! local scalars
+      REAL(wp) ::   zzws          !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwd, zws   ! 3D workspace 
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv    !  -      -
+      REAL(wp)                     ::   ze3uw_A , ze3uw_Ap1   ! local real
+      REAL(wp)                     ::   ze3vw_A , ze3vw_Ap1   ! local real
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsshu_hA, zsshv_hA    ! 2D workspace
       !!---------------------------------------------------------------------
       !
@@ -84,11 +85,7 @@
          IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
-         !
-         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
-         ELSE                   ;    r_vvl = 1._wp
-         ENDIF
-      ENDIF
-      !
-      z2dt_2 = rDt * 0.5_wp        !* =rn_Dt except in 1st Euler time step where it is equal to rn_Dt/2
+      ENDIF
+      !
+      zDt_2 = rDt * 0.5_wp          !* =rn_Dt except in 1st Euler time step where it is equal to rn_Dt/2
       !
       !
@@ -131,8 +128,6 @@
                iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
                ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
-               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
-               ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
-               va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
+               ua(ji,jj,iku) = ua(ji,jj,iku) + zDt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / e3u_a(ji,jj,iku)
+               va(ji,jj,ikv) = va(ji,jj,ikv) + zDt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / e3v_a(ji,jj,ikv)
             END DO
          END DO
@@ -142,8 +137,6 @@
                   iku = miku(ji,jj)         ! top ocean level at u- and v-points 
                   ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
-                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
-                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
-                  ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
-                  va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
+                  ua(ji,jj,iku) = ua(ji,jj,iku) + zDt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / e3u_a(ji,jj,iku)
+                  va(ji,jj,ikv) = va(ji,jj,ikv) + zDt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / e3v_a(ji,jj,ikv)
                END DO
             END DO
@@ -153,4 +146,26 @@
       !              !==  Vertical diffusion on u  ==!
       !
+      !                          !*  multiplicative factors on e3uw(Naa) and e3vw(Naa)
+      !
+      IF( ln_linssh ) THEN             !--  linear ssh case 
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               zsshu_hA(ji,jj) = 0._wp       ! no time variation in e3
+               zsshv_hA(ji,jj) = 0._wp
+            END DO
+         END DO
+      ELSE                             !--  Non linear ssh case
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               zsshu_hA(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji+1,jj  ,Naa) ) * r1_hu_0(ji,jj) * ssumask(ji,jj)
+               zsshv_hA(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji  ,jj+1,Naa) ) * r1_hv_0(ji,jj) * ssvmask(ji,jj)
+            END DO
+         END DO
+      ENDIF
+
+
+
+
+
       SELECT CASE( nldf_dyn )    !* Matrix construction
       !
@@ -159,9 +174,13 @@
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - rDt * ( 0.5 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) + akzu(ji,jj,jk  ) )   &
-                     &            / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                  zzws = - rDt * ( 0.5 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) + akzu(ji,jj,jk+1) )   &
-                     &            / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+!!gm Note that below, since ze3uw_A is used in a expression masked by wumask, 
+!!      one can remove wumask from its expression   (same for ze3uw_Ap1
+                  ze3uw_A   = e3uw_0(ji,jj,jk  ) * ( 1._wp + zsshu_hA(ji,jj) * wumask(ji,jj,jk  ) )
+                  ze3uw_Ap1 = e3uw_0(ji,jj,jk+1) * ( 1._wp + zsshu_hA(ji,jj) * wumask(ji,jj,jk+1) )
+                  !
+                  zzwi = - rDt * ( 0.5 * (   avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) +   akzu(ji,jj,jk  ) )   &
+                     &                 / ( e3u_a(ji  ,jj,jk  ) * ze3uw_A         ) * wumask(ji,jj,jk  )
+                  zzws = - rDt * ( 0.5 * (   avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) +   akzu(ji,jj,jk+1) )   &
+                     &                 / ( e3u_a(ji  ,jj,jk  ) * ze3uw_Ap1       ) * wumask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi
                   zws(ji,jj,jk) = zzws
@@ -174,7 +193,11 @@
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - z2dt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                  zzws = - z2dt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+!!gm Note that below, since ze3uw_A is used in a expression masked by wumask, 
+!!      one can remove wumask from its expression   (same for ze3uw_Ap1
+                  ze3uw_A   = e3uw_0(ji,jj,jk  ) * ( 1._wp + zsshu_hA(ji,jj) * wumask(ji,jj,jk  ) )
+                  ze3uw_Ap1 = e3uw_0(ji,jj,jk+1) * ( 1._wp + zsshu_hA(ji,jj) * wumask(ji,jj,jk+1) )
+                  !
+                  zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( e3u_a(ji,jj,jk) * ze3uw_A   ) * wumask(ji,jj,jk  )
+                  zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( e3u_a(ji,jj,jk) * ze3uw_Ap1 ) * wumask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi
                   zws(ji,jj,jk) = zzws
@@ -202,6 +225,5 @@
             DO ji = 2, jpim1
                iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
-               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
-               zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
+               zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u_a(ji,jj,iku)
             END DO
          END DO
@@ -211,6 +233,5 @@
                   !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
                   iku = miku(ji,jj)       ! ocean top level at u- and v-points 
-                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
-                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
+                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u_a(ji,jj,iku)
                END DO
             END DO
@@ -243,6 +264,5 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
-            ua(ji,jj,1) = ua(ji,jj,1) + z2dt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( ze3ua * rho0 ) * umask(ji,jj,1)
+            ua(ji,jj,1) = ua(ji,jj,1) + zDt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( e3u_a(ji,jj,1) * rho0 ) * umask(ji,jj,1)
          END DO
       END DO
@@ -276,9 +296,11 @@
             DO jj = 2, jpjm1   
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - rDt * ( 0.5 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) + akzv(ji,jj,jk  ) )   &
-                     &            / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                  zzws = - rDt * ( 0.5 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) + akzv(ji,jj,jk+1) )   &
-                     &            / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                  ze3vw_A   = e3vw_0(ji,jj,jk  ) * ( 1._wp + zsshv_hA(ji,jj) * wvmask(ji,jj,jk  ) )
+                  ze3vw_Ap1 = e3vw_0(ji,jj,jk+1) * ( 1._wp + zsshv_hA(ji,jj) * wvmask(ji,jj,jk+1) )
+                  !
+                  zzwi = - rDt * ( 0.5 * (   avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) +   akzv(ji,jj,jk  ) )   &
+                     &                 / ( e3v_a(ji,jj  ,jk  ) * ze3vw_A         ) * wvmask(ji,jj,jk  )
+                  zzws = - rDt * ( 0.5 * (   avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) +   akzv(ji,jj,jk+1) )   &
+                     &                 / ( e3v_a(ji,jj  ,jk  ) * ze3vw_Ap1       ) * wvmask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
                   zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
@@ -291,7 +313,9 @@
             DO jj = 2, jpjm1   
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - z2dt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                  zzws = - z2dt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                  ze3vw_A   = e3vw_0(ji,jj,jk  ) * ( 1._wp + zsshv_hA(ji,jj) * wvmask(ji,jj,jk  ) )
+                  ze3vw_Ap1 = e3vw_0(ji,jj,jk+1) * ( 1._wp + zsshv_hA(ji,jj) * wvmask(ji,jj,jk+1) )
+                  !
+                  zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( e3v_a(ji,jj,jk) * ze3vw_A   ) * wvmask(ji,jj,jk  )
+                  zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( e3v_a(ji,jj,jk) * ze3vw_Ap1 ) * wvmask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
                   zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
@@ -318,6 +342,5 @@
             DO ji = 2, jpim1
                ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
-               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
-               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
+               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_a(ji,jj,ikv)           
             END DO
          END DO
@@ -326,6 +349,5 @@
                DO ji = 2, jpim1
                   ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
-                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
-                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
+                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3v_a(ji,jj,ikv)
                END DO
             END DO
@@ -358,6 +380,5 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
          DO ji = fs_2, fs_jpim1   ! vector opt.          
-            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
-            va(ji,jj,1) = va(ji,jj,1) + z2dt_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( ze3va * rho0 ) * vmask(ji,jj,1)
+            va(ji,jj,1) = va(ji,jj,1) + zDt_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( e3v_a(ji,jj,1) * rho0 ) * vmask(ji,jj,1)
          END DO
       END DO
@@ -402,107 +423,4 @@
    END SUBROUTINE dyn_zdf
 
-!!gm currently not used : just for memory to be able to add dissipation trend through vertical mixing
-
-   SUBROUTINE zdf_trd( ptrdu, ptrdv ,kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE zdf_trd  ***
-      !!
-      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
-      !!              together with the Leap-Frog time stepping using an 
-      !!              implicit scheme.
-      !!
-      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
-      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
-      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
-      !!               - update the after velocity with the implicit vertical mixing.
-      !!      This requires to solver the following system: 
-      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ]
-      !!      with the following surface/top/bottom boundary condition:
-      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
-      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
-      !!
-      !! ** Action :   (ua,va)   after velocity 
-      !!---------------------------------------------------------------------
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptrdu, ptrdv   ! 3D work arrays use for zdf trends diag
-      INTEGER , INTENT(in   )                         ::   kt             ! ocean time-step index
-      !
-      INTEGER  ::   ji, jj, jk       ! dummy loop indices
-      REAL(wp) ::   zzz              ! local scalar
-      REAL(wp) ::   zavmu, zavmum1   !   -      -
-      REAL(wp) ::   zavmv, zavmvm1   !   -      -
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   z2d    !  -      -
-      !!---------------------------------------------------------------------
-      !
-      CALL lbc_lnk_multi( ua, 'U', -1., va, 'V', -1. )   ! apply lateral boundary condition on (ua,va)
-      !
-      !
-      !                 !==  momentum trend due to vertical diffusion  ==!
-      !
-      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
-         ptrdu(:,:,:) = (              ua(:,:,:) -              ub(:,:,:) )                * r1_Dt - ptrdu(:,:,:)
-         ptrdv(:,:,:) = (              va(:,:,:) -              vb(:,:,:) )                * r1_Dt - ptrdv(:,:,:)
-      ELSE                                      ! applied on thickness weighted velocity
-         ptrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt - ptrdu(:,:,:)
-         ptrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt - ptrdv(:,:,:)
-      ENDIF
-      CALL trd_dyn( ptrdu, ptrdv, jpdyn_zdf, kt )
-      !
-      !
-      !                 !==  KE dissipation trend due to vertical diffusion  ==!
-      !
-      IF( iom_use( 'dispkevfo' ) ) THEN   ! ocean kinetic energy dissipation per unit area
-         !                                ! due to v friction (v=vertical) 
-         !                                ! see NEMO_book appendix C, §C.8 (N.B. here averaged at t-points)
-         !                                ! Note that formally, in a Leap-Frog environment, the shear**2 should be the product of 
-         !                                ! now by before shears, i.e. the source term of TKE (local positivity is not ensured).
-         !                                ! Note also that now e3 scale factors are used as after one are not computed !
-         !
-         ALLOCATE( z2d(jpi,jpj) )
-         z2d(:,:) = 0._wp
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  zavmu   = 0.5 * ( avm(ji+1,jj,jk) + avm(ji  ,jj,jk) )
-                  zavmum1 = 0.5 * ( avm(ji  ,jj,jk) + avm(ji-1,jj,jk) )
-                  zavmv   = 0.5 * ( avm(ji,jj+1,jk) + avm(ji,jj  ,jk) )
-                  zavmvm1 = 0.5 * ( avm(ji,jj  ,jk) + avm(ji,jj-1,jk) )
-               
-                  z2d(ji,jj) = z2d(ji,jj)  +  (                                                                                  &
-                     &   zavmu   * ( ua(ji  ,jj,jk-1) - ua(ji  ,jj,jk) )**2 / e3uw_n(ji  ,jj,jk) * wumask(ji  ,jj,jk)   &
-                     & + zavmum1 * ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) )**2 / e3uw_n(ji-1,jj,jk) * wumask(ji-1,jj,jk)   &
-                     & + zavmv   * ( va(ji,jj  ,jk-1) - va(ji,jj  ,jk) )**2 / e3vw_n(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   &
-                     & + zavmvm1 * ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) )**2 / e3vw_n(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   &
-                     &                        )
-!!gm --- This trends is in fact properly computed in zdf_sh2 but with a backward shift of one time-step  ===>>> use it ?
-!!                                                                                     No since in zdfshé only kz tke (or gls) is used
-!!
-!!gm --- formally, as done below, in a Leap-Frog environment, the shear**2 should be the product of
-!!gm     now by before shears, i.e. the source term of TKE (local positivity is not ensured).
-!!       CAUTION: requires to compute e3uw_a and e3vw_a !!!
-!                  z2d(ji,jj) = z2d(ji,jj)  + (                                                                                 &
-!                     &    avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) / e3uw_n(ji  ,jj,jk)                        &
-!                     &                     * ( ua(ji  ,jj,jk-1) - ua(ji  ,jj,jk) ) / e3uw_a(ji  ,jj,jk) * wumask(ji  ,jj,jk)   &
-!                     &  + avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) / e3uw_n(ji-1,jj,jk)                        &
-!                     &                       ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) ) / e3uw_a(ji-1,jj,jk) * wumask(ji-1,jj,jk)   &
-!                     &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) / e3vw_n(ji,jj  ,jk)                        &
-!                     &                       ( va(ji,jj  ,jk-1) - va(ji,jj  ,jk) ) / e3vw_a(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   &
-!                     &  + avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) / e3vw_n(ji,jj-1,jk)                        &
-!                     &                       ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) ) / e3vw_a(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   &
-!                     &                       )
-!!gm end
-               END DO
-            END DO
-         END DO
-         zzz= - 0.5_wp* rho0           ! caution sign minus here
-         z2d(:,:) = zzz * z2d(:,:) 
-         CALL lbc_lnk( z2d,'T', 1. )
-         CALL iom_put( 'dispkevfo', z2d )
-         DEALLOCATE( z2d )
-      ENDIF
-      !
-   END SUBROUTINE zdf_trd
-   
-!!gm end not used
-
    !!==============================================================================
 END MODULE dynzdf
Index: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trazdf.F90
===================================================================
--- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trazdf.F90	(revision 10023)
+++ NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trazdf.F90	(revision 10030)
@@ -112,5 +112,5 @@
       !!
       !! ** Method  :  The vertical diffusion of a tracer ,t , is given by:
-      !!          difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
+      !!          difft = dz( avt dz(t) ) = 1/e3t(Naa) dk+1( avt/e3w(Naa) dk(t) )
       !!      It is computed using a backward time scheme (t=after field)
       !!      which provide directly the after tracer field.
@@ -175,7 +175,6 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-!!gm BUG   here e3w_a  should be used !!!!!   but then should be added in the system 
-                     zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )
-                     zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
+                     zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_a(ji,jj,jk  )
+                     zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_a(ji,jj,jk+1)
                      zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)
                  END DO
Index: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfosm.F90	(revision 10023)
+++ NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfosm.F90	(revision 10030)
@@ -247,4 +247,8 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity
+      
+      REAL(wp)                     ::   ze3uw_BN, ze3vw_BN   ! use for e3uw, e3vw computation
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsshu_hB, zsshv_hB   ! at Before and Now time-step
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsshu_hN, zsshv_hN   
 
       ! For debugging
@@ -1235,18 +1239,45 @@
 
 ! KPP-style Ri# mixing
-       IF( ln_kpprimix) THEN
-          DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
-             DO jj = 1, jpjm1
-                DO ji = 1, jpim1   ! vector opt.
-                   z3du(ji,jj,jk) = 0.5 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   &
-                        &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) &
-                        &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
-                   z3dv(ji,jj,jk) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   &
-                        &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) &
-                        &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
+      IF( ln_kpprimix) THEN
+         ! 
+         IF( ln_linssh ) THEN      !==      linear ssh case  ==!
+            DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1   ! vector opt.
+                     ze3uw_BN = e3uw_0(ji,jj,jk) * e3uw_0(ji,jj,jk)
+                     ze3vw_BN = e3vw_0(ji,jj,jk) * e3vw_0(ji,jj,jk)
+                     z3du(ji,jj,jk) = 0.5 * (  un(ji,jj,jk-1) -  un(ji,jj,jk) )            &
+                        &                 * (  ub(ji,jj,jk-1) -  ub(ji,jj,jk) ) / ze3uw_BN * wumask(ji,jj,jk)
+                     z3dv(ji,jj,jk) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj,jk) )            &
+                        &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj,jk) ) / ze3vw_BN * wvmask(ji,jj,jk)
+                   END DO
                 END DO
-             END DO
-          END DO
-      !
+            END DO
+         ELSE                       !==  Non linear ssh case  ==!
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  zsshu_hB(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * r1_hu_0(ji,jj) * ssumask(ji,jj)
+                  zsshv_hB(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * r1_hv_0(ji,jj) * ssvmask(ji,jj)
+                  zsshu_hN(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * r1_hu_0(ji,jj) * ssumask(ji,jj)
+                  zsshv_hN(ji,jj) = 0.5_wp * ( ssh(ji,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * r1_hv_0(ji,jj) * ssvmask(ji,jj)
+               END DO
+            END DO
+            !
+            DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1   ! vector opt.
+                     ze3uw_BN = e3uw_0(ji,jj,jk) * e3uw_0(ji,jj,jk) * ( 1._wp + zsshu_hB(ji,jj) * wumask(ji,jj,jk) )   &
+                        &                                           * ( 1._wp + zsshu_hN(ji,jj) * wumask(ji,jj,jk) )
+                     ze3vw_BN = e3vw_0(ji,jj,jk) * e3vw_0(ji,jj,jk) * ( 1._wp + zsshu_hB(ji,jj) * wvmask(ji,jj,jk) )   &
+                        &                                           * ( 1._wp + zsshu_hN(ji,jj) * wvmask(ji,jj,jk) )
+                     z3du(ji,jj,jk) = 0.5 * (  un(ji,jj,jk-1) -  un(ji,jj,jk) )            &
+                        &                 * (  ub(ji,jj,jk-1) -  ub(ji,jj,jk) ) / ze3uw_BN * wumask(ji,jj,jk)
+                     z3dv(ji,jj,jk) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj,jk) )            &
+                        &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj,jk) ) / ze3vw_BN * wvmask(ji,jj,jk)
+                   END DO
+                END DO
+            END DO
+         ENDIF
+         !
          DO jk = 2, jpkm1
             DO jj = 2, jpjm1
@@ -1262,15 +1293,15 @@
                 END DO
              END DO
-          END DO
-
-          DO jj = 2, jpjm1
-             DO ji = 2, jpim1
-                DO jk = ibld(ji,jj) + 1, jpkm1
-                   zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
-                   zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
-                END DO
-             END DO
-          END DO
-
+         END DO
+         !
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               DO jk = ibld(ji,jj) + 1, jpkm1
+                  zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
+                  zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
+               END DO
+            END DO
+         END DO
+         !
        END IF ! ln_kpprimix = .true.
 
@@ -1708,10 +1739,6 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               ua(ji,jj,jk) =  ua(ji,jj,jk)                      &
-                  &                 - (  ghamu(ji,jj,jk  )  &
-                  &                    - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
-               va(ji,jj,jk) =  va(ji,jj,jk)                      &
-                  &                 - (  ghamv(ji,jj,jk  )  &
-                  &                    - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
+               ua(ji,jj,jk) = ua(ji,jj,jk) - (  ghamu(ji,jj,jk) - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) - (  ghamv(ji,jj,jk) - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
             END DO
          END DO
Index: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfphy.F90
===================================================================
--- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfphy.F90	(revision 10023)
+++ NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfphy.F90	(revision 10030)
@@ -246,6 +246,6 @@
       !
       IF( l_zdfsh2 )   &         !* shear production at w-points (energy conserving form)
-         CALL zdf_sh2( ub, vb, un, vn, avm_k,   &     ! <<== in
-            &                           zsh2    )     ! ==>> out : shear production
+         CALL zdf_sh2( ssh, ub, vb, un, vn, avm_k,   &     ! <<== in
+            &                                zsh2    )     ! ==>> out : shear production
       !
       SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points
Index: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfsh2.F90
===================================================================
--- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfsh2.F90	(revision 10023)
+++ NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ZDF/zdfsh2.F90	(revision 10030)
@@ -6,4 +6,5 @@
    !! History :   -   !  2014-10  (A. Barthelemy, G. Madec)  original code
    !!   NEMO     4.0  !  2017-04  (G. Madec)  remove u-,v-pts avm
+   !!            5.0  !  2018-08  (G. Madec)  local calculation of e3uw e3vw
    !!----------------------------------------------------------------------
 
@@ -28,5 +29,5 @@
 CONTAINS
 
-   SUBROUTINE zdf_sh2( pub, pvb, pun, pvn, p_avm, p_sh2  ) 
+   SUBROUTINE zdf_sh2( pssh, pub, pvb, pun, pvn, p_avm, p_sh2  ) 
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zdf_sh2  ***
@@ -47,30 +48,72 @@
       !! References :   Bruchard, OM 2002
       !! ---------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pssh                 ! sea surface height
       REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pub, pvb, pun, pvn   ! before, now horizontal velocities
       REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm                ! vertical eddy viscosity (w-points)
       REAL(wp), DIMENSION(:,:,:) , INTENT(  out) ::   p_sh2                ! shear production of TKE (w-points)
       !
-      INTEGER  ::   ji, jj, jk   ! dummy loop arguments
-      REAL(wp), DIMENSION(jpi,jpj) ::   zsh2u, zsh2v   ! 2D workspace
+      INTEGER ::   ji, jj, jk           ! dummy loop arguments
+      REAL(wp)::   ze3uw_BN, ze3vw_BN   ! local real
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsh2u   , zsh2v     ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsshu_hB, zsshv_hB  ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsshu_hN, zsshv_hN
       !!--------------------------------------------------------------------
       !
-      DO jk = 2, jpkm1
-         DO jj = 1, jpjm1        !* 2 x shear production at uw- and vw-points (energy conserving form)
+      IF( ln_linssh ) THEN       !==  linear ssh case  ==!
+         DO jk = 2, jpkm1
+            DO jj = 1, jpjm1        !* 2 x shear production at uw- and vw-points (energy conserving form)
+               DO ji = 1, jpim1
+                  ze3uw_BN = e3uw_0(ji,jj,jk) * e3uw_0(ji,jj,jk)
+                  ze3vw_BN = e3vw_0(ji,jj,jk) * e3vw_0(ji,jj,jk)
+                  !
+                  zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) )            &
+                     &         * (   pun(ji,jj,jk-1) -   pun(ji,jj,jk) )            &
+                     &         * (   pub(ji,jj,jk-1) -   pub(ji,jj,jk) ) / ze3uw_BN * wumask(ji,jj,jk)
+                  zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) )            &
+                     &         * (   pvn(ji,jj,jk-1) -   pvn(ji,jj,jk) )            &
+                     &         * (   pvb(ji,jj,jk-1) -   pvb(ji,jj,jk) ) / ze3vw_BN * wvmask(ji,jj,jk)
+               END DO
+            END DO
+            DO jj = 2, jpjm1        !* shear production at w-point
+               DO ji = 2, jpim1           ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)
+                  p_sh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   &
+                     &                       + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   )
+               END DO
+            END DO
+         END DO 
+         !
+      ELSE                       !==  Non linear ssh case  ==!
+         DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) &
-                  &         * (   pun(ji,jj,jk-1) -   pun(ji,jj,jk) ) &
-                  &         * (   pub(ji,jj,jk-1) -   pub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) * wumask(ji,jj,jk)
-               zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) &
-                  &         * (   pvn(ji,jj,jk-1) -   pvn(ji,jj,jk) ) &
-                  &         * (   pvb(ji,jj,jk-1) -   pvb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) * wvmask(ji,jj,jk)
+               zsshu_hB(ji,jj) = 0.5_wp * ( pssh(ji,jj,Nbb) + pssh(ji+1,jj  ,Nbb) ) * r1_hu_0(ji,jj) * ssumask(ji,jj)
+               zsshv_hB(ji,jj) = 0.5_wp * ( pssh(ji,jj,Nbb) + pssh(ji  ,jj+1,Nbb) ) * r1_hv_0(ji,jj) * ssvmask(ji,jj)
+               zsshu_hN(ji,jj) = 0.5_wp * ( pssh(ji,jj,Nnn) + pssh(ji+1,jj  ,Nnn) ) * r1_hu_0(ji,jj) * ssumask(ji,jj)
+               zsshv_hN(ji,jj) = 0.5_wp * ( pssh(ji,jj,Nnn) + pssh(ji  ,jj+1,Nnn) ) * r1_hv_0(ji,jj) * ssvmask(ji,jj)
             END DO
          END DO
-         DO jj = 2, jpjm1        !* shear production at w-point
-            DO ji = 2, jpim1           ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)
-               p_sh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   &
-                  &                       + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   )
+         !
+         DO jk = 2, jpkm1
+            DO jj = 1, jpjm1        !* 2 x shear production at uw- and vw-points (energy conserving form)
+               DO ji = 1, jpim1
+                  ze3uw_BN = e3uw_0(ji,jj,jk) * e3uw_0(ji,jj,jk) * ( 1._wp + zsshu_hB(ji,jj) * wumask(ji,jj,jk) )   &
+                     &                                           * ( 1._wp + zsshu_hN(ji,jj) * wumask(ji,jj,jk) )
+                  ze3vw_BN = e3vw_0(ji,jj,jk) * e3vw_0(ji,jj,jk) * ( 1._wp + zsshu_hB(ji,jj) * wvmask(ji,jj,jk) )   &
+                     &                                           * ( 1._wp + zsshu_hN(ji,jj) * wvmask(ji,jj,jk) )
+                  zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) &
+                     &         * (   pun(ji,jj,jk-1) -   pun(ji,jj,jk) ) &
+                     &         * (   pub(ji,jj,jk-1) -   pub(ji,jj,jk) ) / ze3uw_BN * wumask(ji,jj,jk)
+                  zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) &
+                     &         * (   pvn(ji,jj,jk-1) -   pvn(ji,jj,jk) ) &
+                     &         * (   pvb(ji,jj,jk-1) -   pvb(ji,jj,jk) ) / ze3vw_BN * wvmask(ji,jj,jk)
+               END DO
             END DO
-         END DO
-      END DO 
+            DO jj = 2, jpjm1        !* shear production at w-point
+               DO ji = 2, jpim1           ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)
+                  p_sh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   &
+                     &                       + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   )
+               END DO
+            END DO
+         END DO 
+      ENDIF
       !
    END SUBROUTINE zdf_sh2
