Index: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 7648)
+++ branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 7649)
@@ -35,9 +35,11 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE iom             ! I/O Manager for dyn_vrt_dia
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   div_cur    ! routine called by step.F90 and istate.F90
+   PUBLIC   div_cur     ! routine called by step.F90 and istate.F90
+   PUBLIC   dyn_vrt_dia ! routine called by various modules
 
    !! * Substitutions
@@ -338,4 +340,143 @@
    
 #endif
+
+
+   SUBROUTINE dyn_vrt_dia( utend, vtend, id_dia_vor_int, id_dia_vor_mn)
+
+      !!----------------------------------------------------------------------
+      !!
+      !! ** Purpose :  compute the integral and mean vorticity tendencies.
+      !!
+      !! ** Action : a) Calculate the vertical integrals of utend & of vtend
+      !!                (u_int & v_int)
+      !!             b) Calculate the vorticity tendencies for the vertical
+      !!                integrals.
+      !!             c) Calculate the vertical means, u_mn, v_mn from u_int
+      !!                and v_int by dividing by the depth
+      !!             d) Calculate the vorticity tendencies for the vertical
+      !!                means
+      !!             e) Call iom_put for the vertical integral vorticity
+      !!                tendencies (using id_dia_vor_int)
+      !!             f) Call iom_put for the vertical mean vorticity
+      !!                tendencies (using id_dia_vor_mn)
+      !!
+      !!----------------------------------------------------------------------
+      REAL :: utend(jpi,jpj,jpk) ! contribution to du/dt
+      REAL :: vtend(jpi,jpj,jpk) ! contribution to dv/dt
+      INTEGER :: id_dia_vor_int  ! identifier for the vertical integral vorticity diagnostic
+      INTEGER :: id_dia_vor_mn   ! identifier for the vertical mean vorticity diagnostic
+      !
+      !!----------------------------------------------------------------------
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ikbu, ikbv   ! dummy loop indices
+      !
+      REAL(wp), POINTER, DIMENSION(:,:) :: u_int   ! u vertical integral
+      REAL(wp), POINTER, DIMENSION(:,:) :: v_int   ! v vertical integral
+      REAL(wp), POINTER, DIMENSION(:,:) :: u_mn    ! u vertical means
+      REAL(wp), POINTER, DIMENSION(:,:) :: v_mn    ! u vertical means
+      REAL(wp), POINTER, DIMENSION(:,:) :: vor_int ! vort trend of vert integrals
+      REAL(wp), POINTER, DIMENSION(:,:) :: vor_mn  ! vort trend of vert means
+
+      CALL wrk_alloc(jpi, jpj, u_int)
+      CALL wrk_alloc(jpi, jpj, v_int)
+      CALL wrk_alloc(jpi, jpj, u_mn)
+      CALL wrk_alloc(jpi, jpj, v_mn)
+      CALL wrk_alloc(jpi, jpj, vor_int)
+      CALL wrk_alloc(jpi, jpj, vor_mn)
+
+      u_int(:,:) = 0.0_wp
+      v_int(:,:) = 0.0_wp
+
+      !
+      ! Calculate the vertical integrals of utend & of vtend
+      !
+      
+      DO jk = 1,jpk
+          DO jj = 1,jpj
+              DO ji = 1,jpi
+                  u_int(ji,jj) = u_int(ji,jj) + utend(ji,jj,jk)*fse3u(ji,jj,jk)
+                  v_int(ji,jj) = v_int(ji,jj) + vtend(ji,jj,jk)*fse3v(ji,jj,jk)
+              END DO
+          END DO
+      END DO
+
+      !
+      ! Calculate the vorticity tendencies for the vertical integrals.
+      ! 1/e1e2 * ((e2*d(vtend)/dx) - (e1*d(utend)/dy))
+      !
+
+      DO jj = 1,jpjm1
+          DO ji = 1,jpim1
+              vor_int(ji,jj) = ( v_int(ji+1,jj) * e2v(ji+1,jj)     &
+                  &            - v_int(ji,jj)   * e2v(ji,jj)       &
+                  &            + u_int(ji,jj)   * e1u(ji,jj)       &
+                  &            - u_int(ji,jj+1) * e1u(ji,jj+1) )   &
+                  &           / ( e1f(ji,jj)    * e2f(ji,jj) )
+          END DO
+      END DO
+
+
+      !
+      ! Calculate the vertical means, u_mn, v_mn from u_int & v_int by dividing
+      ! by the depth
+      ! mbku & mbkv - vertical index of the bottom last U- & W- ocean level
+      !
+
+      DO jj = 1, jpj
+          DO ji = 1, jpi
+              ikbu = mbku(ji,jj)
+              ikbv = mbkv(ji,jj)
+
+              IF (ikbu .ne. 0.0_wp) THEN      ! Don't divide by 0!
+                  u_mn(ji,jj) = u_int(ji,jj) / ikbu
+              ELSE
+                  u_mn(ji,jj) = 0.0_wp
+              END IF
+
+              IF (ikbv .ne. 0.0_wp) THEN      ! Don't divide by 0!
+                  v_mn(ji,jj) = v_int(ji,jj) / ikbv
+              ELSE
+                  v_mn(ji,jj) = 0.0_wp
+              END IF
+          END DO
+      END DO
+
+      !
+      ! Calculate the vorticity tendencies for the vertical means
+      ! 1/e1e2 * ((e2*d(v_mn)/dx) - (e1*d(u_mn)/dy))
+      !
+
+      DO jj = 1,jpjm1
+          DO ji = 1,jpim1
+              vor_mn(ji,jj) = ( v_mn(ji+1,jj) * e2v(ji+1,jj)     &
+                  &           - v_mn(ji,jj)   * e2v(ji,jj)       &
+                  &           + u_mn(ji,jj)   * e1u(ji,jj)       &
+                  &           - u_mn(ji,jj+1) * e1u(ji,jj+1) )   &
+                  &          / ( e1f(ji,jj)   * e2f(ji,jj) )
+          END DO
+      END DO
+
+
+      ! Call iom_put for the vertical integral vorticity tendencies
+      IF (id_dia_vor_int == 1) THEN
+          CALL iom_put( "dia_vor_int", vor_int(:,:))
+      ENDIF
+
+      ! Call iom_put for the vertical mean vorticity tendencies
+      IF (id_dia_vor_int == 1) THEN
+          CALL iom_put( "dia_vor_mn", vor_mn(:,:))
+      ENDIF
+
+      CALL wrk_dealloc(jpi, jpj, u_int)
+      CALL wrk_dealloc(jpi, jpj, v_int)
+      CALL wrk_dealloc(jpi, jpj, u_mn)
+      CALL wrk_dealloc(jpi, jpj, v_mn)
+      CALL wrk_dealloc(jpi, jpj, vor_int)
+      CALL wrk_dealloc(jpi, jpj, vor_mn)
+
+   END SUBROUTINE dyn_vrt_dia
+
+
    !!======================================================================
 END MODULE divcur
Index: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 7648)
+++ branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 7649)
@@ -36,4 +36,5 @@
    USE trd_oce         ! trends: ocean variables
    USE trddyn          ! trend manager: dynamics
+   USE divcur          ! for dyn_vrt_dia
    !
    USE in_out_manager  ! I/O manager
@@ -378,4 +379,6 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
+      INTEGER  ::   id_dia_vrt_hpg_int  = 1    ! TODO remove once flags set properly
+      INTEGER  ::   id_dia_vrt_hpg_mean = 1    ! TODO remove once flags set properly
       INTEGER  ::   ji, jj, jk                 ! dummy loop indices
       REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars
@@ -401,9 +404,4 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ! hydrostatic pressure gradient along s-surfaces
-            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   &
-               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
-            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   &
-               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
             ! s-coordinate pressure gradient correction
             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
@@ -411,7 +409,12 @@
             zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
                &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
+            ! hydrostatic pressure gradient along s-surfaces
+            zhpi(ji,jj,1) = zuap + zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   &
+               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
+            zhpj(ji,jj,1) = zvap + zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   &
+               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
             ! add to the general momentum trend
-            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
-            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
+            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
+            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
          END DO
       END DO
@@ -421,11 +424,4 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! hydrostatic pressure gradient along s-surfaces
-               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
-                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
-                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
-               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
-                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
-                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
                ! s-coordinate pressure gradient correction
                zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
@@ -433,10 +429,22 @@
                zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
                   &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
+               ! hydrostatic pressure gradient along s-surfaces with correction
+               zhpi(ji,jj,jk) = zuap + zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
+                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
+                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
+               zhpj(ji,jj,jk) = zvap + zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
+                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
+                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
                ! add to the general momentum trend
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
-               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
             END DO
          END DO
       END DO
+      !
+      ! calculate dia_vor_int & dia_vor_mn if required
+      IF ( ( id_dia_vrt_hpg_int == 1 ) .or. ( id_dia_vrt_hpg_mean == 1 ) ) THEN
+          CALL dyn_vrt_dia(zhpi, zhpj, id_dia_vrt_hpg_int, id_dia_vrt_hpg_mean)
+      END IF
       !
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
Index: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 7648)
+++ branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 7649)
@@ -17,4 +17,5 @@
    USE trd_oce         ! trends: ocean variables
    USE trddyn          ! trend manager: dynamics
+   USE divcur          ! for dyn_vrt_dia
    !
    USE in_out_manager  ! I/O manager
@@ -71,10 +72,12 @@
       !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
-      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme 
+      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kscheme    ! =0/1   type of KEG scheme 
       !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zu, zv       ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
+      INTEGER  ::   id_dia_vrt_keg_int  = 1 ! TODO remove once flags set properly
+      INTEGER  ::   id_dia_vrt_keg_mean = 1 ! TODO remove once flags set properly
+      INTEGER  ::   ji, jj, jk              ! dummy loop indices
+      REAL(wp) ::   zu, zv                  ! temporary scalars
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke, zhkei, zhkej
       REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
       !!----------------------------------------------------------------------
@@ -82,5 +85,5 @@
       IF( nn_timing == 1 )   CALL timing_start('dyn_keg')
       !
-      CALL wrk_alloc( jpi,jpj,jpk,   zhke )
+      CALL wrk_alloc( jpi,jpj,jpk,   zhke, zhkei, zhkej )
       !
       IF( kt == nit000 ) THEN
@@ -137,9 +140,15 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
-               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
+               zhkei(ji,jj,jk) = - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
+               zhkej(ji,jj,jk) = - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zhkei(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + zhkej(ji,jj,jk)
             END DO 
          END DO
       END DO
+      !
+      IF ( ( id_dia_vrt_keg_int == 1 ) .or. ( id_dia_vrt_keg_mean == 1 ) ) THEN
+          CALL dyn_vrt_dia(zhkei, zhkej, id_dia_vrt_keg_int, id_dia_vrt_keg_mean)
+      END IF
       !
       IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
@@ -153,5 +162,5 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zhke )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zhke, zhkei, zhkej )
       !
       IF( nn_timing == 1 )   CALL timing_stop('dyn_keg')
Index: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
===================================================================
--- branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 7648)
+++ branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 7649)
@@ -19,4 +19,5 @@
    USE dom_oce         ! ocean space and time domain
    USE ldfdyn_oce      ! ocean dynamics: lateral physics
+   USE divcur          ! for dyn_vrt_dia
    !
    USE in_out_manager  ! I/O manager
@@ -75,8 +76,10 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
+      INTEGER  :: id_dia_vrt_ldf_int  = 1      ! TODO remove once flags set properly
+      INTEGER  :: id_dia_vrt_ldf_mean = 1      ! TODO remove once flags set properly
       INTEGER  ::   ji, jj, jk                  ! dummy loop indices
-      REAL(wp) ::   zua, zva, zbt, ze2u, ze2v   ! temporary scalar
+      REAL(wp) ::   zbt, ze2u, ze2v             ! temporary scalar
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zcu, zcv
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zuf, zut, zlu, zlv
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zuf, zut, zlu, zlv, zua, zva
       !!----------------------------------------------------------------------
       !
@@ -84,5 +87,5 @@
       !
       CALL wrk_alloc( jpi, jpj,      zcu, zcv           )
-      CALL wrk_alloc( jpi, jpj, jpk, zuf, zut, zlu, zlv ) 
+      CALL wrk_alloc( jpi, jpj, jpk, zuf, zut, zlu, zlv, zua, zva ) 
       !
       IF( kt == nit000 .AND. lwp ) THEN
@@ -191,12 +194,13 @@
                ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
                ! horizontal biharmonic diffusive trends
-               zua = - ( zuf(ji  ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u   &
-                  &  + ( zut(ji+1,jj,jk) - zut(ji,jj  ,jk) ) / e1u(ji,jj)
-
-               zva = + ( zuf(ji,jj  ,jk) - zuf(ji-1,jj,jk) ) / ze2v   &
-                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj)
+               zua(ji,jj,jk) = - ( zuf(ji  ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u   &
+                  &  + ( zut(ji+1,jj,jk) - zut(ji,jj  ,jk) ) / e1u(ji,jj)       &
+                  &  * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag ))
+               zva(ji,jj,jk) = + ( zuf(ji,jj  ,jk) - zuf(ji-1,jj,jk) ) / ze2v   &
+                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj)       &
+                  &  * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag ))
                ! add it to the general momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zua * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag ))
-               va(ji,jj,jk) = va(ji,jj,jk) + zva * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag ))
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zua(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + zva(ji,jj,jk)
             END DO
          END DO
@@ -205,6 +209,10 @@
       END DO                                           !   End of slab
       !                                                ! ===============
+      IF ( ( id_dia_vrt_ldf_int == 1 ) .or. ( id_dia_vrt_ldf_mean == 1 ) ) THEN
+          CALL dyn_vrt_dia(zua, zva, id_dia_vrt_ldf_int, id_dia_vrt_ldf_mean)
+      END IF
+      !
       CALL wrk_dealloc( jpi, jpj,      zcu, zcv           )
-      CALL wrk_dealloc( jpi, jpj, jpk, zuf, zut, zlu, zlv ) 
+      CALL wrk_dealloc( jpi, jpj, jpk, zuf, zut, zlu, zlv, zua, zva ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_bilap')
Index: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 7648)
+++ branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 7649)
@@ -38,4 +38,5 @@
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
+   USE divcur         ! For dyn_vrt_dia
 
 
@@ -560,11 +561,14 @@
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend
       !!
+      INTEGER  :: id_dia_vrt_vor_int = 1                          ! TODO remove once flags set properly
+      INTEGER  :: id_dia_vrt_vor_mn  = 1                          ! TODO remove once flags set properly
       INTEGER  ::   ji, jj, jk                                    ! dummy loop indices
       INTEGER  ::   ierr                                          ! local integer
-      REAL(wp) ::   zfac12, zua, zva                              ! local scalars
+      REAL(wp) ::   zfac12                                        ! 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
+      REAL(wp), POINTER    , DIMENSION(:,:,:)         :: zua, zva
 #if defined key_vvl
       REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T)
@@ -578,4 +582,5 @@
       CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        ) 
       CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
+      CALL wrk_alloc( jpi, jpj, jpk, zua, zva               ) 
 #if defined key_vvl
       CALL wrk_alloc( jpi, jpj, jpk, ze3f                   )
@@ -687,10 +692,14 @@
          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  ) )
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
+               zua(ji,jj,jk) = + 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(ji,jj,jk) = - 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  ) )
+               pua(ji,jj,jk) = pua(ji,jj,jk) + zua(ji,jj,jk)
+               pva(ji,jj,jk) = pva(ji,jj,jk) + zva(ji,jj,jk)
             END DO  
          END DO  
@@ -698,6 +707,11 @@
       END DO                                           !   End of slab
       !                                                ! ===============
+      IF ( ( id_dia_vrt_vor_int == 1 ) .or. ( id_dia_vrt_vor_mn == 1 ) ) THEN
+          CALL dyn_vrt_dia(zua, zva, id_dia_vrt_vor_int, id_dia_vrt_vor_mn)
+      END IF
+      !
       CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        ) 
       CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse ) 
+      CALL wrk_dealloc( jpi, jpj, jpk, zua, zva               ) 
 #if defined key_vvl
       CALL wrk_dealloc( jpi, jpj, jpk, ze3f                   )
Index: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 7648)
+++ branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 7649)
@@ -18,4 +18,5 @@
    USE trd_oce        ! trends: ocean variables
    USE trddyn         ! trend manager: dynamics
+   USE divcur         ! for dyn_vrt_dia
    !
    USE in_out_manager ! I/O manager
@@ -59,9 +60,11 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
       !
-      INTEGER  ::   ji, jj, jk      ! dummy loop indices
-      REAL(wp) ::   zua, zva        ! temporary scalars
+      INTEGER  ::   id_dia_vrt_zad_int  = 1  ! TODO remove once flags set properly
+      INTEGER  ::   id_dia_vrt_zad_mean = 1  ! TODO remove once flags set properly
+      INTEGER  ::   ji, jj, jk               ! dummy loop indices
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwuw , zwvw
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  zww
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zua, zva
       !!----------------------------------------------------------------------
       !
@@ -69,5 +72,5 @@
       !
       CALL wrk_alloc( jpi,jpj, zww ) 
-      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw ) 
+      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zua, zva ) 
       !
       IF( kt == nit000 ) THEN
@@ -121,9 +124,11 @@
             DO ji = fs_2, fs_jpim1       ! vector opt.
                !                         ! vertical momentum advective trends
-               zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
-               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
+               zua(ji,jj,jk) = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) /            &
+                   &                ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
+               zva(ji,jj,jk) = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) /            &
+                   &                ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
                !                         ! add the trends to the general momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
-               va(ji,jj,jk) = va(ji,jj,jk) + zva
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zua(ji,jj,jk)
+               va(ji,jj,jk) = va(ji,jj,jk) + zva(ji,jj,jk)
             END DO  
          END DO  
@@ -140,6 +145,10 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
+      IF ( ( id_dia_vrt_zad_int == 1 ) .or. ( id_dia_vrt_zad_mean == 1 ) ) THEN
+          CALL dyn_vrt_dia(zua, zva, id_dia_vrt_zad_int, id_dia_vrt_zad_mean)
+      END IF
+      !
       CALL wrk_dealloc( jpi,jpj, zww ) 
-      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw ) 
+      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zua, zva ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_zad')
