Index: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 5967)
+++ branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 5972)
@@ -16,4 +16,5 @@
    !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
    !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 
+   !!             -   ! 2014-06  (G. Madec) suppression of velocity curl from in-core memory
    !!----------------------------------------------------------------------
 
@@ -22,5 +23,4 @@
    !!       vor_ens  : enstrophy conserving scheme       (ln_dynvor_ens=T)
    !!       vor_ene  : energy conserving scheme          (ln_dynvor_ene=T)
-   !!       vor_mix  : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
    !!       vor_een  : energy and enstrophy conserving   (ln_dynvor_een=T)
    !!   dyn_vor_init : set and control of the different vorticity option
@@ -32,4 +32,6 @@
    USE trd_oce        ! trends: ocean variables
    USE trddyn         ! trend manager: dynamics
+   USE c1d            ! 1D vertical configuration
+   !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE prtctl         ! Print control
@@ -44,23 +46,38 @@
 
    PUBLIC   dyn_vor        ! routine called by step.F90
-   PUBLIC   dyn_vor_init   ! routine called by opa.F90
+   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90
 
    !                                   !!* Namelist namdyn_vor: vorticity term
-   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: energy conserving scheme
-   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme
-   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme
-   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and enstrophy conserving scheme
-   LOGICAL, PUBLIC ::   ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation)
-
-   INTEGER ::   nvor = 0   ! type of vorticity trend used
-   INTEGER ::   ncor = 1   ! coriolis
-   INTEGER ::   nrvm = 2   ! =2 relative vorticity ; =3 metric term
-   INTEGER ::   ntot = 4   ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
-
+   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: energy conserving scheme    (ENE)
+   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme (ENS)
+   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                (MIX)
+   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and 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_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
+
+   INTEGER ::   nvor_scheme        ! choice of the type of advection scheme
+   !                               ! associated indices:
+   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme
+   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 2   ! ENS scheme
+   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 3   ! MIX scheme
+   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
+
+   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity 
+   !                               ! associated indices:
+   INTEGER, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
+   INTEGER, PARAMETER ::   np_RVO = 2         ! relative vorticity
+   INTEGER, PARAMETER ::   np_MET = 3         ! metric term
+   INTEGER, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity)
+   INTEGER, PARAMETER ::   np_CME = 5         ! Coriolis + metric term
+   
+   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
+   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
+   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
+   
    !! * Substitutions
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -87,26 +104,11 @@
       IF( l_trddyn )   CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
       !
-      !                                          ! vorticity term 
-      SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend
-      !
-      CASE ( -1 )                                      ! esopa: test all possibility with control print
-         CALL vor_ene( kt, ntot, ua, va )
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL vor_ens( kt, ntot, ua, va )
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL vor_mix( kt )
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         CALL vor_een( kt, ntot, ua, va )
-         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
-            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         !
-      CASE ( 0 )                                       ! energy conserving scheme
-         IF( l_trddyn )   THEN
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, nrvm, ua, va )                ! relative vorticity or metric trend
+      SELECT CASE ( nvor_scheme )               !==  vorticity trend added to the general trend  ==!
+      !
+      CASE ( np_ENE )                                 !* energy conserving scheme
+         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
+            ztrdu(:,:,:) = ua(:,:,:)
+            ztrdv(:,:,:) = va(:,:,:)
+            CALL vor_ene( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -114,17 +116,17 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend
+            CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity
-         ENDIF
-         !
-      CASE ( 1 )                                       ! enstrophy conserving scheme
-         IF( l_trddyn )   THEN    
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend
+            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity trend
+         ENDIF
+         !
+      CASE ( np_ENS )                                 !* enstrophy conserving scheme
+         IF( l_trddyn ) THEN                                ! trend diagnostics: splitthe trend in two    
+            ztrdu(:,:,:) = ua(:,:,:)
+            ztrdv(:,:,:) = va(:,:,:)
+            CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -132,17 +134,17 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, ncor, ua, va )                ! planetary vorticity trend
+            CALL vor_ens( kt, ncor, ua, va )                      ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity
-         ENDIF
-         !
-      CASE ( 2 )                                       ! mixed ene-ens scheme
-         IF( l_trddyn )   THEN
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens)
+            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity trend
+         ENDIF
+         !
+      CASE ( np_MIX )                                 !* mixed ene-ens scheme
+         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
+            ztrdu(:,:,:) = ua(:,:,:)
+            ztrdv(:,:,:) = va(:,:,:)
+            CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend (ens)
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -150,17 +152,18 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene)
+            CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend (ene)
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene)
-         ENDIF
-         !
-      CASE ( 3 )                                       ! energy and enstrophy conserving scheme
-         IF( l_trddyn )   THEN
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_een( kt, nrvm, ua, va )                ! relative vorticity or metric trend
+            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens)
+            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene)
+        ENDIF
+         !
+      CASE ( np_EEN )                                 !* energy and enstrophy conserving scheme
+         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two
+            ztrdu(:,:,:) = ua(:,:,:)
+            ztrdv(:,:,:) = va(:,:,:)
+            CALL vor_een( kt, nrvm, ua, va )                      ! relative vorticity or metric trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
@@ -168,10 +171,10 @@
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-            CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend
+            CALL vor_een( kt, ncor, ua, va )                      ! planetary vorticity trend
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
-            CALL vor_een( kt, ntot, ua, va )                ! total vorticity
+            CALL vor_een( kt, ntot, ua, va )                ! total vorticity trend
          ENDIF
          !
@@ -197,15 +200,11 @@
       !!
       !! ** Method  :   Trend evaluated using now fields (centered in time) 
-      !!      and the Sadourny (1975) flux form formulation : conserves the
-      !!      horizontal kinetic energy.
-      !!      The trend of the vorticity term is given by:
-      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
-      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ]
-      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ]
-      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
-      !!          voru = 1/e1u  mj-1[ (rotn+f)  mi(e1v vn) ]
-      !!          vorv = 1/e2v  mi-1[ (rotn+f)  mj(e2u un) ]
-      !!      Add this trend to the general momentum trend (ua,va):
-      !!          (ua,va) = (ua,va) + ( voru , vorv )
+      !!       and the Sadourny (1975) flux form formulation : conserves the
+      !!       horizontal kinetic energy.
+      !!         The general trend of momentum is increased due to the vorticity 
+      !!       term which is given by:
+      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v vn) ]
+      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u un) ]
+      !!       where rvor is the relative vorticity
       !!
       !! ** Action : - Update (ua,va) with the now vorticity term trend
@@ -213,5 +212,4 @@
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
       !!----------------------------------------------------------------------
-      !
       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
@@ -220,7 +218,7 @@
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
       !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zx1, zy1, zfact2, zx2, zy2   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz
+      INTEGER  ::   ji, jj, jk           ! dummy loop indices
+      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz   ! 2D workspace
       !!----------------------------------------------------------------------
       !
@@ -234,37 +232,55 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF
-
-      zfact2 = 0.5 * 0.5      ! Local constant initialization
-
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
+      !
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
          !
-         ! Potential vorticity and horizontal fluxes
-         ! -----------------------------------------
-         SELECT CASE( kvor )      ! vorticity considered
-         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
-         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
-         CASE ( 3 )                                                ! metric term
+         SELECT CASE( kvor )                 !==  vorticity considered  ==!
+         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
+            zwz(:,:) = ff(:,:) 
+         CASE ( np_RVO )                           !* relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_MET )                           !* metric term
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
-               END DO
-            END DO
-         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
-         CASE ( 5 )                                                ! total (coriolis + metric)
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
-                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               &
-                       &       )
-               END DO
-            END DO
+                       &     * 0.5 * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CRV )                           !* Coriolis + relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                     &                   * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CME )                           !* Coriolis + metric
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = ff(ji,jj)                                                                        &
+                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                       &     * 0.5 * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE DEFAULT                                             ! error
+            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
          END SELECT
+         !
+         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
+               END DO
+            END DO
+         ENDIF
 
          IF( ln_sco ) THEN
@@ -276,7 +292,5 @@
             zwy(:,:) = e1v(:,:) * vn(:,:,jk)
          ENDIF
-
-         ! Compute and add the vorticity term trend
-         ! ----------------------------------------
+         !                                   !==  compute and add the vorticity term trend  =!
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -285,6 +299,6 @@
                zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
                zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
+               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
+               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
             END DO  
          END DO  
@@ -299,123 +313,4 @@
 
 
-   SUBROUTINE vor_mix( kt )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE vor_mix  ***
-      !!
-      !! ** Purpose :   Compute the now total vorticity trend and add it to
-      !!      the general trend of the momentum equation.
-      !!
-      !! ** Method  :   Trend evaluated using now fields (centered in time)
-      !!      Mixte formulation : conserves the potential enstrophy of a hori-
-      !!      zontally non-divergent flow for (rotzu x uh), the relative vor-
-      !!      ticity term and the horizontal kinetic energy for (f x uh), the
-      !!      coriolis term. the now trend of the vorticity term is given by:
-      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
-      !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
-      !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ]
-      !!          vorv = 1/e2v  mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
-      !!              +1/e2v  mi-1[ f/e3f          mj(e2u*e3u un) ]
-      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
-      !!          voru = 1/e1u  mj-1(rotn) mj-1[ mi(e1v vn) ]
-      !!              +1/e1u  mj-1[ f          mi(e1v vn) ]
-      !!          vorv = 1/e2v  mi-1(rotn) mi-1[ mj(e2u un) ]
-      !!              +1/e2v  mi-1[ f          mj(e2u un) ]
-      !!      Add this now trend to the general momentum trend (ua,va):
-      !!          (ua,va) = (ua,va) + ( voru , vorv )
-      !!
-      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
-      !!
-      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
-      !!----------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt   ! ocean timestep index
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! local scalars
-      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('vor_mix')
-      !
-      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww ) 
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-
-      zfact1 = 0.5 * 0.25      ! Local constant initialization
-      zfact2 = 0.5 * 0.5
-
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         !
-         ! Relative and planetary potential vorticity and horizontal fluxes
-         ! ----------------------------------------------------------------
-         IF( ln_sco ) THEN        
-            IF( ln_dynadv_vec ) THEN
-               zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
-            ELSE                       
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1   ! vector opt.
-                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
-                  END DO
-               END DO
-            ENDIF
-            zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk)
-            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
-            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
-         ELSE
-            IF( ln_dynadv_vec ) THEN
-               zww(:,:) = rotn(:,:,jk)
-            ELSE                       
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1   ! vector opt.
-                     zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                        &       * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
-                  END DO
-               END DO
-            ENDIF
-            zwz(:,:) = ff (:,:)
-            zwx(:,:) = e2u(:,:) * un(:,:,jk)
-            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
-         ENDIF
-
-         ! Compute and add the vorticity term trend
-         ! ----------------------------------------
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
-               zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
-               zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
-               zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
-               ! enstrophy conserving formulation for relative vorticity term
-               zua = zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
-               zva =-zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1 + zx2 )
-               ! energy conserving formulation for planetary vorticity term
-               zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-               zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
-               ! mixed vorticity trend added to the momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
-               va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
-            END DO  
-         END DO  
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('vor_mix')
-      !
-   END SUBROUTINE vor_mix
-
-
    SUBROUTINE vor_ens( kt, kvor, pua, pva )
       !!----------------------------------------------------------------------
@@ -429,10 +324,6 @@
       !!      potential enstrophy of a horizontally non-divergent flow. the
       !!      trend of the vorticity term is given by:
-      !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative:
-      !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
-      !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
-      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
-      !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ]
-      !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ]
+      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
+      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
       !!      Add this trend to the general momentum trend (ua,va):
       !!          (ua,va) = (ua,va) + ( voru , vorv )
@@ -442,5 +333,4 @@
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
       !!----------------------------------------------------------------------
-      !
       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
@@ -449,7 +339,7 @@
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
       !
-      INTEGER  ::   ji, jj, jk           ! dummy loop indices
-      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zuav, zvau   ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz, zww   ! 2D workspace
       !!----------------------------------------------------------------------
       !
@@ -463,63 +353,70 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF
-
-      zfact1 = 0.5 * 0.25      ! Local constant initialization
-
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
          !
-         ! Potential vorticity and horizontal fluxes
-         ! -----------------------------------------
-         SELECT CASE( kvor )      ! vorticity considered
-         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis)
-         CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity
-         CASE ( 3 )                                                ! metric term
+         SELECT CASE( kvor )                 !==  vorticity considered  ==!
+         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
+            zwz(:,:) = ff(:,:) 
+         CASE ( np_RVO )                           !* relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_MET )                           !* metric term
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
-               END DO
-            END DO
-         CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity)
-         CASE ( 5 )                                                ! total (coriolis + metric)
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
-                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
-                       &       )
-               END DO
-            END DO
+                       &     * 0.5 * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CRV )                           !* Coriolis + relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                     &                   * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CME )                           !* Coriolis + metric
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = ff(ji,jj)                                                                       &
+                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                       &     * 0.5 * r1_e1e2f(ji,jj)
+               END DO
+            END DO
+         CASE DEFAULT                                             ! error
+            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
          END SELECT
          !
-         IF( ln_sco ) THEN
-            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop 
-               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
-                  zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
-                  zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
-                  zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
-               END DO
-            END DO
+         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==!
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
+               END DO
+            END DO
+         ENDIF
+         !
+         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
+            zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
+            zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
+            zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
          ELSE
-            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop 
-               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
-                  zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
-                  zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         ! Compute and add the vorticity term trend
-         ! ----------------------------------------
+            zwx(:,:) = e2u(:,:) * un(:,:,jk)
+            zwy(:,:) = e1v(:,:) * vn(:,:,jk)
+         ENDIF
+         !                                   !==  compute and add the vorticity term trend  =!
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zuav = zfact1 / e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
-                  &                         + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
-               zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
-                  &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
+               zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
+                  &                       + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) )
+               zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
+                  &                       + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) )
                pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
                pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
@@ -553,32 +450,22 @@
       !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
       !!----------------------------------------------------------------------
-      !
       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index
-      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ;
-      !                                                           ! =nrvm (relative vorticity or metric)
+      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric)
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
-      !!
-      INTEGER  ::   ji, jj, jk                                    ! dummy loop indices
-      INTEGER  ::   ierr                                          ! local integer
-      REAL(wp) ::   zfac12, zua, zva                              ! local scalars
-      REAL(wp) ::   zmsk, ze3                                     ! local scalars
-      !                                                           !  3D workspace 
-      REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz
-      REAL(wp), POINTER    , DIMENSION(:,:  )         :: ztnw, ztne, ztsw, ztse
-#if defined key_vvl
-      REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T)
-#else
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all
-#endif
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ierr         ! local integer
+      REAL(wp) ::   zua, zva     ! local scalars
+      REAL(wp) ::   zmsk, ze3    ! local scalars
+      !
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zwx, zwy, zwz, z1_e3f
+      REAL(wp), POINTER, DIMENSION(:,:)   :: ztnw, ztne, ztsw, ztse
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('vor_een')
       !
-      CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        ) 
-      CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
-#if defined key_vvl
-      CALL wrk_alloc( jpi, jpj, jpk, ze3f                   )
-#endif
+      CALL wrk_alloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
+      CALL wrk_alloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
       !
       IF( kt == nit000 ) THEN
@@ -586,90 +473,97 @@
          IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-#if ! defined key_vvl
-         IF( .NOT.ALLOCATED(ze3f) ) THEN
-            ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )
-            IF( lk_mpp    )   CALL mpp_sum ( ierr )
-            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
-         ENDIF
-         ze3f(:,:,:) = 0._wp
-#endif
       ENDIF
-
-      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points)
-
-         IF( ln_dynvor_een_old ) THEN ! original formulation
-            DO jk = 1, jpk
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
-                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = 4.0_wp / ze3
-                  END DO
-               END DO
-            END DO
-         ELSE ! new formulation from NEMO 3.6
-            DO jk = 1, jpk
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*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( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3
-                  END DO
-               END DO
-            END DO
-         ENDIF
-
-         CALL lbc_lnk( ze3f, 'F', 1. )
-      ENDIF
-
-      zfac12 = 1._wp / 12._wp    ! Local constant initialization
-
-      
-!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
+      !
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
-         
-         ! Potential vorticity and horizontal fluxes
-         ! -----------------------------------------
-         SELECT CASE( kvor )      ! vorticity considered
-         CASE ( 1 )                                                ! planetary vorticity (Coriolis)
-            zwz(:,:) = ff(:,:)      * ze3f(:,:,jk)
-         CASE ( 2 )                                                ! relative  vorticity
-            zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
-         CASE ( 3 )                                                ! metric term
+         !
+         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point
+         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
+                     &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
+                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4.0_wp / ze3
+                  ELSE                      ;   z1_e3f(ji,jj) = 0.0_wp
+                  ENDIF
+               END DO
+            END DO
+         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
+                     &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*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( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3
+                  ELSE                      ;   z1_e3f(ji,jj) = 0.0_wp
+                  ENDIF
+               END DO
+            END DO
+         END SELECT
+         !
+         SELECT CASE( kvor )                 !==  vorticity considered  ==!
+         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_RVO )                           !* relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                     &       * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_MET )                           !* metric term
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
-               END DO
-            END DO
-            CALL lbc_lnk( zwz, 'F', 1. )
-        CASE ( 4 )                                                ! total (relative + planetary vorticity)
-            zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
-         CASE ( 5 )                                                ! total (coriolis + metric)
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( ff (ji,jj)                                                                       &
-                       &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
-                       &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
-                       &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                &
-                       &       ) * ze3f(ji,jj,jk)
-               END DO
-            END DO
-            CALL lbc_lnk( zwz, 'F', 1. )
+                       &     * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CRV )                           !* Coriolis + relative vorticity
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
+                     &                         - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) &
+                     &                      * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj)
+               END DO
+            END DO
+         CASE ( np_CME )                           !* Coriolis + metric
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = (  ff(ji,jj)                                                                        &
+                       &        + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &
+                       &            - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &
+                       &        * 0.5 * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
+               END DO
+            END DO
+         CASE DEFAULT                                             ! error
+            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
          END SELECT
-
+         !
+         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==!
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
+               END DO
+            END DO
+         ENDIF
+         !
+         CALL lbc_lnk( zwz, 'F', 1. )
+         !
+         !                                   !==  horizontal fluxes  ==!
          zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
          zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
 
-         ! Compute and add the vorticity term trend
-         ! ----------------------------------------
+         !                                   !==  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   
+         DO ji = 2, jpi          ! split in 2 parts due to vector opt.
                ztne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
                ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
@@ -687,8 +581,8 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zua = + zfac12 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
-                  &                           + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
-               zva = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
-                  &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
+               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
+                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
+               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
+                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
                pua(ji,jj,jk) = pua(ji,jj,jk) + zua
                pva(ji,jj,jk) = pva(ji,jj,jk) + zva
@@ -698,9 +592,7 @@
       END DO                                           !   End of slab
       !                                                ! ===============
-      CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        ) 
-      CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
-#if defined key_vvl
-      CALL wrk_dealloc( jpi, jpj, jpk, ze3f                   )
-#endif
+      !
+      CALL wrk_dealloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f ) 
+      CALL wrk_dealloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('vor_een')
@@ -720,5 +612,5 @@
       INTEGER ::   ios             ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old
+      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk
       !!----------------------------------------------------------------------
 
@@ -737,13 +629,17 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
-         WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
-         WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens
-         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
-         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
-         WRITE(numout,*) '           enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old
+         WRITE(numout,*) '           energy    conserving scheme                    ln_dynvor_ene = ', ln_dynvor_ene
+         WRITE(numout,*) '           enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
+         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
+         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,*) '           masked (=1) or unmasked(=0) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
       ENDIF
 
+!!gm  this should be removed when choosing a unique strategy for fmask at the coast
       ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
       ! at angles with three ocean points and one land point
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '           namlbc: change fmask value in the angles (T)   ln_vorlat = ', ln_vorlat
       IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
          DO jk = 1, jpk
@@ -759,42 +655,32 @@
           !
       ENDIF
-
-      ioptio = 0                     ! Control of vorticity scheme options
-      IF( ln_dynvor_ene )   ioptio = ioptio + 1
-      IF( ln_dynvor_ens )   ioptio = ioptio + 1
-      IF( ln_dynvor_mix )   ioptio = ioptio + 1
-      IF( ln_dynvor_een )   ioptio = ioptio + 1
-      IF( ln_dynvor_een_old )   ioptio = ioptio + 1
-      IF( lk_esopa      )   ioptio =          1
-
-      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
-
-      !                              ! Set nvor (type of scheme for vorticity)
-      IF( ln_dynvor_ene )   nvor =  0
-      IF( ln_dynvor_ens )   nvor =  1
-      IF( ln_dynvor_mix )   nvor =  2
-      IF( ln_dynvor_een .or. ln_dynvor_een_old )   nvor =  3
-      IF( lk_esopa      )   nvor = -1
-      
-      !                              ! Set ncor, nrvm, ntot (type of vorticity)
-      IF(lwp) WRITE(numout,*)
-      ncor = 1
+!!gm end
+
+      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
+      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENE   ;   ENDIF
+      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENS   ;   ENDIF
+      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_MIX   ;   ENDIF
+      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_EEN   ;   ENDIF
+      !
+      IF( ( ioptio /= 1).AND.( .NOT.lk_c1d ) ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
+      !                      
+      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
+      ncor = np_COR
       IF( ln_dynadv_vec ) THEN     
          IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
-         nrvm = 2
-         ntot = 4
+         nrvm = np_RVO        ! relative vorticity
+         ntot = np_CRV        ! relative + planetary vorticity
       ELSE                        
          IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
-         nrvm = 3
-         ntot = 5
+         nrvm = np_MET        ! metric term
+         ntot = np_CME        ! Coriolis + metric term
       ENDIF
       
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
-         IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme'
-         IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme'
-         IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme'
-         IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme'
-         IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options'
+         IF( nvor_scheme ==  np_ENE )   WRITE(numout,*) '         vorticity scheme ==>> energy conserving scheme'
+         IF( nvor_scheme ==  np_ENS )   WRITE(numout,*) '         vorticity scheme ==>> enstrophy conserving scheme'
+         IF( nvor_scheme ==  np_MIX )   WRITE(numout,*) '         vorticity scheme ==>> mixed enstrophy/energy conserving scheme'
+         IF( nvor_scheme ==  np_EEN )   WRITE(numout,*) '         vorticity scheme ==>> energy and enstrophy conserving scheme'
       ENDIF
       !
