Index: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor.F90
===================================================================
--- branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor.F90	(revision 4763)
+++ branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor.F90	(revision 5208)
@@ -4,10 +4,9 @@
    !! Ocean diagnostics:  momentum trends
    !!=====================================================================
-   !! History :  1.0  !  04-2006  (L. Brunier, A-M. Treguier) Original code 
-   !!            2.0  !  04-2008  (C. Talandier) New trends organization
+   !! History :  1.0  !  2006-01  (L. Brunier, A-M. Treguier) Original code 
+   !!            2.0  !  2008-04  (C. Talandier) New trends organization
+   !!            3.5  !  2012-02  (G. Madec) regroup beta.V computation with pvo trend
    !!----------------------------------------------------------------------
-#if defined key_trdvor   ||   defined key_esopa
-   !!----------------------------------------------------------------------
-   !!   'key_trdvor'   : momentum trend diagnostics
+
    !!----------------------------------------------------------------------
    !!   trd_vor      : momentum trends averaged over the depth
@@ -17,16 +16,16 @@
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables
-   USE trdmod_oce      ! ocean variables trends
+   USE trd_oce         ! trends: ocean variables
    USE zdf_oce         ! ocean vertical physics
-   USE in_out_manager  ! I/O manager
+   USE sbc_oce         ! surface boundary condition: ocean
    USE phycst          ! Define parameters for the routines
    USE ldfdyn_oce      ! ocean active tracers: lateral physics
    USE dianam          ! build the name of file (routine)
    USE zdfmxl          ! mixed layer depth
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager  ! I/O manager
    USE ioipsl          ! NetCDF library
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! MPP library
    USE wrk_nemo        ! Memory allocation
-
 
    IMPLICIT NONE
@@ -37,6 +36,5 @@
    END INTERFACE
 
-   PUBLIC   trd_vor        ! routine called by step.F90
-   PUBLIC   trd_vor_zint   ! routine called by dynamics routines
+   PUBLIC   trd_vor        ! routine called by trddyn.F90
    PUBLIC   trd_vor_init   ! routine called by opa.F90
    PUBLIC   trd_vor_alloc  ! routine called by nemogcm.F90
@@ -80,4 +78,52 @@
       IF( trd_vor_alloc /= 0 )   CALL ctl_warn('trd_vor_alloc: failed to allocate arrays')
    END FUNCTION trd_vor_alloc
+
+
+   SUBROUTINE trd_vor( putrd, pvtrd, ktrd, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_vor  ***
+      !! 
+      !! ** Purpose :  computation of cumulated trends over analysis period
+      !!               and make outputs (NetCDF or DIMG format)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends 
+      INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
+      INTEGER                   , INTENT(in   ) ::   kt             ! time step
+      !
+      INTEGER ::   ji, jj   ! dummy loop indices
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztswu, ztswv    ! 2D workspace 
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi, jpj, ztswu, ztswv )
+
+      SELECT CASE( ktrd ) 
+      CASE( jpdyn_hpg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_prg )   ! Hydrostatique Pressure Gradient 
+      CASE( jpdyn_keg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_keg )   ! KE Gradient 
+      CASE( jpdyn_rvo )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_rvo )   ! Relative Vorticity 
+      CASE( jpdyn_pvo )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_pvo )   ! Planetary Vorticity Term 
+      CASE( jpdyn_ldf )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_ldf )   ! Horizontal Diffusion 
+      CASE( jpdyn_zad )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_zad )   ! Vertical Advection 
+      CASE( jpdyn_spg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_spg )   ! Surface Pressure Grad. 
+      CASE( jpdyn_zdf )                                                      ! Vertical Diffusion 
+         ztswu(:,:) = 0.e0   ;   ztswv(:,:) = 0.e0
+         DO jj = 2, jpjm1                                                             ! wind stress trends
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 )
+               ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 )
+            END DO
+         END DO
+         !
+         CALL trd_vor_zint( putrd, pvtrd, jpvor_zdf )                             ! zdf trend including surf./bot. stresses 
+         CALL trd_vor_zint( ztswu, ztswv, jpvor_swf )                             ! surface wind stress 
+      CASE( jpdyn_bfr )
+         CALL trd_vor_zint( putrd, pvtrd, jpvor_bfr )                             ! Bottom stress
+         !
+      CASE( jpdyn_atf )       ! last trends: perform the output of 2D vorticity trends
+         CALL trd_vor_iom( kt )
+      END SELECT
+      !
+      CALL wrk_dealloc( jpi, jpj, ztswu, ztswv )
+      !
+   END SUBROUTINE trd_vor
 
 
@@ -109,5 +155,4 @@
       !!      trends output in netCDF format using ioipsl
       !!----------------------------------------------------------------------
-      !
       INTEGER                     , INTENT(in   ) ::   ktrd       ! ocean trend index
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   putrdvor   ! u vorticity trend 
@@ -131,7 +176,7 @@
       !  =====================================
 
-      SELECT CASE (ktrd) 
-      !
-      CASE (jpvor_bfr)        ! bottom friction
+      SELECT CASE( ktrd ) 
+      !
+      CASE( jpvor_bfr )        ! bottom friction
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1 
@@ -143,5 +188,5 @@
          END DO
          !
-      CASE (jpvor_swf)        ! wind stress
+      CASE( jpvor_swf )        ! wind stress
          zudpvor(:,:) = putrdvor(:,:) * fse3u(:,:,1) * e1u(:,:) * umask(:,:,1)
          zvdpvor(:,:) = pvtrdvor(:,:) * fse3v(:,:,1) * e2v(:,:) * vmask(:,:,1)
@@ -154,6 +199,6 @@
    
       ! Curl
-      DO ji=1,jpim1
-         DO jj=1,jpjm1
+      DO ji = 1, jpim1
+         DO jj = 1, jpjm1
             vortrd(ji,jj,ktrd) = (    zvdpvor(ji+1,jj) - zvdpvor(ji,jj)       &
                  &                - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) )   ) / ( e1f(ji,jj) * e2f(ji,jj) )
@@ -229,15 +274,23 @@
       END DO
 
-      ! Save Beta.V term to avoid average before Curl
-      ! Beta.V : intergration, no average
-      IF( ktrd == jpvor_bev ) THEN 
+      ! Planetary vorticity: 2nd computation (Beta.V term) store the vertical sum
+      ! as Beta.V term need intergration, not average
+      IF( ktrd == jpvor_pvo ) THEN 
          zubet(:,:) = zudpvor(:,:)
          zvbet(:,:) = zvdpvor(:,:)
-      ENDIF
-
-      ! Average except for Beta.V
+         DO ji = 1, jpim1
+            DO jj = 1, jpjm1
+               vortrd(ji,jj,jpvor_bev) = (    zvbet(ji+1,jj) - zvbet(ji,jj)     &
+                  &                       - ( zubet(ji,jj+1) - zubet(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) )
+            END DO
+         END DO
+         ! Average of the Curl and Surface mask
+         vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * hur(:,:) * fmask(:,:,1)
+      ENDIF
+      !
+      ! Average 
       zudpvor(:,:) = zudpvor(:,:) * hur(:,:)
       zvdpvor(:,:) = zvdpvor(:,:) * hvr(:,:)
-   
+      !
       ! Curl
       DO ji=1,jpim1
@@ -247,24 +300,6 @@
          END DO
       END DO
-
       ! Surface mask
       vortrd(:,:,ktrd) = vortrd(:,:,ktrd) * fmask(:,:,1)
-
-      ! Special treatement for the Beta.V term
-      ! Compute the Curl of the Beta.V term which is not averaged
-      IF( ktrd == jpvor_bev ) THEN
-         DO ji=1,jpim1
-            DO jj=1,jpjm1
-               vortrd(ji,jj,jpvor_bev) = (    zvbet(ji+1,jj) - zvbet(ji,jj)     &
-                  &                       - ( zubet(ji,jj+1) - zubet(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) )
-            END DO
-         END DO
-
-         ! Average on the Curl
-         vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * hur(:,:)
-
-         ! Surface mask
-         vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * fmask(:,:,1)
-      ENDIF
    
       IF( ndebug /= 0 ) THEN
@@ -278,5 +313,5 @@
 
 
-   SUBROUTINE trd_vor( kt )
+   SUBROUTINE trd_vor_iom( kt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trd_vor  ***
@@ -285,6 +320,5 @@
       !!               and make outputs (NetCDF or DIMG format)
       !!----------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER                   , INTENT(in   ) ::   kt             ! time step
       !
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
@@ -305,9 +339,4 @@
 
       IF( kt > nit000 )   vor_avrb(:,:) = vor_avr(:,:)
-
-      IF( ndebug /= 0 ) THEN
-          WRITE(numout,*) ' debuging trd_vor: I.1 done '
-          CALL FLUSH(numout)
-      ENDIF
 
       ! I.2 vertically integrated vorticity
@@ -330,6 +359,6 @@
 
       ! Curl
-      DO ji=1,jpim1
-         DO jj=1,jpjm1
+      DO ji = 1, jpim1
+         DO jj = 1, jpjm1
             vor_avr(ji,jj) = (  ( zvn(ji+1,jj) - zvn(ji,jj) )    &
                &              - ( zun(ji,jj+1) - zun(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,1)
@@ -337,9 +366,4 @@
       END DO
       
-      IF( ndebug /= 0 ) THEN
-         WRITE(numout,*) ' debuging trd_vor: I.2 done'
-         CALL FLUSH(numout)
-      ENDIF
-
       !  =================================
       !   II. Cumulated trends
@@ -351,9 +375,4 @@
          vor_avrbb(:,:) = vor_avrb(:,:)
          vor_avrbn(:,:) = vor_avr (:,:)
-      ENDIF
-
-      IF( ndebug /= 0 ) THEN
-         WRITE(numout,*) ' debuging trd_vor: I1.1 done'
-         CALL FLUSH(numout)
       ENDIF
 
@@ -371,9 +390,4 @@
       ENDIF
 
-      IF( ndebug /= 0 ) THEN
-         WRITE(numout,*) ' debuging trd_vor: II.2 done'
-         CALL FLUSH(numout)
-      ENDIF
-
       !  =============================================
       !   III. Output in netCDF + residual computation
@@ -391,9 +405,4 @@
          vor_avrtot(:,:) = (  vor_avr(:,:) - vor_avrbn(:,:) + vor_avrb(:,:) - vor_avrbb(:,:) ) * zmean
 
-         IF( ndebug /= 0 ) THEN
-             WRITE(numout,*) ' zmean = ',zmean
-             WRITE(numout,*) ' debuging trd_vor: III.1 done'
-             CALL FLUSH(numout)
-         ENDIF
 
          ! III.2 compute residual
@@ -406,8 +415,4 @@
          CALL lbc_lnk( vor_avrres, 'F', 1. )
 
-         IF( ndebug /= 0 ) THEN
-            WRITE(numout,*) ' debuging trd_vor: III.2 done'
-            CALL FLUSH(numout)
-         ENDIF
 
          ! III.3 time evolution array swap
@@ -415,9 +420,4 @@
          vor_avrbb(:,:) = vor_avrb(:,:)
          vor_avrbn(:,:) = vor_avr (:,:)
-
-         IF( ndebug /= 0 ) THEN
-            WRITE(numout,*) ' debuging trd_vor: III.3 done'
-            CALL FLUSH(numout)
-         ENDIF
          !
          nmoydpvor = 0
@@ -463,5 +463,5 @@
       CALL wrk_dealloc( jpi, jpj, zun, zvn )                                   
       !
-   END SUBROUTINE trd_vor
+   END SUBROUTINE trd_vor_iom
 
 
@@ -587,29 +587,4 @@
    END SUBROUTINE trd_vor_init
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                       Empty module
-   !!----------------------------------------------------------------------
-   INTERFACE trd_vor_zint
-      MODULE PROCEDURE trd_vor_zint_2d, trd_vor_zint_3d
-   END INTERFACE
-CONTAINS
-   SUBROUTINE trd_vor( kt )        ! Empty routine
-      WRITE(*,*) 'trd_vor: You should not have seen this print! error?', kt
-   END SUBROUTINE trd_vor
-   SUBROUTINE trd_vor_zint_2d( putrdvor, pvtrdvor, ktrd )
-      REAL, DIMENSION(:,:), INTENT( inout ) ::   putrdvor, pvtrdvor
-      INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index
-      WRITE(*,*) 'trd_vor_zint_2d: You should not have seen this print! error?', putrdvor(1,1), pvtrdvor(1,1), ktrd
-   END SUBROUTINE trd_vor_zint_2d
-   SUBROUTINE trd_vor_zint_3d( putrdvor, pvtrdvor, ktrd )
-      REAL, DIMENSION(:,:,:), INTENT( inout ) ::   putrdvor, pvtrdvor
-      INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index
-      WRITE(*,*) 'trd_vor_zint_3d: You should not have seen this print! error?', putrdvor(1,1,1), pvtrdvor(1,1,1), ktrd
-   END SUBROUTINE trd_vor_zint_3d
-   SUBROUTINE trd_vor_init              ! Empty routine
-      WRITE(*,*) 'trd_vor_init: You should not have seen this print! error?'
-   END SUBROUTINE trd_vor_init
-#endif
    !!======================================================================
 END MODULE trdvor
