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 13620)
+++ NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynvor.F90	(revision 13621)
@@ -21,11 +21,14 @@
    !!             -   ! 2018-03  (G. Madec)  add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
    !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation
+   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  make ln_dynvor_msk truly efficient on relative vorticity
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   dyn_vor       : Update the momentum trend with the vorticity trend
+   !!       vor_enT   : energy conserving scheme at T-pt  (ln_dynvor_enT=T)
+   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
    !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
-   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
    !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
+   !!       vor_eeT   : energy conserving at T-pt         (ln_dynvor_eeT=T)
    !!   dyn_vor_init  : set and control of the different vorticity option
    !!----------------------------------------------------------------------
@@ -80,6 +83,6 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       - 
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2v)/(2*e1e2f)  used in F-point metric term calculation
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1u)/(2*e1e2f)   -        -      -       - 
+   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) ::   r1_4  = 0.250_wp         ! =1/4
@@ -218,5 +221,5 @@
       REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
       REAL(wp), DIMENSION(jpi,jpj)     ::   zwx, zwy, zwt   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz             ! 3D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zwz      ! 3D workspace
       !!----------------------------------------------------------------------
       !
@@ -228,12 +231,14 @@
       !
       !
-      SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
-      CASE ( np_RVO )                           !* relative vorticity
-         DO jk = 1, jpkm1                                 ! Horizontal slab
+      SELECT CASE( kvor )                 !== relative vorticity considered  ==!
+      !
+      CASE ( np_RVO , np_CRV )                  !* relative vorticity at f-point is used
+         ALLOCATE( zwz(jpi,jpj,jpk) )
+         DO jk = 1, jpkm1                                ! Horizontal slab
             DO_2D( 1, 0, 1, 0 )
                zwz(ji,jj,jk) = (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
                   &             - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
-            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
+            IF( ln_dynvor_msk ) THEN                     ! mask relative vorticity 
                DO_2D( 1, 0, 1, 0 )
                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
@@ -241,58 +246,43 @@
             ENDIF
          END DO
-
-         CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
-
-      CASE ( np_CRV )                           !* Coriolis + relative vorticity
-         DO jk = 1, jpkm1                                 ! Horizontal slab
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   &
-                  &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
-            END_2D
-            IF( ln_dynvor_msk ) THEN                     ! mask/unmask relative vorticity 
-               DO_2D( 1, 0, 1, 0 )
-                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
-               END_2D
-            ENDIF
-         END DO
-
-         CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
-
+         CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
+         !
       END SELECT
 
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
-      !                                                ! ===============
-
+         !                                             ! ===============
+         !
          SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
+         !
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
             zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
          CASE ( np_RVO )                           !* relative vorticity
             DO_2D( 0, 1, 0, 1 )
-               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
-                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
-                  &                  * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
+               zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)       &
+                  &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk)   )   &
+                  &              * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
             END_2D
          CASE ( np_MET )                           !* metric term
             DO_2D( 0, 1, 0, 1 )
-               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
-                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) &
-                  &             * e3t(ji,jj,jk,Kmm)
+               zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)       &
+                  &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   )   &
+                  &       * e3t(ji,jj,jk,Kmm)
             END_2D
          CASE ( np_CRV )                           !* Coriolis + relative vorticity
             DO_2D( 0, 1, 0, 1 )
-               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    &
-                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) &
-                  &                                 * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
+               zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)        &
+                  &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  )   &
+                  &       * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
             END_2D
          CASE ( np_CME )                           !* Coriolis + metric
             DO_2D( 0, 1, 0, 1 )
-               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           &
-                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  &
-                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) &
-                    &          * e3t(ji,jj,jk,Kmm)
+               zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                               &
+                    &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)      &
+                    &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  )   &
+                    &     * e3t(ji,jj,jk,Kmm)
             END_2D
          CASE DEFAULT                                             ! error
-            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
+            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor')
          END SELECT
          !
@@ -310,4 +300,9 @@
       END DO                                           !   End of slab
       !                                                ! ===============
+      !
+      SELECT CASE( kvor )        ! deallocate zwz if necessary
+      CASE ( np_RVO , np_CRV )   ;   DEALLOCATE( zwz )
+      END SELECT
+      !
    END SUBROUTINE vor_enT
 
@@ -362,4 +357,9 @@
                   &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
+            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
+               DO_2D( 1, 0, 1, 0 )
+                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
+               END_2D
+            ENDIF
          CASE ( np_MET )                           !* metric term
             DO_2D( 1, 0, 1, 0 )
@@ -372,4 +372,9 @@
                   &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
+            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
+               DO_2D( 1, 0, 1, 0 )
+                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
+               END_2D
+            ENDIF
          CASE ( np_CME )                           !* Coriolis + metric
             DO_2D( 1, 0, 1, 0 )
@@ -381,18 +386,9 @@
          END SELECT
          !
-         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
-            END_2D
-         ENDIF
-
-         IF( ln_sco ) THEN
-            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
-         ELSE
-            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
-         ENDIF
+         !                                   !==  horizontal fluxes and potential vorticity ==!
+         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
+         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
+         zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
+         !
          !                                   !==  compute and add the vorticity term trend  =!
          DO_2D( 0, 0, 0, 0 )
@@ -458,4 +454,9 @@
                   &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
+            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
+               DO_2D( 1, 0, 1, 0 )
+                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
+               END_2D
+            ENDIF
          CASE ( np_MET )                           !* metric term
             DO_2D( 1, 0, 1, 0 )
@@ -468,4 +469,9 @@
                   &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END_2D
+            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity (NOT the Coriolis term)
+               DO_2D( 1, 0, 1, 0 )
+                  zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
+               END_2D
+            ENDIF
          CASE ( np_CME )                           !* Coriolis + metric
             DO_2D( 1, 0, 1, 0 )
@@ -477,18 +483,10 @@
          END SELECT
          !
-         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
-            END_2D
-         ENDIF
-         !
-         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
-            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
-         ELSE
-            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
-         ENDIF
+         !
+         !                                   !==  horizontal fluxes and potential vorticity ==!
+         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
+         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
+         zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
+         !
          !                                   !==  compute and add the vorticity term trend  =!
          DO_2D( 0, 0, 0, 0 )
@@ -574,4 +572,5 @@
          !
          SELECT CASE( kvor )                 !==  vorticity considered  ==!
+         !
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
             DO_2D( 1, 0, 1, 0 )
@@ -583,4 +582,9 @@
                   &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
             END_2D
+            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
+               DO_2D( 1, 0, 1, 0 )
+                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
+               END_2D
+            ENDIF
          CASE ( np_MET )                           !* metric term
             DO_2D( 1, 0, 1, 0 )
@@ -594,4 +598,9 @@
                   &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
             END_2D
+            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
+               DO_2D( 1, 0, 1, 0 )
+                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 
+               END_2D
+            ENDIF
          CASE ( np_CME )                           !* Coriolis + metric
             DO_2D( 1, 0, 1, 0 )
@@ -602,37 +611,26 @@
             CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
          END SELECT
-         !
-         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
-            END_2D
-         ENDIF
+         !                                             ! ===============
       END DO                                           !   End of slab
-         !
+      !                                                ! ===============
+      !
       CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
-
+      !
+      !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
+         !                                             ! ===============
          !
          !                                   !==  horizontal fluxes  ==!
          zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
          zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
-
+         !
          !                                   !==  compute and add the vorticity term trend  =!
-         jj = 2
-         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
-         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
-               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
-               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
-               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
-               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
-         END DO
-         DO jj = 3, jpj
-            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
-               ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
-               ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
-               ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
-               ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
-            END DO
-         END DO
+         DO_2D( 0, 1, 0, 1 )
+            ztne(ji,jj) = zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk)
+            ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk)
+            ztse(ji,jj) = zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk)
+            ztsw(ji,jj) = zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk)
+         END_2D
+         !
          DO_2D( 0, 0, 0, 0 )
             zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
@@ -704,4 +702,9 @@
                   &          * r1_e1e2f(ji,jj)
             END_2D
+            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
+               DO_2D( 1, 0, 1, 0 )
+                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
+               END_2D
+            ENDIF
          CASE ( np_MET )                           !* metric term
             DO_2D( 1, 0, 1, 0 )
@@ -715,4 +718,9 @@
                   &                         * r1_e1e2f(ji,jj)    )
             END_2D
+            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity
+               DO_2D( 1, 0, 1, 0 )
+                  zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 
+               END_2D
+            ENDIF
          CASE ( np_CME )                           !* Coriolis + metric
             DO_2D( 1, 0, 1, 0 )
@@ -724,38 +732,27 @@
          END SELECT
          !
-         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
-            DO_2D( 1, 0, 1, 0 )
-               zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
-            END_2D
-         ENDIF
-      END DO
+         !                                             ! ===============
+      END DO                                           !   End of slab
+      !                                                ! ===============
       !
       CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
       !
+      !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
-
-      !                                   !==  horizontal fluxes  ==!
+         !                                             ! ===============
+         !
+         !                                   !==  horizontal fluxes  ==!
          zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
          zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
 
          !                                   !==  compute and add the vorticity term trend  =!
-         jj = 2
-         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
-         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
-               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
-               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
-               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
-               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
-               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
-         END DO
-         DO jj = 3, jpj
-            DO ji = 2, jpi   ! vector opt. ok because we start at jj = 3
-               z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
-               ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
-               ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
-               ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
-               ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
-            END DO
-         END DO
+         DO_2D( 0, 1, 0, 1 )
+            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
+            ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
+            ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
+            ztse(ji,jj) = ( zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
+            ztsw(ji,jj) = ( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) ) * z1_e3t
+         END_2D
+         !
          DO_2D( 0, 0, 0, 0 )
             zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
@@ -809,6 +806,4 @@
          WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
       ENDIF
-
-      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
 
 !!gm  this should be removed when choosing a unique strategy for fmask at the coast
