Index: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 9506)
+++ branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 9528)
@@ -36,4 +36,5 @@
    USE sbcapr          ! surface boundary condition: atmospheric pressure
    USE dynadv    , ONLY: ln_dynadv_vec
+   USE dynvor          ! vortivity scheme indicators
    USE phycst          ! physical constants
    USE dynvor          ! vorticity term
@@ -103,6 +104,7 @@
       ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) )
       !
-      IF( ln_dynvor_een )   ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
-         &                            ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) )
+      IF( ln_dynvor_een .OR. ln_dynvor_eeT )   &
+         &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
+         &               ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) )
          !
       ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj)                    , STAT=ierr(3) )
@@ -149,13 +151,13 @@
       INTEGER  ::   ikbu, iktu, noffset   ! local integers
       INTEGER  ::   ikbv, iktv            !   -      -
-      REAL(wp) ::   r1_2dt_b, z2dt_bf  	   ! local scalars
-      REAL(wp) ::   zx1, zx2, zu_spg, zhura  !   -      -
-      REAL(wp) ::   zy1, zy2, zv_spg, zhvra  !   -      -
-      REAL(wp) ::   za0, za1, za2, za3       !   -      -
-      REAL(wp) ::   zmdi, zztmp              !   -      -
+      REAL(wp) ::   r1_2dt_b, z2dt_bf  	          ! local scalars
+      REAL(wp) ::   zx1, zx2, zu_spg, zhura, z1_hu  !   -      -
+      REAL(wp) ::   zy1, zy2, zv_spg, zhvra, z1_hv  !   -      -
+      REAL(wp) ::   za0, za1, za2, za3              !   -      -
+      REAL(wp) ::   zmdi, zztmp            , z1_ht  !   -      -
       REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf
       REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc
       REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv
-      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e
+      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e
       REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e
       REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
@@ -237,5 +239,7 @@
       !
       IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN
-         IF( ln_dynvor_een ) THEN               !==  EEN scheme  ==!
+         !
+         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 ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
@@ -250,8 +254,15 @@
                DO jj = 1, jpjm1
                   DO ji = 1, jpim1
-                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                     &
-                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   )                   &
-                        &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
-                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
+!!gm bug for ocean cavities
+!                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                     &
+!                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   )                   &
+!                        &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    &
+!                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) )
+!! replace by :
+                     zwz(ji,jj) =             (  ht_n  (ji  ,jj+1) + ht_n  (ji+1,jj+1)      &
+                        &                      + ht_n  (ji  ,jj  ) + ht_n  (ji+1,jj  )  )   &
+                        &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      &
+                        &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   )
+!!gm end
                      IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
                   END DO
@@ -270,5 +281,18 @@
             END DO
             !
-         ELSE                                !== all other schemes (ENE, ENS, MIX)
+         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 jj = 2, jpj
+               DO ji = 2, jpi
+                  z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )
+                  ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
+                  ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht
+                  ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
+                  ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht
+               END DO
+            END DO
+            !
+         CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT !
+            !
             zwz(:,:) = 0._wp
             zhf(:,:) = 0._wp
@@ -280,5 +304,5 @@
 !!gm 
 !!            
-              IF( .NOT.ln_sco ) THEN
+            IF( .NOT.ln_sco ) THEN
   
    !!gm  agree the JC comment  : this should be done in a much clear way
@@ -290,25 +314,33 @@
    !              ENDIF
    !              zhf(:,:) = gdepw_0(:,:,jk+1)
-               ELSE
-                 !zhf(:,:) = hbatf(:,:)
-                 DO jj = 1, jpjm1
-                   DO ji = 1, jpim1
-                     zhf(ji,jj) = MAX( 0._wp,                                &
-                                & ( ht_0(ji  ,jj  )*tmask(ji  ,jj  ,1) +     &
-                                &   ht_0(ji+1,jj  )*tmask(ji+1,jj  ,1) +     &
-                                &   ht_0(ji  ,jj+1)*tmask(ji  ,jj+1,1) +     &
-                                &   ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) /   &
-                                & ( tmask(ji  ,jj  ,1) + tmask(ji+1,jj  ,1) +&
-                                &   tmask(ji  ,jj+1,1) + tmask(ji+1,jj+1,1) +&
-                                &   rsmall  ) )
-                   END DO
-                 END DO
-              END IF
-  
-              DO jj = 1, jpjm1
-                 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
-              END DO
+               !
+            ELSE
+               !
+               !zhf(:,:) = hbatf(:,:)
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1
+!!gm Bug here in case of ocean cavities, further more ht_0 is masked by definition ==>> remove mask multiplication 
+!                     zhf(ji,jj) = MAX( 0._wp,                                &
+!                                & ( ht_0(ji  ,jj  )*tmask(ji  ,jj  ,1) +     &
+!                                &   ht_0(ji+1,jj  )*tmask(ji+1,jj  ,1) +     &
+!                                &   ht_0(ji  ,jj+1)*tmask(ji  ,jj+1,1) +     &
+!                                &   ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) /   &
+!                                & ( tmask(ji  ,jj  ,1) + tmask(ji+1,jj  ,1) +&
+!                                &   tmask(ji  ,jj+1,1) + tmask(ji+1,jj+1,1) +&
+!                                &   rsmall  ) )
+!!gm replaced by :
+                     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  )
 !!gm end
-
+                  END DO
+               END DO
+            ENDIF
+            !
+            DO jj = 1, jpjm1
+               zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))
+            END DO
+            !
             DO jk = 1, jpkm1
                DO jj = 1, jpjm1
@@ -324,5 +356,5 @@
             END DO
             zwz(:,:) = ff_f(:,:) * zwz(:,:)
-         ENDIF
+         END SELECT
       ENDIF
       !
@@ -369,8 +401,23 @@
       !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
       !                                   ! --------------------------------------------------------
+      !
       zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes 
       zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:)
       !
-      IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
+      SELECT CASE( nvor_scheme )
+      CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1   ! vector opt.
+               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj)                    &
+                  &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   &
+                  &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   )
+                  !
+               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj)                    &
+                  &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   & 
+                  &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   ) 
+            END DO  
+         END DO  
+         !         
+      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -385,5 +432,5 @@
          END DO
          !
-      ELSEIF ( ln_dynvor_ens ) THEN                    ! enstrophy conserving scheme
+      CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -397,5 +444,5 @@
          END DO
          !
-      ELSEIF ( ln_dynvor_een ) THEN  ! enstrophy and energy conserving scheme
+      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -411,64 +458,60 @@
          END DO
          !
-      ENDIF 
+      END SELECT
       !
       !                                   !* Right-Hand-Side of the barotropic momentum equation
       !                                   ! ----------------------------------------------------
       IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient
-        IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters
-           DO jj = 2, jpjm1
-              DO ji = 2, jpim1 
-                ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
-                     &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
-                     &    MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
+         IF( ln_wd_il ) THEN                        ! Calculating and applying W/D gravity filters
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1 
+                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
+                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
+                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
+                     &                                                         > rn_wdmin1 + rn_wdmin2
+                  ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( &
+                     &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
+                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
+                  IF(ll_tmp1) THEN
+                     zcpx(ji,jj) = 1.0_wp
+                  ELSEIF(ll_tmp2) THEN
+                     ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
+                     zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
+                                 &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
+                     zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
+                  ELSE
+                     zcpx(ji,jj) = 0._wp
+                  ENDIF
+                  !
+                  ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                &
+                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
+                     &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
                      &                                                       > rn_wdmin1 + rn_wdmin2
-                ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( &
-                     &    MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
-                     &    MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
-   
-                IF(ll_tmp1) THEN
-                  zcpx(ji,jj) = 1.0_wp
-                ELSE IF(ll_tmp2) THEN
-                  ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-                  zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
-                              &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
-                  zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
-
-                ELSE
-                  zcpx(ji,jj) = 0._wp
-                END IF
-         
-                ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                &
-                     &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
-                     &    MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
-                     &                                                       > rn_wdmin1 + rn_wdmin2
-                ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( &
-                     &    MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
-                     &    MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
+                  ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( &
+                     &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
+                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
   
-                IF(ll_tmp1) THEN
-                  zcpy(ji,jj) = 1.0_wp
-                ELSE IF(ll_tmp2) THEN
-                  ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-                  zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
-                              &    / (sshn(ji,jj+1) - sshn(ji,jj  )) )
-                  zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp)
-
-                ELSE
-                  zcpy(ji,jj) = 0._wp
-                END IF
-              END DO
-           END DO
- 
-           DO jj = 2, jpjm1
-              DO ji = 2, jpim1
-                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
-                        &                        * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth
-                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   &
-                        &                        * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth
-
-              END DO
-           END DO
-
+                  IF(ll_tmp1) THEN
+                     zcpy(ji,jj) = 1.0_wp
+                  ELSE IF(ll_tmp2) THEN
+                     ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
+                     zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
+                        &             / (sshn(ji,jj+1) - sshn(ji,jj  )) )
+                     zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
+                  ELSE
+                     zcpy(ji,jj) = 0._wp
+                  ENDIF
+               END DO
+            END DO
+            !
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
+                     &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth
+                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   &
+                     &                          * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth
+               END DO
+            END DO
+            !
          ELSE
             !
@@ -675,14 +718,13 @@
       un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop
       vn_adv(:,:) = 0._wp
-      !                                             ! ==================== !
-
-      IF (ln_wd_dl) THEN
+      !
+      IF( ln_wd_dl ) THEN
          zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary) 
          zvwdmask(:,:) = 0._wp  ! 
-         zuwdav2(:,:) =  0._wp 
-         zvwdav2(:,:) =  0._wp   
+         zuwdav2 (:,:) = 0._wp 
+         zvwdav2 (:,:) = 0._wp   
       END IF 
 
-
+      !                                             ! ==================== !
       DO jn = 1, icycle                             !  sub-time-step loop  !
          !                                          ! ==================== !
@@ -715,5 +757,5 @@
             ! set wetting & drying mask at tracer points for this barotropic sub-step 
             IF ( ln_wd_dl ) THEN 
-
+               !
                IF ( ln_wd_dl_rmp ) THEN 
                   DO jj = 1, jpj                                 
@@ -736,12 +778,11 @@
                         ELSE 
                            ztwdmask(ji,jj) = 0._wp
-                        END IF
+                        ENDIF
                      END DO
                   END DO
-               END IF 
-
-            END IF 
-           
-
+               ENDIF 
+               !
+            ENDIF 
+            !
             DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
                DO ji = 2, fs_jpim1   ! Vector opt.
@@ -756,9 +797,11 @@
             CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp )
             !
-            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
-            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
+            zhup2_e(:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
+            zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:)
+            zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:)
          ELSE
-            zhup2_e (:,:) = hu_n(:,:)
-            zhvp2_e (:,:) = hv_n(:,:)
+            zhup2_e(:,:) = hu_n(:,:)
+            zhvp2_e(:,:) = hv_n(:,:)
+            zhtp2_e(:,:) = ht_n(:,:)
          ENDIF
          !                                                !* after ssh
@@ -795,11 +838,10 @@
          ENDIF
 #endif
-         IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
+         IF( ln_wd_il )   CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
 
          IF ( ln_wd_dl ) THEN 
-
-
-! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 
-
+            !
+            ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 
+            !
             DO jj = 1, jpjm1                                 
                DO ji = 1, jpim1   
@@ -821,7 +863,6 @@
                END DO
             END DO
-
-
-         END IF    
+            !
+         ENDIF    
          
          ! Sum over sub-time-steps to compute advective velocities
@@ -896,6 +937,7 @@
          ENDIF
          !
-         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) &
-          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
+         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:)   &
+            &          + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:)
+         
          IF( ln_wd_il ) THEN                   ! Calculating and applying W/D gravity filters
            DO jj = 2, jpjm1
@@ -917,6 +959,6 @@
                 ELSE
                   zcpx(ji,jj) = 0._wp
-                END IF
-         
+                ENDIF
+                !
                 ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                &
                      &    MAX(    -ht_0(ji,jj)               ,    -ht_0(ji,jj+1) ) .AND.            &
@@ -929,5 +971,5 @@
                 IF(ll_tmp1) THEN
                   zcpy(ji,jj) = 1.0_wp
-                ELSE IF(ll_tmp2) THEN
+                ELSEIF(ll_tmp2) THEN
                   ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here
                   zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +     ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) &
@@ -935,8 +977,8 @@
                 ELSE
                   zcpy(ji,jj) = 0._wp
-                END IF
+                ENDIF
               END DO
            END DO
-         END IF
+         ENDIF
          !
          ! Compute associated depths at U and V points:
@@ -955,5 +997,5 @@
                END DO
             END DO
-
+            !
          ENDIF
          !
@@ -965,5 +1007,23 @@
          ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)
          !
-         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN     !==  energy conserving or mixed scheme  ==!
+         SELECT CASE( nvor_scheme )
+         CASE( np_ENT )             ! energy conserving scheme (t-point)
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1   ! vector opt.
+
+               z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) )
+               z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) )
+            
+               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                   &
+                  &               * (  e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) )   &
+                  &                  + e1e2t(ji  ,jj)*zhtp2_e(ji  ,jj)*ff_t(ji  ,jj) * ( va_e(ji  ,jj) + va_e(ji  ,jj-1) )   )
+                  !
+               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
+                  &               * (  e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) )   & 
+                  &                  + e1e2t(ji,jj  )*zhtp2_e(ji,jj  )*ff_t(ji,jj  ) * ( ua_e(ji,jj  ) + ua_e(ji-1,jj  ) )   ) 
+            END DO  
+         END DO  
+         !         
+         CASE( np_ENE, np_MIX )     ! energy conserving scheme (f-point)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -977,5 +1037,5 @@
             END DO
             !
-         ELSEIF ( ln_dynvor_ens ) THEN                   !==  enstrophy conserving scheme  ==!
+         CASE( np_ENS )             ! enstrophy conserving scheme (f-point)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -989,19 +1049,19 @@
             END DO
             !
-         ELSEIF ( ln_dynvor_een ) THEN                   !==  energy and enstrophy conserving scheme  ==!
+         CASE( np_EET , np_EEN )   ! energy & enstrophy scheme (using e3t or e3f)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
-                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
-                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
-                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
-                  zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
-                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
-                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
-                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
+                  zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  )  &
+                     &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  )  &
+                     &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1)  & 
+                     &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1)  )
+                  zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1)  & 
+                     &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1)  &
+                     &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  )  & 
+                     &                                       + ftne(ji,jj  ) * zwx(ji  ,jj  )  )
                END DO
             END DO
             ! 
-         ENDIF
+         END SELECT
          !
          ! Add tidal astronomical forcing if defined
@@ -1341,8 +1401,6 @@
       !! ** Purpose : Read or write time-splitting arrays in restart file
       !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
-      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
-      !
-      INTEGER ::   id1, id2   ! local integers
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
       !!----------------------------------------------------------------------
       !
