Index: /trunk/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 107)
+++ /trunk/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 108)
@@ -10,4 +10,5 @@
    !!   dyn_vor_energy   : energy conserving scheme          (ln_dynvor_ene=T)
    !!   dyn_vor_mixed    : mixed enstrophy/energy conserving (ln_dynvor_mix=T)
+   !!   dyn_vor_ene_ens  : energy and enstrophy conserving   (ln_dynvor_een=T)
    !!   dyn_vor_ctl      : control of the different vorticity option
    !!----------------------------------------------------------------------
@@ -17,4 +18,5 @@
    USE in_out_manager ! I/O manager
    USE trddyn_oce     ! ocean momentum trends
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
 
    IMPLICIT NONE
@@ -25,4 +27,5 @@
    PUBLIC dyn_vor_energy         ! routine called by step.F90
    PUBLIC dyn_vor_mixed          ! routine called by step.F90
+   PUBLIC dyn_vor_ene_ens        ! routine called by step.F90
    PUBLIC dyn_vor_ctl            ! routine called by step.F90
 
@@ -31,4 +34,5 @@
    LOGICAL, PUBLIC ::   ln_dynvor_ens = .TRUE.    !: enstrophy conserving scheme
    LOGICAL, PUBLIC ::   ln_dynvor_mix = .FALSE.   !: mixed scheme
+   LOGICAL, PUBLIC ::   ln_dynvor_een = .TRUE.    !: energy and enstrophy conserving scheme
 
    !! * Substitutions
@@ -82,5 +86,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   &
          zwx, zwy, zwz                      ! temporary workspace
-#if defined key_trddyn
+#if defined key_trddyn || defined key_trd_vor
       REAL(wp) ::   &
          zcu, zcv, zce3                     !    "         "
@@ -125,5 +129,5 @@
                ua(ji,jj,jk) = ua(ji,jj,jk) + zua
                va(ji,jj,jk) = va(ji,jj,jk) + zva
-#   if defined key_trddyn
+#   if defined key_trddyn || defined key_trd_vor
 #      if defined key_s_coord
                zce3= ff(ji,jj) / fse3f(ji,jj,jk)
@@ -248,5 +252,5 @@
                ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
                va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
-#   if defined key_trddyn
+#   if defined key_trddyn || defined key_trd_vor
                utrd(ji,jj,jk,3) = zua
                vtrd(ji,jj,jk,3) = zva
@@ -313,5 +317,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
          zwz                           ! temporary workspace
-#   if defined key_trddyn
+#   if defined key_trddyn || defined key_trd_vor
       REAL(wp) ::   &
          zcu, zcv, zce3                ! temporary scalars 
@@ -368,5 +372,5 @@
                va(ji,jj,jk) = va(ji,jj,jk) + zva
 
-#   if defined key_trddyn
+#   if defined key_trddyn || defined key_trd_vor
 #      if defined key_s_coord
                zce3 = ff(ji,jj) / fse3f(ji,jj,jk)
@@ -404,4 +408,138 @@
 
 
+   SUBROUTINE dyn_vor_ene_ens( kt )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dyn_vor_ene_ens  ***
+      !!
+      !! ** 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) 
+      !!      and the Arakawa and Lamb (19XX) flux form formulation : conserves 
+      !!      both the horizontal kinetic energy and the potential enstrophy
+      !!      when horizontal divergence is zero.
+      !!      The trend of the vorticity term is given by:
+      !!       * s-coordinate (lk_sco=T), the e3. are inside the derivatives:
+      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
+      !!      Add this trend to the general momentum trend (ua,va):
+      !!          (ua,va) = (ua,va) + ( voru , vorv )
+      !!
+      !! ** Action : - Update (ua,va) with the now vorticity term trend
+      !!             - save the trends in (utrd,vtrd) in 2 parts (relative
+      !!               and planetary vorticity trends) ('key_trddyn')
+      !!
+      !! References :
+      !!      Arakawa and Lamb 19XX, ???
+      !! History :
+      !!   5.0  !  04-02  (G. Madec)  Original code
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
+
+      !! * Local declarations
+      INTEGER ::   ji, jj, jk              ! dummy loop indices
+      REAL(wp) ::   &
+         zfac12, zua, zva                  ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   &
+         zwx, zwy, zwz,                 &  ! temporary workspace
+         ztnw, ztne, ztsw, ztse            !    "           "
+#if defined key_trddyn || defined key_trd_vor
+      REAL(wp) ::   &
+         zcu, zcv                          !    "         "
+      REAL(wp), DIMENSION(jpi,jpj) ::   &
+         zcor                              ! potential planetary vorticity (f/e3)
+#endif
+      REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE ::   &
+         ze3f
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_vor_ene_ens : vorticity term: energy and enstrophy conserving scheme'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
+
+         DO jk = 1, jpk
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  ze3f(ji,jj,jk) = ( 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) ) * 0.25
+!!!               ze3f(ji,jj,jk) = MAX( ze3f(ji,jj,jk) , 1.e-20)
+                  IF( ze3f(ji,jj,jk) /= 0.e0 )   ze3f(ji,jj,jk) = 1.e0 / ze3f(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         CALL lbc_lnk( ze3f, 'F', 1. )
+      ENDIF
+
+      ! Local constant initialization
+      zfac12 = 1.e0 / 12.e0
+      
+      !                                                ! ===============
+      DO jk = 1, jpkm1                                 ! Horizontal slab
+         !                                             ! ===============
+         
+         ! Potential vorticity and horizontal fluxes
+         ! -----------------------------------------
+!!!bug   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) / fse3f(:,:,jk)
+         zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
+         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
+         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
+#if defined key_trddyn || defined key_trd_vor
+         zcor(:,:) = ff(:,:) * ze3f(:,:,jk)
+#endif
+
+         ! 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   
+               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  )
+               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
+               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
+         END DO
+         DO jj = 3, jpj
+            DO ji = fs_2, jpi   ! 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  )
+               ztse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
+               ztsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
+            END DO
+         END DO
+         
+         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  ) )
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
+               va(ji,jj,jk) = va(ji,jj,jk) + zva
+#   if defined key_trddyn || defined key_trd_vor
+               zcu = + zfac12 / e1u(ji,jj) * (  zcor(ji,jj  ) * zwy(ji  ,jj  ) + zcor(ji+1,jj) * zwy(ji+1,jj  )   &
+                  &                           + zcor(ji,jj  ) * zwy(ji  ,jj-1) + zcor(ji+1,jj) * zwy(ji+1,jj-1) )
+               zcv = - zfac12 / e2v(ji,jj) * (  zcor(ji,jj+1) * zwx(ji-1,jj+1) + zcor(ji,jj+1) * zwx(ji  ,jj+1)   &
+                  &                           + zcor(ji,jj  ) * zwx(ji-1,jj  ) + zcor(ji,jj  ) * zwx(ji  ,jj  ) )
+               utrd(ji,jj,jk,3) = zua - zcu
+               vtrd(ji,jj,jk,3) = zva - zcv
+               utrd(ji,jj,jk,4) = zcu
+               vtrd(ji,jj,jk,4) = zcv
+#   endif
+            END DO  
+         END DO  
+         !                                             ! ===============
+      END DO                                           !   End of slab
+      !                                                ! ===============
+
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
+         zua = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
+         zva = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
+         WRITE(numout,*) ' vor  een - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
+         u_ctl = zua   ;   v_ctl = zva
+      ENDIF
+
+   END SUBROUTINE dyn_vor_ene_ens
+
+
    SUBROUTINE dyn_vor_ctl
       !!---------------------------------------------------------------------
@@ -417,5 +555,5 @@
       INTEGER ::   ioptio = 0      ! temporary integer
 
-      NAMELIST/nam_dynvor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix
+      NAMELIST/nam_dynvor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
       !!----------------------------------------------------------------------
 
@@ -436,4 +574,5 @@
          WRITE(numout,*) '                 energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene
          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
       ENDIF
 
@@ -451,4 +590,9 @@
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '              vorticity term : mixed enstrophy/energy conserving scheme'
+         ioptio = ioptio + 1
+      ENDIF
+      IF( ln_dynvor_een ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '              vorticity term : energy and enstrophy conserving scheme'
          ioptio = ioptio + 1
       ENDIF
Index: /trunk/NEMO/OPA_SRC/step.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/step.F90	(revision 107)
+++ /trunk/NEMO/OPA_SRC/step.F90	(revision 108)
@@ -107,4 +107,5 @@
    USE diagap          ! hor. mean model-data gap         (dia_gap routine)
    USE trdmld          ! mixed-layer trends               (trd_mld routine)
+   USE trdvor          ! vorticity budget                 (trd_vor routine)
    USE diahdy          ! dynamic height                   (dia_hdy routine)
    USE ptr             ! poleward transports         (dia_ptr_init routine)
@@ -373,4 +374,5 @@
       IF( ln_dynvor_ene    )   CALL dyn_vor_energy   ( kstp )        ! energy conserving scheme
       IF( ln_dynvor_mix    )   CALL dyn_vor_mixed    ( kstp )        ! mixed energy/enstrophy conserving scheme
+      IF( ln_dynvor_een    )   CALL dyn_vor_ene_ens  ( kstp )        ! combined energy/enstrophy conserving scheme
 
       !                                                       ! lateral mixing 
@@ -439,4 +441,5 @@
          IF( lk_trdtra  )   CALL trd_tra( kstp )                 ! trends: active tracers
          IF( lk_trdmld  )   CALL trd_mld( kstp )                 ! trends: Mixed-layer 
+         IF( lk_trdvor  )   CALL trd_vor( kstp )                 ! trends: vorticity budget
          IF( lk_diaspr  )   CALL dia_spr( kstp )                 ! Surface pressure diagnostics
          IF( lk_diahth  )   CALL dia_hth( kstp )                 ! Thermocline depth (20 degres isotherm depth)
