Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12491)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90	(revision 12492)
@@ -23,4 +23,6 @@
    USE phycst          ! physical constant
    USE dom_oce         ! ocean space and time domain
+   USE dynadv  , ONLY : ln_dynadv_vec
+   USE isf_oce         ! iceshelf cavities
    USE sbc_oce         ! ocean surface boundary condition
    USE wet_dry         ! wetting and drying
@@ -109,5 +111,6 @@
       !
    END SUBROUTINE dom_qe_init
-   !
+
+
    SUBROUTINE dom_qe_zgr(Kbb, Kmm, Kaa)
       !!----------------------------------------------------------------------
@@ -139,17 +142,27 @@
       !                    !== Set of all other vertical scale factors  ==!  (now and before)
       !                                ! Horizontal interpolation of e3t
-      CALL dom_qe_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U
-      CALL dom_qe_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
-      CALL dom_qe_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V
-      CALL dom_qe_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
-      CALL dom_qe_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F
-      !                                ! Vertical interpolation of e3t,u,v
-      CALL dom_qe_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W
-      CALL dom_qe_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  )
-      CALL dom_qe_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW
-      CALL dom_qe_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
-      CALL dom_qe_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW
-      CALL dom_qe_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
-
+      CALL dom_qe_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) )
+      CALL dom_qe_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
+      !
+      DO jk = 1, jpkm1                    ! Horizontal interpolation of e3t
+         e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) * tmask(:,:,jk) )   ! Kbb time level
+         e3u(:,:,jk,Kbb) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) * umask(:,:,jk) )
+         e3v(:,:,jk,Kbb) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) * vmask(:,:,jk) )
+         e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) * tmask(:,:,jk) )   ! Kmm time level
+         e3u(:,:,jk,Kmm) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) * umask(:,:,jk) )
+         e3v(:,:,jk,Kmm) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) * vmask(:,:,jk) )
+         e3f(:,:,jk)     = e3f_0(:,:,jk) * ( 1._wp + r3f(:,:)     * fmask(:,:,jk) )
+      END DO
+      !
+      DO jk = 1, jpk                      ! Vertical interpolation of e3t,u,v
+         !                                   ! The ratio does not have to be masked at w-level
+         e3w (:,:,jk,Kbb) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )   ! Kbb time level
+         e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) )
+         e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) )
+         e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
+         e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
+         e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
+      END DO
+      !
       ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
       e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
@@ -158,37 +171,43 @@
       !
       !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
-      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)       ! reference to the ocean surface (used for MLD and light penetration)
-      gdepw(:,:,1,Kmm) = 0.0_wp
-      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
-      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb)
-      gdepw(:,:,1,Kbb) = 0.0_wp
-      DO_3D_11_11( 2, jpk )
-         !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
-         !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
-         !                             ! 0.5 where jk = mikt
-!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ??
-         zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
-         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
-         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
-            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))
-         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
-         gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
-         gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  &
-            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))
-      END_3D
+      IF( ln_isf ) THEN          !** IceShelF cavities
+      !                             ! to be created depending of the new names in isf
+      !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
+      !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
+!!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
+         gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
+         gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
+         gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
+         DO jk = 2, jpk
+            gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
+                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
+            gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
+                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
+            gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
+            gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
+                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
+            gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
+                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
+         END DO
+         !
+      ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
+         !
+!!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
+         DO jk = 1, jpk
+            gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
+            gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
+            gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
+            gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
+            gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
+         END DO
+         !
+      ENDIF
       !
       !                    !==  thickness of the water column  !!   (ocean portion only)
-      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
-      hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1)
-      hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1)
-      hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1)
-      hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1)
-      DO jk = 2, jpkm1
-         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
-         hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
-         hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
-         hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
-         hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
-      END DO
+      ht(:,:)     = ht_0(:,:)           + ssh(:,:,Kmm)
+      hu(:,:,Kbb) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kbb) )
+      hu(:,:,Kmm) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kmm) )
+      hv(:,:,Kbb) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kbb) )
+      hv(:,:,Kmm) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kmm) )
       !
       !                    !==  inverse of water column thickness   ==!   (u- and v- points)
@@ -251,29 +270,23 @@
       !                                                ! --------------------------------------------- !
       !
-      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
-      DO jk = 1, jpkm1
-         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
-         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
-      END DO
       !
       ! *********************************** !
       ! After scale factors at u- v- points !
       ! *********************************** !
-
-      CALL dom_qe_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
-      CALL dom_qe_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
-
+      !
+      CALL dom_qe_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa) )
+      !
+      DO jk = 1, jpkm1
+         e3t(:,:,jk,Kaa) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kaa) * tmask(:,:,jk) )
+         e3u(:,:,jk,Kaa) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kaa) * umask(:,:,jk) )
+         e3v(:,:,jk,Kaa) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kaa) * vmask(:,:,jk) )
+      END DO
+      !
       ! *********************************** !
       ! After depths at u- v points         !
       ! *********************************** !
-
-      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
-      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
-      DO jk = 2, jpkm1
-         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
-         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
-      END DO
+      hu(:,:,Kaa) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kaa) )
+      hv(:,:,Kaa) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kaa) )
       !                                        ! Inverse of the local depth
-!!gm BUG ?  don't understand the use of umask_i here .....
       r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
       r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
@@ -329,8 +342,19 @@
       ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
       ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
-
-      CALL dom_qe_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
+      !!st  ! r3t/u/v should be unchanged
+      CALL dom_qe_r3c( ssh(:,:,Kmm), r3t_f(:,:), r3u_f(:,:), r3v_f(:,:), r3f(:,:) )
+      !
+      DO jk = 1, jpkm1                    ! Horizontal interpolation of e3t
+         e3f(:,:,jk)     = e3f_0(:,:,jk) * ( 1._wp + r3f(:,:)     * fmask(:,:,jk) )   ! Kmm time level
+      END DO
+      !CALL dom_qe_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
 
       ! Vertical scale factor interpolations
+      ! DO jk = 1, jpk                      ! Vertical interpolation of e3t,u,v
+      !    !                                   ! The ratio does not have to be masked at w-level
+      !    e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
+      !    e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
+      !    e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
+      ! END DO
       CALL dom_qe_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
       CALL dom_qe_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
@@ -353,12 +377,35 @@
          gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
       END_3D
+!       IF( ln_isf ) THEN          !** IceShelF cavities
+!       !                             ! to be created depending of the new names in isf
+!       !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
+!       !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
+! !!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
+!          gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
+!          gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
+!          gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
+!          DO jk = 2, jpk
+!             gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
+!                               + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
+!             gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
+!                               + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
+!             gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
+!          END DO
+!          !
+!       ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
+!          !
+! !!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
+!          DO jk = 1, jpk
+!             gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
+!             gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
+!             gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
+!          END DO
+!          !
+!       ENDIF
 
       ! Local depth and Inverse of the local depth of the water
       ! -------------------------------------------------------
       !
-      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
-      DO jk = 2, jpkm1
-         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
-      END DO
+      ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm)
 
       ! write restart file
@@ -464,4 +511,69 @@
       !
    END SUBROUTINE dom_qe_interpol
+
+
+   SUBROUTINE dom_qe_r3c( pssh, pr3t, pr3u, pr3v, pr3f )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE r3c  ***
+      !!
+      !! ** Purpose :   compute the filtered ratio ssh/h_0 at t-,u-,v-,f-points
+      !!
+      !! ** Method  : - compute the ssh at u- and v-points (f-point optional)
+      !!                   Vector Form : surface weighted averaging
+      !!                   Flux   Form : simple           averaging
+      !!              - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:)          , INTENT(in   )  ::   pssh               ! sea surface height   [m]
+      REAL(wp), DIMENSION(:,:)          , INTENT(  out)  ::   pr3t, pr3u, pr3v   ! ssh/h0 ratio at t-, u-, v-,points  [-]
+      REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(  out)  ::   pr3f               ! ssh/h0 ratio at f-point   [-]
+      !
+      INTEGER ::   ji, jj   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      !
+      pr3t(:,:) = pssh(:,:) * r1_ht_0(:,:)   !==  ratio at t-point  ==!
+      !
+      !
+      !                                      !==  ratio at u-,v-point  ==!
+      !
+      IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging)
+         DO_2D_11_11
+            pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
+               &                    + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
+            pr3v(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
+               &                    + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
+         END_2D
+      ELSE                                         !- Flux Form   (simple averaging)
+         DO_2D_11_11
+            pr3u(ji,jj) = 0.5_wp * (  pssh(ji  ,jj) + pssh(ji+1,jj)  ) * r1_hu_0(ji,jj)
+            pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj  ) + pssh(ji,jj+1)  ) * r1_hv_0(ji,jj)
+         END_2D
+      ENDIF
+      !
+      IF( .NOT.PRESENT( pr3f ) ) THEN              !- lbc on ratio at u-, v-points only
+         CALL lbc_lnk_multi( 'dom_qe_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp )
+         !
+         !
+      ELSE                                   !==  ratio at f-point  ==!
+         !
+         IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging)
+            DO_2D_01_01                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
+               pr3f(ji,jj) = 0.25_wp * (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  &
+                  &                     + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  &
+                  &                     + e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  &
+                  &                     + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
+            END_2D
+         ELSE                                      !- Flux Form   (simple averaging)
+            DO_2D_01_01                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
+               pr3f(ji,jj) = 0.25_wp * (  pssh(ji  ,jj  ) + pssh(ji+1,jj  )  &
+                  &                     + pssh(ji  ,jj+1) + pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj)
+            END_2D
+         ENDIF
+         !                                                 ! lbc on ratio at u-,v-,f-points
+         CALL lbc_lnk_multi( 'dom_qe_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
+         !
+      ENDIF
+      !
+   END SUBROUTINE dom_qe_r3c
 
 
Index: /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90
===================================================================
--- /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90	(revision 12491)
+++ /NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/steplf.F90	(revision 12492)
@@ -126,4 +126,14 @@
                          CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition (including sea-ice)
 
+      !             !!st      !!!!!!!!!!!!!!!!!!!!!!!
+      !
+      ! emp = 0._wp
+      ! emp_b = 0._wp
+      ! qns = 0._wp
+      ! qsr = 0._wp
+      ! qns_b = 0._wp
+      !
+      !             !!st
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update stochastic parameters and random T/S fluctuations
