Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 4622)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 5038)
@@ -18,4 +18,5 @@
    !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
    !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
+   !!            3.7  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
    !!-------------------------------------------------------------------------
   
@@ -34,14 +35,18 @@
    USE bdydyn          ! ocean open boundary conditions
    USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
+   USE trd_oce         ! trends: ocean variables
+   USE trddyn          ! trend manager: dynamics
+   USE trdken          ! trend manager: kinetic energy
+   !
    USE in_out_manager  ! I/O manager
+   USE iom             ! I/O manager library
    USE lbclnk          ! lateral boundary condition (or mpp link)
    USE lib_mpp         ! MPP library
    USE wrk_nemo        ! Memory Allocation
    USE prtctl          ! Print control
-
+   USE timing          ! Timing
 #if defined key_agrif
    USE agrif_opa_interp
 #endif
-   USE timing          ! Timing
 
    IMPLICIT NONE
@@ -79,5 +84,5 @@
       !!             at the local domain boundaries through lbc_lnk call,
       !!             at the one-way open boundaries (lk_bdy=T),
-      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
+      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
       !!
       !!              * Apply the time filter applied and swap of the dynamics
@@ -99,14 +104,14 @@
       REAL(wp) ::   z2dt         ! temporary scalar
 #endif
-      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
-      REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
+      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars
+      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_nxt')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
-      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva )
+      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
+      !
+      CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
+      IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve )
       !
       IF( kt == nit000 ) THEN
@@ -152,15 +157,16 @@
 
 # if defined key_dynspg_ts
+!!gm IF ( lk_dynspg_ts ) THEN ....
       ! Ensure below that barotropic velocities match time splitting estimate
       ! Compute actual transport and replace it with ts estimate at "after" time step
-      zua(:,:) = 0._wp
-      zva(:,:) = 0._wp
-      DO jk = 1, jpkm1
-         zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-         zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
+      zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
+      zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
+         zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
       END DO
       DO jk = 1, jpkm1
-         ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
-         va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
+         ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
+         va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
       END DO
 
@@ -175,4 +181,5 @@
          END DO  
       ENDIF
+!!gm ENDIF
 # endif
 
@@ -195,4 +202,23 @@
 # endif
 #endif
+
+      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
+         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step 
+         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
+         !
+         !                                  ! Kinetic energy and Conversion
+         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
+         !
+         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
+            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
+            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
+            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
+            CALL iom_put( "vtrd_tot", zva )
+         ENDIF
+         !
+         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
+         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
+         !                                  !  computation of the asselin filter trends)
+      ENDIF
 
       ! Time filter and swap of dynamics arrays
@@ -217,6 +243,6 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi    
-                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) )
-                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) )
+                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
+                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
                      !
                      ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
@@ -301,13 +327,13 @@
             ! Revert "before" velocities to time split estimate
             ! Doing it here also means that asselin filter contribution is removed  
-            zua(:,:) = 0._wp
-            zva(:,:) = 0._wp
-            DO jk = 1, jpkm1
-               zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
-               zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)    
+            zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
+            zve(:,:) = fse3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)    
+            DO jk = 2, jpkm1
+               zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
+               zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)    
             END DO
             DO jk = 1, jpkm1
-               ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
-               vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
+               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
+               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
             END DO
          ENDIF
@@ -327,6 +353,6 @@
             hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
          END DO
-         hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )
-         hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )
+         hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
+         hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
       ENDIF
       !
@@ -335,11 +361,6 @@
       !
       DO jk = 1, jpkm1
-#if defined key_vectopt_loop
-         DO jj = 1, 1         !Vector opt. => forced unrolling
-            DO ji = 1, jpij
-#else 
          DO jj = 1, jpj
             DO ji = 1, jpi
-#endif                  
                un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
                vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
@@ -358,9 +379,16 @@
       !
       !
+
+      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
+         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
+         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
+         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
+      ENDIF
+      !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
       ! 
-      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
-      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva )
+      CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
+      IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
@@ -370,3 +398,2 @@
    !!=========================================================================
 END MODULE dynnxt
-
