Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domzgr_substitute.h90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 13694)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domzgr_substitute.h90	(revision 13696)
@@ -16,4 +16,5 @@
 #   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1._wp+r3v(i,j,t)*vmask(i,j,k)))
 #   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1._wp+r3f(i,j)*fmask(i,j,k)))
+#   define  e3f_vor(i,j,k) (e3f_0vor(i,j,k)*(1._wp+r3f(i,j)*fmask(i,j,k)))
 #   define  e3w(i,j,k,t)   (e3w_0(i,j,k)*(1._wp+r3t(i,j,t)))
 #   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1._wp+r3u(i,j,t)))
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynspg_ts.F90	(revision 13694)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynspg_ts.F90	(revision 13696)
@@ -1119,5 +1119,5 @@
       !! 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
+      !! and update depths at T- points (ht) at each barotropic time step
       !!
       !! Compute zwz = f / ( height of the water colomn )
@@ -1126,12 +1126,11 @@
       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( np_EEN, np_ENE, np_ENS , np_MIX )   !=  schemes using the same e3f definition
+         SELECT CASE( nn_e3f_typ )                  !* ff_f/e3 at F-point
          CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 0, 0, 0, 0 )
                zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1)   &
                     &       + ht(ji,jj  ) + ht(ji+1,jj  ) ) * 0.25_wp  
@@ -1139,5 +1138,5 @@
             END_2D
          CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 0, 0, 0, 0 )
                zwz(ji,jj) =     (    ht(ji,jj+1) +     ht(ji+1,jj+1)      &
                     &            +   ht(ji,jj  ) +     ht(ji+1,jj  )  )   &
@@ -1148,4 +1147,8 @@
          END SELECT
          CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
+      END SELECT
+      !
+      SELECT CASE( nvor_scheme )
+      CASE( np_EEN )
          !
          ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp
@@ -1157,5 +1160,5 @@
          END_2D
          !
-      CASE( np_EET )                  != EEN scheme using e3t energy conserving scheme
+      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_2D( 0, 1, 0, 1 )
@@ -1167,56 +1170,7 @@
          END_2D
          !
-      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_2D( 1, 0, 1, 0 )
-               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_2D
-         ENDIF
-         !
-         DO jj = 1, jpjm1   ! keep only the value at the coast if sco
-            zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
-         END DO
-         !
-         DO jk = 1, jpkm1   ! ocean point : sum of masked e3f
-            DO jj = 1, jpjm1
-               zhf(:,jj) = zhf(:,jj) + e3f(:,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_2D( 1, 1, 1, 1 )
-            IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj)
-         END_2D
-         zwz(:,:) = ff_f(:,:) * zwz(:,:)
       END SELECT
       
    END SUBROUTINE dyn_cor_2d_init
-
 
 
@@ -1412,5 +1366,4 @@
    END SUBROUTINE wad_spg
      
-
 
    SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
Index: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynvor.F90
===================================================================
--- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynvor.F90	(revision 13694)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynvor.F90	(revision 13696)
@@ -60,7 +60,7 @@
    LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET)
    LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN)
-   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
    LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX)
    LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
+   INTEGER, PUBLIC ::   nn_e3f_typ      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
 
    INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme
@@ -84,5 +84,7 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       - 
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       - 
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -
+   !
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   e3f_0vor   ! e3f used in EEN, ENE and ENS cases (key_qco only)
    
    REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
@@ -335,5 +337,5 @@
       !
       INTEGER  ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
+      REAL(wp) ::   zx1, zy1, zx2, zy2, ze3f, zmsk   ! local scalars
       REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace
       !!----------------------------------------------------------------------
@@ -387,19 +389,36 @@
          !
 #if defined key_qco
-         !                                   !==  horizontal fluxes and potential vorticity  ==!
+         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco)
+            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk)
+         END_2D
+#else
+         SELECT CASE( nn_e3f_typ  )           !==  potential vorticity  ==!
+         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
+            DO_2D( 1, 0, 1, 0 )
+               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
+                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
+                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
+                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
+               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
+               ELSE                       ;   zwz(ji,jj) = 0._wp
+               ENDIF
+            END_2D
+         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
+            DO_2D( 1, 0, 1, 0 )
+               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
+                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
+                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
+                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   )
+               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   &
+                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   )
+               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
+               ELSE                       ;   zwz(ji,jj) = 0._wp
+               ENDIF
+            END_2D
+         END SELECT
+#endif
+         !                                   !==  horizontal fluxes  ==!
          zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
          zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
-         zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
-#else
-         !                                   !==  horizontal fluxes and potential vorticity  ==!
-         IF( ln_sco ) THEN                        ! weighted by e3
-            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
-         ELSE                                     ! ugly fix: unweighted by e3
-            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
-         ENDIF
-#endif
          !
          !                                   !==  compute and add the vorticity term trend  =!
@@ -445,5 +464,5 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zuav, zvau   ! local scalars
+      REAL(wp) ::   zuav, zvau, ze3f, zmsk   ! local scalars
       REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace
       !!----------------------------------------------------------------------
@@ -497,19 +516,36 @@
          !
 #if defined key_qco
-         !                                   !==  horizontal fluxes and potential vorticity  ==!
+         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco)
+            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk)
+         END_2D
+#else
+         SELECT CASE( nn_e3f_typ )           !==  potential vorticity  ==!
+         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
+            DO_2D( 1, 0, 1, 0 )
+               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
+                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
+                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
+                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
+               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f
+               ELSE                       ;   zwz(ji,jj) = 0._wp
+               ENDIF
+            END_2D
+         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
+            DO_2D( 1, 0, 1, 0 )
+               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   &
+                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
+                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   &
+                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   )
+               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   &
+                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   )
+               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f
+               ELSE                       ;   zwz(ji,jj) = 0._wp
+               ENDIF
+            END_2D
+         END SELECT
+#endif
+         !                                   !==  horizontal fluxes  ==!
          zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
          zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
-         zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
-#else
-         !                                   !==  horizontal fluxes and potential vorticity  ==!
-         IF( ln_sco ) THEN                        ! weighted by e3
-            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
-         ELSE                                     ! ugly fix: unweighted by e3
-            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
-         ENDIF
-#endif
          !
          !                                   !==  compute and add the vorticity term trend  =!
@@ -570,5 +606,10 @@
          !                                             ! ===============
          !
-         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
+#if defined key_qco
+         DO_2D( 1, 0, 1, 0 )                 ! == reciprocal of e3 at F-point (key_qco)
+            z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk)
+         END_2D
+#else
+         SELECT CASE( nn_e3f_typ )           ! == reciprocal of e3 at F-point
          CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
             DO_2D( 1, 0, 1, 0 )
@@ -594,4 +635,5 @@
             END_2D
          END SELECT
+#endif
          !
          SELECT CASE( kvor )                 !==  vorticity considered  ==!
@@ -801,7 +843,8 @@
       INTEGER ::   ji, jj, jk    ! dummy loop indices
       INTEGER ::   ioptio, ios   ! local integer
+      REAL(wp) ::   zmsk    ! local scalars
       !!
       NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
-         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk
+         &                 ln_dynvor_een, nn_e3f_typ   , ln_dynvor_mix, ln_dynvor_msk
       !!----------------------------------------------------------------------
       !
@@ -825,5 +868,5 @@
          WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
          WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
-         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
+         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_e3f_typ = ', nn_e3f_typ
          WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
          WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
@@ -891,5 +934,41 @@
          !
       END SELECT
-      
+!!st ADD !! pense a recalculer le hf_0
+#if defined key_qco
+      SELECT CASE( nvor_scheme )    ! qco case: pre-computed a specific e3f_0 for some vorticity schemes
+      CASE( np_ENS , np_ENE , np_EEN , np_MIX )
+         !
+         ALLOCATE( e3f_0vor(jpi,jpj,jpk) )
+         !
+         SELECT CASE( nn_e3f_typ )
+         CASE ( 0 )                        ! original formulation  (masked averaging of e3t divided by 4)
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   &
+                  &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
+                  &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   &
+                  &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) * 0.25_wp
+            END_3D
+         CASE ( 1 )                        ! new formulation  (masked averaging of e3t divided by the sum of mask)
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               zmsk = (tmask(ji,jj+1,jk) +tmask(ji+1,jj+1,jk)   &
+                  &  + tmask(ji,jj  ,jk) +tmask(ji+1,jj  ,jk)  )
+               !
+               IF( zmsk /= 0._wp ) THEN 
+                  e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   &
+                     &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
+                     &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   &
+                     &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) / zmsk
+               ENDIF
+            END_3D
+         END SELECT
+         !
+         CALL lbc_lnk( 'dynvor', e3f_0vor, 'F', 1._wp )
+         !                                 ! insure e3f_0vor /= 0
+         WHERE( e3f_0vor(:,:,:) == 0._wp )   e3f_0vor(:,:,:) = e3f_0(:,:,:)
+         !
+      END SELECT
+      !
+#endif
+!!st END      
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
