Index: trunk/NEMO/OPA_SRC/TRD/trddyn.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trddyn.F90	(revision 129)
+++ 	(revision )
@@ -1,456 +1,0 @@
-MODULE trddyn
-   !!======================================================================
-   !!                       ***  MODULE  trddyn  ***
-   !! Ocean trends :  ocean momentum trends
-   !!=====================================================================
-#if   defined key_trddyn   ||   defined key_esopa
-   !!----------------------------------------------------------------------
-   !!   'key_trddyn'                                   momentum trend diag.
-   !!----------------------------------------------------------------------
-   !!   trd_dyn      : verify the basin averaged properties for momentum
-   !!   trd_dyn_init :
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE oce             ! ocean dynamics and tracers 
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE trdtra_oce      ! ocean active tracer trend 
-   USE trddyn_oce      ! ocean dynamics trend
-   USE zdf_oce         ! ocean vertical physics
-   USE ldftra_oce      ! ocean active tracers: lateral physics
-   USE ldfdyn_oce      ! ocean dynamics: lateral physics
-   USE sol_oce         ! solver variables
-   USE obc_oce         ! ocean lateral open boundary condition
-   USE eosbn2          ! equation of state (eos routine)
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distributed memory computing library
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! * Routine accessibility
-   PUBLIC trd_dyn      ! called by step.F90
-   PUBLIC trd_dyn_init ! called by step.F90
-
-   !! * Shared module vaiables
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .TRUE.    !: momentum trend flag
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LODYC-IPSL  (2003)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE trd_dyn( kt )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_dyn  ***
-      !!
-      !! ** Purpose :   verify the basin averaged properties of the momentum 
-      !!      equation at every time step frequency ntrd.
-      !!      momentum equation:
-      !!       * non linear momentum trend (relative vorticity trend + hori-
-      !!      zontal kinetic energy gradient trend + vertical advection
-      !!      trend) : 
-      !!      1. conserve the momentum
-      !!      2. conserve the potential relative vorticity (rotn/e3f)
-      !!         (except for the vertical advection)
-      !!      3. conserve the potential enstrophy (except for the
-      !!         vertical advection trend) IF no cpp key activated
-      !!        or #defined key_vorcombined
-      !!      4. conserve the horizontal kinetic energy if #defined
-      !!        key_vorenergy
-      !!
-      !! ** Method :
-      !!      damping trend
-      !!      horizontal pressure gradient trend
-      !!      horizontal kinetic energy gradient trend
-      !!      relative vorticity term trend
-      !!      coriolis trend
-      !!
-      !! History :
-      !!        !  91-12  (G. Madec)  Original code
-      !!        !  92-06  (M. Imbard)  add time step frequency
-      !!        !  96-01  (G. Madec)  Statement function for e3
-      !!                              suppression of common work arrays
-      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
-      !!---------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk, jn        ! dummy loop indices
-
-      REAL(wp) ::   ze1e2w,zbe1ru,zbe2rv,zbtr,zth,ztz,zpeke,zcof
-      REAL(wp) ::   zmsku, zmskv
-      REAL(wp) ::   zumo(11),zvmo(11),zhke(10)
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         zkepe, zkx, zky, zkz            ! workspace
-      
-      NAMELIST/namtrd/ ntrd, nctls
-      !!----------------------------------------------------------------------
-
-
-      ! 1. forcing trend and mask trends
-      ! --------------------------------
-      
-      IF( MOD(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
-
-         ! Mask surface forcing and bottom friction fluxes 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
-               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
-               tautrd(ji,jj,1) = tautrd(ji,jj,1) * zmsku
-               tautrd(ji,jj,2) = tautrd(ji,jj,2) * zmskv
-               tautrd(ji,jj,3) = tautrd(ji,jj,3) * zmsku
-               tautrd(ji,jj,4) = tautrd(ji,jj,4) * zmskv
-            END DO
-         END DO
-         tautrd( : ,jpj,1:4) = 0.e0
-         tautrd(jpi, : ,1:4) = 0.e0
-         
-         ! Mask the trends
-         DO jn = 1,9
-            DO jk = 1,jpk
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
-                     zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
-                     utrd(ji,jj,jk,jn) = utrd(ji,jj,jk,jn) * zmsku
-                     vtrd(ji,jj,jk,jn) = vtrd(ji,jj,jk,jn) * zmskv
-                  END DO
-               END DO
-            END DO
-         END DO
-         utrd( : ,jpj,:,1:9) = 0.e0      ;      vtrd( : ,jpj,:,1:9) = 0.e0
-         utrd(jpi, : ,:,1:9) = 0.e0      ;      vtrd(jpi, : ,:,1:9) = 0.e0
-         
-         ! Conversion potential energy - kinetic energy
-         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
-
-         zkx(:,:,:) = 0.e0
-         zky(:,:,:) = 0.e0
-         zkz(:,:,:) = 0.e0
-      
-         CALL eos( tn, sn, rhd, rhop )       ! now potential and in situ densities
-
-         ! Density flux at w-point
-         DO jk = 2, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ze1e2w = 0.5 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) * tmask_i(ji,jj)
-                  zkz(ji,jj,jk) = ze1e2w / rau0 * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) )
-               END DO
-            END DO
-         END DO
-         zkz  (:,:, 1 ) = 0.e0
-         
-         ! Density flux at u and v-points
-         DO jk = 1, jpk
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  zcof   = 0.5 / rau0
-                  zbe1ru = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
-                  zbe2rv = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
-                  zkx(ji,jj,jk) = zbe1ru * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
-                  zky(ji,jj,jk) = zbe2rv * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
-               END DO
-            END DO
-         END DO
-         
-         ! Density flux divergence at t-point
-         DO jk = 1, jpkm1
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
-                  ztz = - zbtr * (    zkz(ji,jj,jk) - zkz(ji,jj,jk+1) )
-                  zth = - zbtr * (  ( zkx(ji,jj,jk) - zkx(ji-1,jj,jk) )   &
-                                  + ( zky(ji,jj,jk) - zky(ji,jj-1,jk) )  )
-                  zkepe(ji,jj,jk) = (zth + ztz) * tmask(ji,jj,jk) * tmask_i(ji,jj)
-               END DO
-            END DO
-         END DO
-         zkepe( : , : ,jpk) = 0.e0
-         zkepe( : ,jpj, : ) = 0.e0
-         zkepe(jpi, : , : ) = 0.e0
-
-
-         ! 2. Basin averaged momentum trend
-         ! --------------------------------
-         
-         DO jn = 1,9
-            zumo(jn) = 0.
-            zvmo(jn) = 0.
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zumo(jn) = zumo(jn) + utrd(ji,jj,jk,jn) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-                     zvmo(jn) = zvmo(jn) + vtrd(ji,jj,jk,jn) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         END DO
-
-         zumo(10) = 0.
-         zvmo(10) = 0.
-         zumo(11) = 0.
-         zvmo(11) = 0.
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zumo(10) = zumo(10) + tautrd(ji,jj,1) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
-               zvmo(10) = zvmo(10) + tautrd(ji,jj,2) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
-               zumo(11) = zumo(11) + tautrd(ji,jj,3)
-               zvmo(11) = zvmo(11) + tautrd(ji,jj,4)
-            END DO
-         END DO
-         
-         
-         ! 3. Basin averaged kinetic energy trend
-         ! --------------------------------------
-         
-         ! Field before, because after the array swap
-         
-         DO jn = 1,9
-            zhke(jn) = 0.e0
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zhke(jn) = zhke(jn)   &
-                     &   + ub(ji,jj,jk) * utrd(ji,jj,jk,jn) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)   &
-                     &   + vb(ji,jj,jk) * vtrd(ji,jj,jk,jn) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         END DO
-         
-         zhke(10) = 0.e0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zhke(10) = zhke(10)   &
-                  &     + ub(ji,jj,1) * tautrd(ji,jj,1) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)   &
-                  &     + vb(ji,jj,1) * tautrd(ji,jj,2) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
-            END DO
-         END DO
-         
-         zpeke = 0.e0
-         DO jk = 1,jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zpeke    = zpeke + zkepe(ji,jj,jk) * grav * fsdept(ji,jj,jk)   &
-                     &                     * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         
-         IF( lk_mpp ) THEN
-            CALL mpp_sum( zpeke )
-            CALL mpp_sum( zumo , 11 )
-            CALL mpp_sum( zvmo , 11 )
-            CALL mpp_sum( zhke , 10 )
-         ENDIF
-
-
-         ! 4. Print
-         ! --------
-
-         IF(lwp) THEN
-            WRITE (numout,*)
-            WRITE (numout,*)
-            WRITE (numout,9400) kt
-            WRITE (numout,9401) zumo( 1) / tvolu, zvmo( 1) / tvolv
-            WRITE (numout,9402) zumo( 2) / tvolu, zvmo( 2) / tvolv
-            WRITE (numout,9403) zumo( 3) / tvolu, zvmo( 3) / tvolv
-            WRITE (numout,9404) zumo( 4) / tvolu, zvmo( 4) / tvolv
-            WRITE (numout,9405) zumo( 5) / tvolu, zvmo( 5) / tvolv
-            WRITE (numout,9406) zumo( 6) / tvolu, zvmo( 6) / tvolv
-            WRITE (numout,9407) zumo( 7) / tvolu, zvmo( 7) / tvolv
-            WRITE (numout,9408) zumo( 8) / tvolu, zvmo( 8) / tvolv
-            WRITE (numout,9409) zumo(10) / tvolu, zvmo(10) / tvolv
-            WRITE (numout,9410) zumo( 9) / tvolu, zvmo( 9) / tvolv
-            WRITE (numout,9411) zumo(11) / tvolu, zvmo(11) / tvolv
-            WRITE (numout,9412)
-            WRITE (numout,9413)                                                     &
-            &     (  zumo(1) + zumo(2) + zumo(3) + zumo(4) + zumo(5) + zumo(6)      &
-            &      + zumo(7) + zumo(8) + zumo(9) + zumo(10) + zumo(11) ) / tvolu,   &
-            &     (  zvmo(1) + zvmo(2) + zvmo(3) + zvmo(4) + zvmo(5) + zvmo(6)      &
-            &      + zvmo(7) + zvmo(8) + zvmo(9) + zvmo(10) + zvmo(11) ) / tvolv
-         ENDIF
-
- 9400    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
- 9401    FORMAT(' pressure gradient          u= ', e20.13, '    v= ', e20.13)
- 9402    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
- 9403    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
- 9404    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
- 9405    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
- 9406    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
- 9407    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
- 9408    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
- 9409    FORMAT(' forcing term               u= ', e20.13, '    v= ', e20.13)
- 9410    FORMAT(' dampimg term               u= ', e20.13, '    v= ', e20.13)
- 9411    FORMAT(' bottom flux                u= ', e20.13, '    v= ', e20.13)
- 9412    FORMAT(' -----------------------------------------------------------------------------')
- 9413    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
-
-         IF(lwp) THEN
-            WRITE (numout,*)
-            WRITE (numout,*)
-            WRITE (numout,9420) kt
-            WRITE (numout,9421) zhke( 1) / tvols
-            WRITE (numout,9422) zhke( 2) / tvols
-            WRITE (numout,9423) zhke( 3) / tvols
-            WRITE (numout,9424) zhke( 4) / tvols
-            WRITE (numout,9425) zhke( 5) / tvols
-            WRITE (numout,9426) zhke( 6) / tvols
-            WRITE (numout,9427) zhke( 7) / tvols
-            WRITE (numout,9428) zhke( 8) / tvols
-            WRITE (numout,9429) zhke(10) / tvols
-            WRITE (numout,9430) zhke( 9) / tvols
-            WRITE (numout,9431)
-            WRITE (numout,9432)   &
-            &     (  zhke(1) + zhke(2) + zhke(3) + zhke(4) + zhke(5) + zhke(6)   &
-            &      + zhke(7) + zhke(8) + zhke(9) + zhke(10) ) / tvols
-         ENDIF
-
- 9420    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
- 9421    FORMAT(' pressure gradient         u2= ', e20.13)
- 9422    FORMAT(' ke gradient               u2= ', e20.13)
- 9423    FORMAT(' relative vorticity term   u2= ', e20.13)
- 9424    FORMAT(' coriolis term             u2= ', e20.13)
- 9425    FORMAT(' horizontal diffusion      u2= ', e20.13)
- 9426    FORMAT(' vertical advection        u2= ', e20.13)
- 9427    FORMAT(' vertical diffusion        u2= ', e20.13)
- 9428    FORMAT(' surface pressure gradient u2= ', e20.13)
- 9429    FORMAT(' forcing term              u2= ', e20.13)
- 9430    FORMAT(' dampimg term              u2= ', e20.13)
- 9431    FORMAT(' --------------------------------------------------')
- 9432    FORMAT(' total trend               u2= ', e20.13)
-
-         IF(lwp) THEN
-            WRITE (numout,*)
-            WRITE (numout,*)
-            WRITE (numout,9440) kt
-            WRITE (numout,9441) ( zhke(2) + zhke(3) + zhke(6) ) / tvols
-            WRITE (numout,9442) ( zhke(2) + zhke(6) ) / tvols
-            WRITE (numout,9443) ( zhke(4) ) / tvols
-            WRITE (numout,9444) ( zhke(3) ) / tvols
-            WRITE (numout,9445) ( zhke(8) ) / tvols
-            WRITE (numout,9446) ( zhke(5) ) / tvols
-            WRITE (numout,9447) ( zhke(7) ) / tvols
-            WRITE (numout,9448) ( zhke(1) ) / tvols, rpktrd / tvols
-         ENDIF
-
- 9440    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
- 9441    FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
- 9442    FORMAT(' 0 = ke gradient + vertical advection              : ', e20.13)
- 9443    FORMAT(' 0 = coriolis term  (true if key_vorenergy or key_combined): ', e20.13)
- 9444    FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.)    : ', e20.13)
- 9445    FORMAT(' 0 = surface pressure gradient                     : ', e20.13)
- 9446    FORMAT(' 0 > horizontal diffusion                          : ', e20.13)
- 9447    FORMAT(' 0 > vertical diffusion                            : ', e20.13)
- 9448    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)        : ', e20.13, '  u.dz(rhop) =', e20.13)
-
-         ! Save potential to kinetic energy conversion for next time step
-         
-         rpktrd = zpeke
-         
-      ENDIF
-
-   END SUBROUTINE trd_dyn
-
-
-   SUBROUTINE trd_dyn_init
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_dyn_init  ***
-      !!
-      !! ** Purpose :   
-      !!
-      !! ** Method :
-      !!
-      !! History :
-      !!   9.0  !  03-09  (G. Madec)  Original code
-      !!---------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk        ! dummy loop indices
-
-      REAL(wp) ::   zmskt,zmsku,zmskv
-      
-      NAMELIST/namtrd/ ntrd, nctls
-      !!----------------------------------------------------------------------
-
-      ! 0. Initialization
-      ! -----------------
-
-      ! set the trends to zero
-      utrd  (:,:,:,:) = 0.e0
-      vtrd  (:,:,:,:) = 0.e0
-      tautrd(:,:,  :) = 0.e0
-
-      ! namelist namtrd : ocean momentum trend 
-      REWIND( numnam )
-      READ  ( numnam, namtrd )
-
-      IF(lwp) THEN
-         WRITE(numout,*) 'trd_dyn : read namelist namtrd'
-         WRITE(numout,*) '~~~~~~~'
-         WRITE(numout,*) ' time step frequency trend       ntrd  = ', ntrd
-         WRITE(numout,*) ' '
-      ENDIF
-
-#if defined key_s_coord || defined key_partial_steps
-      ! dynamics trend diagnostics not yet implemented
-      IF(lwp) WRITE(numout,*) 'trd_dyn : dynamics trend diagnostics not yet implemented'
-      nstop = nstop + 1
-#endif
-
-      ! total volume at u-, v- and t-points:
-      tvols = 0.
-      tvolu = 0.
-      tvolv = 0.
-
-      DO jk = 1, jpk
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zmskt =                      tmask_i(ji,jj) * tmask(ji,jj,jk)
-               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
-               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
-               tvols = tvols + zmskt * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-               tvolu = tvolu + zmsku * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-               tvolv = tvolv + zmskv * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      IF( lk_mpp )   CALL mpp_sum( tvols )   ! sums over the global domain
-      IF( lk_mpp )   CALL mpp_sum( tvolu )
-      IF( lk_mpp )   CALL mpp_sum( tvolv )
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'trd_dyn : frequency ntrd = ', ntrd
-         WRITE(numout,*) '~~~~~~~   tvols = ', tvols
-         WRITE(numout,*) '          tvolu = ', tvolu,' tvolv = ', tvolv
-      ENDIF
-
-      ! 0.2 Initialization of potential to kinetic energy conversion
-      
-      rpktrd = 0.e0
-
-   END SUBROUTINE trd_dyn_init
-
-#   else 
-   !!----------------------------------------------------------------------
-   !!   Default option :                      NO mementum trend diagnostics
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .FALSE.   !: momentum trend flag
-CONTAINS
-   SUBROUTINE trd_dyn( kt )        ! Empty routine
-      WRITE(*,*) 'trd_dyn: You should not have seen this print! error?', kt
-   END SUBROUTINE trd_dyn
-   SUBROUTINE trd_dyn_init         ! Empty routine
-   END SUBROUTINE trd_dyn_init
-#endif
-
-   !!======================================================================
-END MODULE trddyn
Index: trunk/NEMO/OPA_SRC/TRD/trddyn_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trddyn_oce.F90	(revision 129)
+++ 	(revision )
@@ -1,39 +1,0 @@
-MODULE trddyn_oce
-   !!======================================================================
-   !!                      ***  trddyn_oce  ***
-   !! Ocean trends :  dynamical trends arrays defined in memory 
-   !!======================================================================
-#if   defined key_trddyn   ||   defined key_esopa || defined key_trd_vor
-   !!----------------------------------------------------------------------
-   !!   'key_trddyn'                             momentum trend diagnostics
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE par_oce      ! ocean parameters
-
-   IMPLICIT NONE
-   PUBLIC
-
-    INTEGER  ::      & !!! namdia :  diagnostics on dynamics and/or tracer trends
-      ntrd=10,   &
-      nctls=0
-
-   REAL(wp) ::   &
-      tvols,     &  ! volume of the whole ocean computed at t-points
-      tvolu,     &  ! volume of the whole ocean computed at u-points
-      tvolv,     &  ! volume of the whole ocean computed at v-points
-      rpktrd        ! potential to kinetic energy conversion
-   REAL(wp), DIMENSION(jpi,jpj,jpk,9) ::   & 
-      utrd,      &  ! trends of the momentum equation
-      vtrd
-   REAL(wp), DIMENSION(jpi,jpj,2) ::   & 
-      utrd2,     &  ! trend with the free surface formulation
-      vtrd2
-   REAL(wp), DIMENSION(jpi,jpj,4) ::    &
-      tautrd        ! momentum forcing trends
-#else
-   !!----------------------------------------------------------------------
-   !!   Default case                          NO momentum trend diagnostics
-   !!----------------------------------------------------------------------
-#endif
-   !!----------------------------------------------------------------------
-END MODULE trddyn_oce
Index: trunk/NEMO/OPA_SRC/TRD/trdicp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 215)
+++ trunk/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 215)
@@ -0,0 +1,813 @@
+MODULE trdicp
+   !!======================================================================
+   !!                       ***  MODULE  trdicp  ***
+   !! Ocean diagnostics:  ocean tracers and dynamic trends
+   !!=====================================================================
+#if  defined key_trdtra   ||   defined key_trddyn   ||   defined key_esopa
+   !!----------------------------------------------------------------------
+   !!   'key_trdtra'  or                  active tracers trends diagnostics
+   !!   'key_trddyn'                            momentum trends diagnostics
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   trd              : verify the basin averaged properties for tra/dyn 
+   !!   trd_dwr          : print dynmaic trends in ocean.output file
+   !!   trd_twr          : print tracers trends in ocean.output file
+   !!   trd_icp_init     : initialization step
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables
+   USE trdmod_oce      ! ocean variables trends
+   USE ldftra_oce      ! ocean active tracers: lateral physics
+   USE ldfdyn_oce      ! ocean dynamics: lateral physics
+   USE zdf_oce         ! ocean vertical physics
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distibuted memory computing library
+   USE eosbn2          ! equation of state
+   USE phycst          ! physical constants
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Interfaces
+   INTERFACE trd
+      MODULE PROCEDURE trd_2d, trd_3d
+   END INTERFACE
+
+   !! * Routine accessibility
+   PUBLIC trd                   ! called by step.F90
+   PUBLIC trd_dwr               ! called by step.F90
+   PUBLIC trd_twr               ! called by step.F90
+   PUBLIC trd_icp_init          ! called by opa.F90
+
+   !! * Shared module variables
+#if  defined key_trdtra   &&   defined key_trddyn
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trdtra = .TRUE.    !: tracers  trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .TRUE.    !: momentum trend flag
+#elif  defined key_trdtra
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trdtra = .TRUE.    !: tracers  trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .FALSE.   !: momentum trend flag
+#elif  defined key_trddyn
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trdtra = .FALSE.   !: tracers  trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .TRUE.    !: momentum trend flag
+#endif
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , LODYC-IPSL  (2004)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trd_2d(ptrd2dx, ptrd2dy, ktrd , ctype)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_2d  ***
+      !! 
+      !! ** Purpose : verify the basin averaged properties of tracers and/or
+      !!              momentum equations at every time step frequency ntrd.
+      !!
+      !! ** Method :
+      !!
+      !! History :
+      !!        !  91-12 (G. Madec)
+      !!        !  92-06 (M. Imbard) add time step frequency
+      !!        !  96-01 (G. Madec)  terrain following coordinates
+      !!   8.5  !  02-06 (G. Madec)  F90: Free form and module
+      !!   9.0  !  04-08 (C. Talandier) New trends organization
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   &
+         ptrd2dx,                      &   ! Temperature or U trend 
+         ptrd2dy                           ! Salinity    or V trend
+
+      INTEGER, INTENT( in ) ::   ktrd      ! tracer trend index
+
+      CHARACTER(len=3), INTENT( in ) ::   &
+         ctype                             ! momentum or tracers trends type
+         !                                 ! 'DYN' or 'TRA'
+
+      !! * Local declarations
+      INTEGER ::   ji, jj        ! loop indices
+      REAL(wp) ::   &
+         zbt, zbtu, zbtv,     &  ! temporary scalars
+         zmsku, zmskv            !    "         "
+      !!----------------------------------------------------------------------
+
+      ! 1. Advective trends and forcing trend
+      ! -------------------------------------
+
+      ! 1.1 Mask the forcing trend and substract it from the vertical diffusion trend
+      SELECT CASE (ctype)
+
+      CASE ('DYN')              ! Momentum
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,1)
+               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,1)
+               ptrd2dx(ji,jj) = ptrd2dx(ji,jj) * zmsku
+               ptrd2dy(ji,jj) = ptrd2dy(ji,jj) * zmskv
+            END DO
+         END DO
+         ptrd2dx(jpi, : ) = 0.e0      ;      ptrd2dy(jpi, : ) = 0.e0
+         ptrd2dx( : ,jpj) = 0.e0      ;      ptrd2dy( : ,jpj) = 0.e0
+
+      CASE ('TRA')              ! Tracers
+         ptrd2dx(:,:) = ptrd2dx(:,:) * tmask_i(:,:)
+         ptrd2dy(:,:) = ptrd2dy(:,:) * tmask_i(:,:)
+
+      END SELECT
+      
+      ! 2. Basin averaged tracer trend
+      ! ------------------------------
+
+      SELECT CASE (ctype)
+
+      CASE ('DYN')              ! Momentum
+         umo(ktrd) = 0.e0
+         vmo(ktrd) = 0.e0
+
+         SELECT CASE (ktrd)
+
+         CASE (jpdtdswf)         ! surface forcing
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  umo(ktrd) = umo(ktrd) + ptrd2dx(ji,jj) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
+                  vmo(ktrd) = vmo(ktrd) + ptrd2dy(ji,jj) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
+               END DO
+            END DO
+
+         CASE (jpdtdbfr)         ! bottom friction fluxes
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  umo(ktrd) = umo(ktrd) + ptrd2dx(ji,jj)
+                  vmo(ktrd) = vmo(ktrd) + ptrd2dy(ji,jj)
+               END DO
+            END DO
+
+         END SELECT
+
+      CASE ('TRA')              ! Tracers
+         tmo(ktrd) = 0.e0
+         smo(ktrd) = 0.e0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
+               tmo(ktrd) =  tmo(ktrd) + ptrd2dx(ji,jj) * zbt
+               smo(ktrd) =  smo(ktrd) + ptrd2dy(ji,jj) * zbt
+            END DO
+         END DO
+
+      END SELECT
+      
+      ! 3. Basin averaged tracer square trend
+      ! -------------------------------------
+      ! c a u t i o n: field now
+      
+      SELECT CASE (ctype)
+
+      CASE ('DYN')              ! Momentum
+         hke(ktrd) = 0.e0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
+               zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
+               hke(ktrd) = hke(ktrd)   &
+               &   + un(ji,jj,1) * ptrd2dx(ji,jj) * zbtu &
+               &   + vn(ji,jj,1) * ptrd2dy(ji,jj) * zbtv
+            END DO
+         END DO
+
+      CASE ('TRA')              ! Tracers
+         t2(ktrd) = 0.e0
+         s2(ktrd) = 0.e0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
+               t2(ktrd) = t2(ktrd) + ptrd2dx(ji,jj) * zbt * tn(ji,jj,1)
+               s2(ktrd) = s2(ktrd) + ptrd2dy(ji,jj) * zbt * sn(ji,jj,1)
+            END DO
+         END DO
+      
+      END SELECT
+
+   END SUBROUTINE trd_2d
+
+
+
+   SUBROUTINE trd_3d(ptrd3dx, ptrd3dy, ktrd, ctype)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_3d  ***
+      !! 
+      !! ** Purpose : verify the basin averaged properties of tracers and/or
+      !!              momentum equations at every time step frequency ntrd.
+      !!
+      !! ** Method :
+      !!
+      !! History :
+      !!        !  91-12 (G. Madec)
+      !!        !  92-06 (M. Imbard) add time step frequency
+      !!        !  96-01 (G. Madec)  terrain following coordinates
+      !!   8.5  !  02-06 (G. Madec)  F90: Free form and module
+      !!   9.0  !  04-08 (C. Talandier) New trends organization
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   &
+          ptrd3dx,                     &   ! Temperature or U trend 
+          ptrd3dy                          ! Salinity    or V trend
+
+      INTEGER, INTENT( in ) ::   ktrd      ! momentum or tracer trend index
+
+      CHARACTER(len=3), INTENT( in ) ::   &
+         ctype                             ! momentum or tracers trends type
+         !                                 ! 'DYN' or 'TRA'
+
+      !! * Local declarations
+      INTEGER ::   ji, jj, jk
+      REAL(wp) ::   &
+         zbt, zbtu, zbtv,               &  ! temporary scalars
+         zmsku, zmskv
+      !!----------------------------------------------------------------------
+
+      ! 1. Advective trends and forcing trend
+      ! -------------------------------------
+
+      ! Mask the trends
+      SELECT CASE (ctype)
+
+      CASE ('DYN')              ! Momentum        
+         DO jk = 1, jpk
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
+                  zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
+                  ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * zmsku
+                  ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * zmskv
+               ENDDO
+            ENDDO
+         ENDDO
+
+         ptrd3dx(jpi, : ,:) = 0.e0      ;      ptrd3dy(jpi, : ,:) = 0.e0
+         ptrd3dx( : ,jpj,:) = 0.e0      ;      ptrd3dy( : ,jpj,:) = 0.e0
+
+      CASE ('TRA')              ! Tracers
+         DO jk = 1, jpk
+            ptrd3dx(:,:,jk) = ptrd3dx(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
+            ptrd3dy(:,:,jk) = ptrd3dy(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
+         ENDDO
+
+      END SELECT   
+
+      ! 2. Basin averaged tracer/momentum trend
+      ! ---------------------------------------
+      
+      SELECT CASE (ctype)
+
+      CASE ('DYN')              ! Momentum
+         umo(ktrd) = 0.e0
+         vmo(ktrd) = 0.e0
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
+                  zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
+                  umo(ktrd) = umo(ktrd) + ptrd3dx(ji,jj,jk) * zbtu
+                  vmo(ktrd) = vmo(ktrd) + ptrd3dy(ji,jj,jk) * zbtv
+               END DO
+            END DO
+         END DO
+
+      CASE ('TRA')              ! Tracers
+         tmo(ktrd) = 0.e0
+         smo(ktrd) = 0.e0
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) 
+                  tmo(ktrd) =  tmo(ktrd) + ptrd3dx(ji,jj,jk) * zbt
+                  smo(ktrd) =  smo(ktrd) + ptrd3dy(ji,jj,jk) * zbt
+               END DO
+            END DO
+         END DO
+
+      END SELECT
+
+      ! 3. Basin averaged tracer/momentum square trend
+      ! ----------------------------------------------
+      ! c a u t i o n: field now
+      
+      SELECT CASE (ctype)
+
+      CASE ('DYN')              ! Momentum
+         hke(ktrd) = 0.e0
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
+                  zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
+                  hke(ktrd) = hke(ktrd)   &
+                  &   + un(ji,jj,jk) * ptrd3dx(ji,jj,jk) * zbtu &
+                  &   + vn(ji,jj,jk) * ptrd3dy(ji,jj,jk) * zbtv
+               END DO
+            END DO
+         END DO
+
+      CASE ('TRA')              ! Tracers
+         t2(ktrd) = 0.e0
+         s2(ktrd) = 0.e0
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
+                  t2(ktrd) = t2(ktrd) + ptrd3dx(ji,jj,jk) * zbt * tn(ji,jj,jk)
+                  s2(ktrd) = s2(ktrd) + ptrd3dy(ji,jj,jk) * zbt * sn(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+
+      END SELECT
+
+   END SUBROUTINE trd_3d
+
+
+
+   SUBROUTINE trd_icp_init
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_icp_init  ***
+      !! 
+      !! ** Purpose :   
+      !!
+      !! ** Method  :
+      !!
+      !! History :
+      !!   9.0  !  03-09 (G. Madec)  Original code
+      !!        !  04-08 (C. Talandier) New trends organization
+      !!----------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER :: ji, jj, jk
+
+      REAL(wp) ::   zmskt
+#if  defined key_trddyn
+      REAL(wp) ::   zmsku,zmskv
+#endif
+
+      NAMELIST/namtrd/ ntrd, nctls
+      !!----------------------------------------------------------------------
+
+      ! namelist namtrd : trend diagnostic
+      REWIND( numnam )
+      READ  ( numnam, namtrd )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'trd_icp_init : integral constraints properties trends'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+         WRITE(numout,*) ' '
+         WRITE(numout,*) '          Namelist namtrd : '
+         WRITE(numout,*) '             time step frequency trend       ntrd  = ', ntrd
+      ENDIF
+
+      ! initialisation of BBL tracers lateral diffusion to zero
+      tldfbbl(:,:) = 0.e0   ;   sldfbbl(:,:) = 0.e0  
+      ! initialisation of BBL tracers lateral advection to zero
+      tladbbl(:,:) = 0.e0   ;   sladbbl(:,:) = 0.e0  
+      ! initialisation of workspace
+      tladi(:,:,:) = 0.e0  ;  tladj(:,:,:) = 0.e0
+      sladi(:,:,:) = 0.e0  ;  sladj(:,:,:) = 0.e0
+
+      ! Total volume at t-points:
+      tvolt = 0.e0
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zmskt = tmask(ji,jj,jk) * tmask_i(ji,jj)
+               tvolt = tvolt + zmskt * e1t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '          total ocean volume at T-point   tvolt = ',tvolt
+      ENDIF
+
+#if  defined key_trddyn
+      ! Initialization of potential to kinetic energy conversion
+      rpktrd = 0.e0
+
+      ! Total volume at u-, v- points:
+      tvolu = 0.e0
+      tvolv = 0.e0
+
+      DO jk = 1, jpk
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
+               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
+               tvolu = tvolu + zmsku * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
+               tvolv = tvolv + zmskv * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
+      IF( lk_mpp )   CALL mpp_sum( tvolv )
+
+      IF(lwp) THEN
+         WRITE(numout,*) '          total ocean volume at U-point   tvolu = ',tvolu
+         WRITE(numout,*) '          total ocean volume at V-point   tvolv = ',tvolv
+         WRITE(numout,*) ' '
+      ENDIF
+#endif
+
+   END SUBROUTINE trd_icp_init
+
+
+
+   SUBROUTINE trd_dwr( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_dwr  ***
+      !! 
+      !! ** Purpose :  write dynamic trends in ocean.output 
+      !!
+      !! ** Method  :
+      !!
+      !! History :
+      !!   9.0  !  03-09  (G. Madec)  Original code
+      !!        !  04-08  (C. Talandier)  New trends organization
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER ::   ji, jj, jk
+      REAL(wp) ::   &
+         ze1e2w,zcof,        &  !    "         "
+         zbe1ru, zbe2rv,     &  !    "         "
+         zbtr, ztz, zth 
+
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
+               zkepe, zkx, zky, zkz              ! temporary arrays
+      !!----------------------------------------------------------------------
+
+      ! I. Momentum trends
+      ! -------------------
+
+      IF( MOD(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
+
+         ! I.1 Conversion potential energy - kinetic energy
+         ! --------------------------------------------------
+         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
+
+         zkx(:,:,:) = 0.e0
+         zky(:,:,:) = 0.e0
+         zkz(:,:,:) = 0.e0
+         zkepe(:,:,:) = 0.e0
+   
+         CALL eos( tn, sn, rhd, rhop )       ! now potential and in situ densities
+
+         ! 4.1 Density flux at w-point
+         DO jk = 2, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ze1e2w = 0.5 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) * tmask_i(ji,jj)
+                  zkz(ji,jj,jk) = ze1e2w / rau0 * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) )
+               END DO
+            END DO
+         END DO
+         zkz  (:,:, 1 ) = 0.e0
+         
+         ! Density flux at u and v-points
+         DO jk = 1, jpk
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  zcof   = 0.5 / rau0
+                  zbe1ru = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
+                  zbe2rv = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
+                  zkx(ji,jj,jk) = zbe1ru * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
+                  zky(ji,jj,jk) = zbe2rv * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
+               END DO
+            END DO
+         END DO
+         
+         ! Density flux divergence at t-point
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
+                  ztz = - zbtr * (    zkz(ji,jj,jk) - zkz(ji,jj,jk+1) )
+                  zth = - zbtr * (  ( zkx(ji,jj,jk) - zkx(ji-1,jj,jk) )   &
+                    &             + ( zky(ji,jj,jk) - zky(ji,jj-1,jk) )  )
+                  zkepe(ji,jj,jk) = (zth + ztz) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+               END DO
+            END DO
+         END DO
+         zkepe( : , : ,jpk) = 0.e0
+         zkepe( : ,jpj, : ) = 0.e0
+         zkepe(jpi, : , : ) = 0.e0
+
+         ! I.2 Basin averaged kinetic energy trend
+         ! ----------------------------------------
+         peke = 0.e0
+         DO jk = 1,jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  peke = peke + zkepe(ji,jj,jk) * grav * fsdept(ji,jj,jk)   &
+                     &                     * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+
+         ! I.3 Sums over the global domain
+         ! ---------------------------------
+         IF( lk_mpp ) THEN
+               CALL mpp_sum( peke )
+               CALL mpp_sum( umo , 11 )
+               CALL mpp_sum( vmo , 11 )
+               CALL mpp_sum( hke , 10 )
+         END IF
+
+         ! I.2 Print dynamic trends in the ocean.output file
+         ! --------------------------------------------------
+
+         IF(lwp) THEN
+            WRITE (numout,*)
+            WRITE (numout,*)
+            WRITE (numout,9500) kt
+            WRITE (numout,9501) umo( 1) / tvolu, vmo( 1) / tvolv
+            WRITE (numout,9502) umo( 2) / tvolu, vmo( 2) / tvolv
+            WRITE (numout,9503) umo( 3) / tvolu, vmo( 3) / tvolv
+            WRITE (numout,9504) umo( 4) / tvolu, vmo( 4) / tvolv
+            WRITE (numout,9505) umo( 5) / tvolu, vmo( 5) / tvolv
+            WRITE (numout,9506) umo( 6) / tvolu, vmo( 6) / tvolv
+            WRITE (numout,9507) umo( 7) / tvolu, vmo( 7) / tvolv
+            WRITE (numout,9508) umo( 8) / tvolu, vmo( 8) / tvolv
+            WRITE (numout,9509) umo(10) / tvolu, vmo(10) / tvolv
+            WRITE (numout,9510) umo( 9) / tvolu, vmo( 9) / tvolv
+            WRITE (numout,9511) umo(11) / tvolu, vmo(11) / tvolv
+            WRITE (numout,9512)
+            WRITE (numout,9513)                                                 &
+            &     (  umo(1) + umo(2) + umo(3) + umo( 4) + umo( 5) + umo(6)    &
+            &      + umo(7) + umo(8) + umo(9) + umo(10) + umo(11) ) / tvolu,   &
+            &     (  vmo(1) + vmo(2) + vmo(3) + vmo( 4) + vmo( 5) + vmo(6)    &
+            &      + vmo(7) + vmo(8) + vmo(9) + vmo(10) + vmo(11) ) / tvolv
+         ENDIF
+
+ 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
+ 9501    FORMAT(' pressure gradient          u= ', e20.13, '    v= ', e20.13)
+ 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
+ 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
+ 9504    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
+ 9505    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
+ 9506    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
+ 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
+ 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
+ 9509    FORMAT(' forcing term               u= ', e20.13, '    v= ', e20.13)
+ 9510    FORMAT(' dampimg term               u= ', e20.13, '    v= ', e20.13)
+ 9511    FORMAT(' bottom flux                u= ', e20.13, '    v= ', e20.13)
+ 9512    FORMAT(' -----------------------------------------------------------------------------')
+ 9513    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
+
+         IF(lwp) THEN
+            WRITE (numout,*)
+            WRITE (numout,*)
+            WRITE (numout,9520) kt
+            WRITE (numout,9521) hke( 1) / tvolt
+            WRITE (numout,9522) hke( 2) / tvolt
+            WRITE (numout,9523) hke( 3) / tvolt
+            WRITE (numout,9524) hke( 4) / tvolt
+            WRITE (numout,9525) hke( 5) / tvolt
+            WRITE (numout,9526) hke( 6) / tvolt
+            WRITE (numout,9527) hke( 7) / tvolt
+            WRITE (numout,9528) hke( 8) / tvolt
+            WRITE (numout,9529) hke(10) / tvolt
+            WRITE (numout,9530) hke( 9) / tvolt
+            WRITE (numout,9531)
+            WRITE (numout,9532)   &
+            &     (  hke(1) + hke(2) + hke(3) + hke(4) + hke(5) + hke(6)   &
+            &      + hke(7) + hke(8) + hke(9) + hke(10) ) / tvolt
+         ENDIF
+
+ 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
+ 9521    FORMAT(' pressure gradient         u2= ', e20.13)
+ 9522    FORMAT(' ke gradient               u2= ', e20.13)
+ 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
+ 9524    FORMAT(' coriolis term             u2= ', e20.13)
+ 9525    FORMAT(' horizontal diffusion      u2= ', e20.13)
+ 9526    FORMAT(' vertical advection        u2= ', e20.13)
+ 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
+ 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
+ 9529    FORMAT(' forcing term              u2= ', e20.13)
+ 9530    FORMAT(' dampimg term              u2= ', e20.13)
+ 9531    FORMAT(' --------------------------------------------------')
+ 9532    FORMAT(' total trend               u2= ', e20.13)
+
+         IF(lwp) THEN
+            WRITE (numout,*)
+            WRITE (numout,*)
+            WRITE (numout,9540) kt
+            WRITE (numout,9541) ( hke(2) + hke(3) + hke(6) ) / tvolt
+            WRITE (numout,9542) ( hke(2) + hke(6) ) / tvolt
+            WRITE (numout,9543) ( hke(4) ) / tvolt
+            WRITE (numout,9544) ( hke(3) ) / tvolt
+            WRITE (numout,9545) ( hke(8) ) / tvolt
+            WRITE (numout,9546) ( hke(5) ) / tvolt
+            WRITE (numout,9547) ( hke(7) ) / tvolt
+            WRITE (numout,9548) ( hke(1) ) / tvolt, rpktrd / tvolt
+         ENDIF
+
+ 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
+ 9541    FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
+ 9542    FORMAT(' 0 = ke gradient + vertical advection              : ', e20.13)
+ 9543    FORMAT(' 0 = coriolis term  (true if key_vorenergy or key_combined): ', e20.13)
+ 9544    FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.)    : ', e20.13)
+ 9545    FORMAT(' 0 = surface pressure gradient                     : ', e20.13)
+ 9546    FORMAT(' 0 > horizontal diffusion                          : ', e20.13)
+ 9547    FORMAT(' 0 > vertical diffusion                            : ', e20.13)
+ 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)        : ', e20.13, '  u.dz(rhop) =', e20.13)
+
+         ! Save potential to kinetic energy conversion for next time step
+         rpktrd = peke
+
+      ENDIF
+
+   END SUBROUTINE trd_dwr
+
+
+
+
+   SUBROUTINE trd_twr( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_twr  ***
+      !! 
+      !! ** Purpose :  write active tracers trends in ocean.output 
+      !!
+      !! ** Method  :
+      !!
+      !! History :
+      !!   9.0  !  03-09  (G. Madec)  Original code
+      !!        !  04-08  (C. Talandier)  New trends organization
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+
+      !!----------------------------------------------------------------------
+
+      ! I. Tracers trends
+      ! -----------------
+
+      IF( MOD(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
+
+         ! I.1 Sums over the global domain
+         ! -------------------------------
+         IF( lk_mpp ) THEN
+            CALL mpp_sum( tmo, 10 )   
+            CALL mpp_sum( smo, 10 )
+            CALL mpp_sum( t2 , 10 )
+            CALL mpp_sum( s2 , 10 )
+         ENDIF
+
+         ! I.2 Print tracers trends in the ocean.output file
+         ! --------------------------------------------------
+         
+         IF(lwp) THEN
+            WRITE (numout,*)
+            WRITE (numout,*)
+            WRITE (numout,9400) kt
+            WRITE (numout,9401) tmo(1) / tvolt, smo(1) / tvolt
+            WRITE (numout,9402) tmo(2) / tvolt, smo(2) / tvolt
+            WRITE (numout,9403) tmo(3) / tvolt, smo(3) / tvolt
+            WRITE (numout,9404) tmo(4) / tvolt, smo(4) / tvolt
+            WRITE (numout,9405) tmo(5) / tvolt, smo(5) / tvolt
+            WRITE (numout,9406) tmo(6) / tvolt, smo(6) / tvolt
+            WRITE (numout,9407) tmo(7) / tvolt
+            WRITE (numout,9408) tmo(8) / tvolt, smo(8) / tvolt
+            WRITE (numout,9409)
+            WRITE (numout,9410) (  tmo(1) + tmo(2) + tmo(3) + tmo(4)              &
+            &                    + tmo(5) + tmo(6) + tmo(7) + tmo(8) ) / tvolt,   &
+            &                   (  smo(1) + smo(2) + smo(3) + smo(4)              &
+            &                    + smo(5) + smo(6)           + smo(8) ) / tvolt
+         ENDIF
+
+9400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
+              '              salinity',/' ============================')
+9401     FORMAT(' horizontal advection        ',e20.13,'     ',e20.13)
+9402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
+9403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
+9404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
+9405     FORMAT(' STATIC instability mixing   ',e20.13,'     ',e20.13)
+9406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
+9407     FORMAT(' penetrative qsr             ',e20.13,'     ',e20.13)
+9408     FORMAT(' forcing term                ',e20.13,'     ',e20.13)
+9409     FORMAT(' -------------------------------------------------------------------------')
+9410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
+
+
+         IF(lwp) THEN
+            WRITE (numout,*)
+            WRITE (numout,*)
+            WRITE (numout,9420) kt
+            WRITE (numout,9421) t2(1) / tvolt, s2(1) / tvolt
+            WRITE (numout,9422) t2(2) / tvolt, s2(2) / tvolt
+            WRITE (numout,9423) t2(3) / tvolt, s2(3) / tvolt
+            WRITE (numout,9424) t2(4) / tvolt, s2(4) / tvolt
+            WRITE (numout,9425) t2(5) / tvolt, s2(5) / tvolt
+            WRITE (numout,9426) t2(6) / tvolt, s2(6) / tvolt
+            WRITE (numout,9427) t2(7) / tvolt
+            WRITE (numout,9428) t2(8) / tvolt, s2(8) / tvolt
+            WRITE (numout,9429)
+            WRITE (numout,9430) (  t2(1) + t2(2) + t2(3) + t2(4)              &
+            &                    + t2(5) + t2(6) + t2(7) + t2(8) ) / tvolt,   &
+            &                   (  s2(1) + s2(2) + s2(3) + s2(4)              &
+            &                    + s2(5) + s2(6)          + s2(8) ) / tvolt
+         ENDIF
+
+9420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
+            '               salinity', /, ' ===============================')
+9421     FORMAT(' horizontal advection      * t   ', e20.13, '     ', e20.13)
+9422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
+9423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
+9424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
+9425     FORMAT(' STATIC instability mixing * t   ', e20.13, '     ', e20.13)
+9426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
+9427     FORMAT(' penetrative qsr           * t   ', e20.13, '     ', e20.13)
+9428     FORMAT(' forcing term              * t   ', e20.13, '     ', e20.13)
+9429     FORMAT(' -----------------------------------------------------------------------------')
+9430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
+
+
+         IF(lwp) THEN
+            WRITE (numout,*)
+            WRITE (numout,*)
+            WRITE (numout,9440) kt
+            WRITE (numout,9441) ( tmo(1)+tmo(2) )/tvolt, ( smo(1)+smo(2) )/tvolt
+            WRITE (numout,9442)   tmo(3)/tvolt,  smo(3)/tvolt
+            WRITE (numout,9443)   tmo(4)/tvolt,  smo(4)/tvolt
+            WRITE (numout,9444)   tmo(5)/tvolt,  smo(5)/tvolt
+            WRITE (numout,9445) ( t2(1)+t2(2) )/tvolt, ( s2(1)+s2(2) )/tvolt
+            WRITE (numout,9446)   t2(3)/tvolt,   s2(3)/tvolt
+            WRITE (numout,9447)   t2(4)/tvolt,   s2(4)/tvolt
+            WRITE (numout,9448)   t2(5)/tvolt,   s2(5)/tvolt
+         ENDIF
+
+9440     FORMAT(' tracer consistency at it= ',i6,   &
+            ' :         temperature','                salinity',/,   &
+            ' ==================================')
+9441     FORMAT(' 0 = horizontal+vertical advection      ',e20.13,'       ',e20.13)
+9442     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
+9443     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
+9444     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
+9445     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
+9446     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
+9447     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
+9448     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
+
+      ENDIF
+
+   END SUBROUTINE trd_twr
+
+#   else
+   !!----------------------------------------------------------------------
+   !!   Default case :                                         Empty module
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trdtra = .FALSE.   !: tracers  trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .FALSE.   !: momentum trend flag
+CONTAINS
+   SUBROUTINE trd_2d(ptrd2dx, ptrd2dy, ktrd , ctype)       ! Empty routine
+      REAL, DIMENSION(:,:,:), INTENT( inout ) ::   &
+          ptrd2dx,                     &   ! Temperature or U trend 
+          ptrd2dy                          ! Salinity    or V trend
+      INTEGER, INTENT( in ) ::   ktrd      ! momentum or tracer trend index
+      CHARACTER(len=3), INTENT( in ) ::   &
+         ctype                             ! momentum or tracers trends type
+      WRITE(*,*) 'trd_2d: You should not have seen this print! error ?', ptrd2dx(1,1,1)
+      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ptrd2dy(1,1,1)
+      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ktrd
+      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ctype
+   END SUBROUTINE trd_2d
+   SUBROUTINE trd_3d(ptrd3dx, ptrd3dy, ktrd , ctype)       ! Empty routine
+      REAL, DIMENSION(:,:,:), INTENT( inout ) ::   &
+          ptrd3dx,                     &   ! Temperature or U trend 
+          ptrd3dy                          ! Salinity    or V trend
+      INTEGER, INTENT( in ) ::   ktrd      ! momentum or tracer trend index
+      CHARACTER(len=3), INTENT( in ) ::   &
+         ctype                             ! momentum or tracers trends type
+      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', ptrd3dx(1,1,1)
+      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ptrd3dy(1,1,1)
+      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ktrd
+      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ctype
+   END SUBROUTINE trd_3d
+   SUBROUTINE trd_icp_init               ! Empty routine
+   END SUBROUTINE trd_icp_init
+   SUBROUTINE trd_dwr( kt )          ! Empty routine
+      INTEGER, INTENT(in) :: kt
+      WRITE(*,*) 'trd_dwr: You should not have seen this print! error ?', kt
+   END SUBROUTINE trd_dwr
+   SUBROUTINE trd_twr( kt )          ! Empty routine
+      INTEGER, INTENT(in) :: kt
+      WRITE(*,*) 'trd_twr: You should not have seen this print! error ?', kt
+   END SUBROUTINE trd_twr
+#endif
+
+   !!======================================================================
+END MODULE trdicp
Index: trunk/NEMO/OPA_SRC/TRD/trdicp_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdicp_oce.F90	(revision 215)
+++ trunk/NEMO/OPA_SRC/TRD/trdicp_oce.F90	(revision 215)
@@ -0,0 +1,121 @@
+MODULE trdicp_oce
+   !!======================================================================
+   !!                   ***  MODULE trdicp_oce  ***
+   !! Ocean trends :   set tracer and momentum trend variables
+   !!======================================================================
+
+   !!----------------------------------------------------------------------
+   !!   'key_trdtra'   or                         tracer trends diagnostics
+   !!   'key_trddyn'                            momentum trends diagnostics
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE par_oce                 ! ocean parameters
+
+   IMPLICIT NONE
+   PUBLIC
+
+   !! Namelist parameters
+   !!----------------------------------------------------------------------
+   INTEGER  ::      & !!: namdia :  diagnostics on dynamics and/or tracer trends
+      ntrd  = 10 ,  &  !: time step frequency dynamics and tracers trends
+      nctls =  0       !: control surface type for trends vertical integration
+
+   !! Tracers trends diagnostics parameters
+   !!---------------------------------------------------------------------
+   INTEGER, PARAMETER ::            &  !: trends index
+      jpttdlad = 1,   &  !: tracer horizontal advection
+      jpttdzad = 2,   &  !: tracer vertical advection
+      jpttdldf = 3,   &  !: tracer horizontal diffusion
+      jpttdzdf = 4,   &  !: tracer vertical diffusion
+      jpttdnpc = 5,   &  !: tracer non penetrative convection
+      jpttddoe = 6,   &  !: tracer D.amping O.r vertical E.iv
+      jpttdqsr = 7,   &  !: tracer penetrative solar radiation
+      jpttdnsr = 8       !: tracer non solar radiation
+
+   !! Momentum trends diagnostics parameters
+   !!---------------------------------------------------------------------
+   INTEGER, PARAMETER ::            &  !: trends index
+      jpdtdhpg =  1,   &  !: dynamic hydrostatic pressure gradient 
+      jpdtdkeg =  2,   &  !: dynamic kinetic energy gradient
+      jpdtdrvo =  3,   &  !: dynamic relative vorticity
+      jpdtdpvo =  4,   &  !: dynamic planetary vorticity
+      jpdtdldf =  5,   &  !: dynamic lateral diffusion
+      jpdtdzad =  6,   &  !: dynamic vertical advection
+      jpdtdzdf =  7,   &  !: dynamic vertical diffusion
+      jpdtdspg =  8,   &  !: dynamic surface pressure gradient
+      jpdtddat =  9,   &  !: dynamic damping term
+      jpdtdswf = 10,   &  !: dynamic surface wind forcing
+      jpdtdbfr = 11       !: dynamic bottom friction 
+
+   REAL, DIMENSION(jpi,jpj) ::   &  !:
+      tldfbbl, sldfbbl,          &  ! Temperature/salinity lateral diffusion trends
+      !                             ! in the BBL  
+      tladbbl, sladbbl              ! Temperature/salinity lateral advection trends 
+      !                             ! in the BBL
+
+   REAL, DIMENSION(jpi,jpj,jpk) ::   &  !:
+      tladi, sladi,                  &  ! Temp./sal. MUSCL OR TVD advection fluxes 
+      !                                 ! terms along i- 
+      tladj, sladj                      ! Temp./sal. MUSCL OR TVD advection fluxes 
+      !                                 ! terms along j- 
+#if defined key_ldfslp
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &  !:
+      uldftrd, vldftrd     !: lateral diffusion trend in isopycnal case
+#endif
+#if   defined key_trdtra   ||   defined key_trddyn   ||   defined key_esopa
+
+   !! Variables used for diagnostics
+   !!---------------------------------------------------------------------
+   REAL(wp) ::   &  !:
+      tvolt,     &  !: volume of the whole ocean computed at t-points
+      tvolu,     &  !: volume of the whole ocean computed at u-points
+      tvolv         !: volume of the whole ocean computed at v-points
+
+   !! Tracers trends diagnostics variables
+   !!---------------------------------------------------------------------
+   REAL(wp), DIMENSION(10) ::   &  !:
+      tmo, smo         !: tracers trends average 
+      !                !  tmo(1) : horizontal advection
+      !                !  tmo(2) : vertical advection
+      !                !  tmo(3) : horizontal diffusion
+      !                !  tmo(4) : vertical diffusion
+      !                !  tmo(5) : static instability
+      !                !  tmo(6) : damping OR vertical EIV
+      !                !  tmo(7) : penetrative solar radiation (T only)
+   REAL(wp), DIMENSION(10) ::   &  !:
+      t2, s2           !: tracers square trends average 
+      !                !  t2(1) : horizontal advection
+      !                !  t2(2) : vertical advection
+      !                !  t2(3) : horizontal diffusion
+      !                !  t2(4) : vertical diffusion
+      !                !  t2(5) : static instability
+      !                !  t2(6) : damping OR vertical EIV
+      !                !  t2(7) : penetrative solar radiation (T only)
+   
+   !! Momentum trends diagnostics variables
+   !!---------------------------------------------------------------------
+   REAL(wp), DIMENSION(11) ::   &  !:
+      umo, vmo         !: momentum trends average 
+      !                !  umo(1) : hydrostatic pressure gradient
+      !                !  umo(2) : kinetic energy
+      !                !  umo(3) : lateral diffusion geo-pot
+      !                !  umo(4) : 
+      !                !  umo(5) : lateral diffusion
+      !                !  umo(6) : vertical advection
+      !                !  umo(7) : vertical diffusion
+      !                !  umo(8) : surface pressure gradient
+      !                !  umo(9) : 
+
+   REAL(wp), DIMENSION(10) ::   &  !:
+      hke              !: momentum square trends average 
+      !                !  hke(1) : horizontal advection
+      !                !  hke(2) : vertical advection
+
+   REAL(wp) ::   &  !:
+      rpktrd,    &  !: potential to kinetic energy conversion
+      peke          !: conversion potential energy - kinetic energy trend
+
+#endif
+
+  !!======================================================================
+END MODULE trdicp_oce
Index: trunk/NEMO/OPA_SRC/TRD/trdmld_ncinit.h90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdmld_ncinit.h90	(revision 129)
+++ 	(revision )
@@ -1,121 +1,0 @@
-      !-----------------------------------------
-      ! II.1 Define frequency of output and means
-      ! -----------------------------------------
-
-#if defined key_diainstant
-      zsto = nwrite*rdt
-      clop ="inst(only(x))"
-#else
-      zsto = rdt
-      clop ="ave(only(x))"
-#endif
-      zout = nwrite*rdt
-
-      IF( kt == nit000 ) THEN
-
-         IF(lwp) WRITE (numout,*) ' Diamld: netCDF initialization'
-
-         ! II.2 Compute julian date from starting date of the run
-         ! ------------------------
-
-         CALL ymds2ju( nyear, nmonth, nday, 0.e0, zjulian )
-         IF (lwp) WRITE(numout,*)' '  
-         IF (lwp) WRITE(numout,*)' Date 0 used :',nit000   &
-              ,' YEAR ', nyear,' MONTH ', nmonth,' DAY ', nday   &
-              ,'Julian day : ', zjulian
-
-
-         ! II.3 Define the T grid trend file (nidtrd)
-         ! ---------------------------------
-
-         CALL dia_nam( clhstnam, nwrite, 'trends' )                  ! filename
-         IF(lwp) WRITE(numout,*) ' Name of NETCDF file ', clhstnam
-         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,1, jpi,   &  ! Horizontal grid : glamt and gphit
-              1, jpj, 0, zjulian, rdt, nh_t, nidtrd)
-         IF( isum > 0 ) THEN                                         ! Index of ocean points (2D only)
-            WRITE(numout,*)' Number of invalid points nmld > jpktrd', isum 
-            CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )   ! volume
-         ELSE 
-            CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )    ! surface
-         ENDIF
-         
-
-         ! Declare output fields as netCDF variables
-
-         ! Mixed layer Depth
-         CALL histdef( nidtrd, "somlttml", clmxl//"Depth"              , "m"   ,   &  ! hmlp
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-
-         ! Temperature
-         CALL histdef( nidtrd, "somltemp", clmxl//"Temperature"        , "C"   ,   &  ! ???
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         ! Temperature trends
-         CALL histdef( nidtrd, "somlttto", clmxl//"T Total"             , "C/s",   &  ! total
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlttax", clmxl//"T Zonal Advection", "C/s",       & ! i-adv.
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlttay", clmxl//"T Meridional Advection", "C/s",   & ! j-adv.
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlttaz", clmxl//"T Vertical Advection", "C/s",   & ! vert. adv.
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlttdh", clmxl//"T Horizontal Diffusion ", "C/s",   & ! hor. lateral diff.
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlttfo", clmxl//"T Forcing", "C/s",   & ! forcing
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlbtdz", clmxl//"T Vertical Diffusion", "C/s",   & ! vert. diff.
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlbtdt", clmxl//"T dh/dt Entrainment (Residual)", "C/s",   & ! T * dh/dt 
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         IF( l_traldf_iso ) THEN
-         CALL histdef( nidtrd, "somlbtdv", clmxl//"T Vert. lateral Diffusion","C/s",   & ! vertical diffusion entrainment (ISO)
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         ENDIF
-#if defined key_traldf_eiv
-         CALL histdef( nidtrd, "somlgtax", clmxl//"T Zonal EIV Advection", "C/s",   & ! i-adv. (eiv)
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlgtay", clmxl//"T Meridional EIV Advection", "C/s",   & ! j-adv. (eiv)
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlgtaz", clmxl//"T Vertical EIV Advection", "C/s",   & ! vert. adv. (eiv)
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlgtat", clmxl//"T Total EIV Advection", "C/s",   & ! total advection (eiv)
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-#endif
-         ! Salinity
-         CALL histdef( nidtrd, "somlsalt", clmxl//"Salinity", "PSU",   & ! Mixed-layer salinity
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         ! Salinity trends
-         CALL histdef( nidtrd, "somltsto", clmxl//"S Total", "PSU/s",   & ! total 
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somltsax", clmxl//"S Zonal Advection", "PSU/s",   & ! i-advection
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somltsay", clmxl//"S Meridional Advection", "PSU/s",   & ! j-advection
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somltsaz", clmxl//"S Vertical Advection", "PSU/s",   & ! vertical advection
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somltsdh", clmxl//"S Horizontal Diffusion ", "PSU/s",   & ! hor. lat. diff.
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somltsfo", clmxl//"S Forcing", "PSU/s",   & ! forcing
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-
-         CALL histdef( nidtrd, "somlbsdz", clmxl//"S Vertical Diffusion", "PSU/s",   & ! vert. diff.
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlbsdt", clmxl//"S dh/dt Entrainment (Residual)", "PSU/s",   & ! S * dh/dt 
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         IF( l_traldf_iso ) THEN
-         ! vertical diffusion entrainment (ISO)
-         CALL histdef( nidtrd, "somlbsdv", clmxl//"S Vertical lateral Diffusion", "PSU/s",   & ! vert. lat. diff.
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         ENDIF
-#if defined key_traldf_eiv
-         CALL histdef( nidtrd, "somlgsax", clmxl//"S Zonal EIV Advection", "PSU/s",   & ! i-advection (eiv)
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlgsay", clmxl//"S Meridional EIV Advection", "PSU/s",   & ! j-advection (eiv)
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlgsaz", clmxl//"S Vertical EIV Advection", "PSU/s",   & ! vert. adv. (eiv)
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nidtrd, "somlgsat", clmxl//"S Total EIV Advection", "PSU/s",   & ! total adv. (eiv)
-            &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-#endif
-         CALL histend( nidtrd )
-
-      ENDIF
Index: trunk/NEMO/OPA_SRC/TRD/trdmld_ncwrite.h90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdmld_ncwrite.h90	(revision 129)
+++ 	(revision )
@@ -1,54 +1,0 @@
-      ! define time axis
-      it= kt-nit000+1
-      IF( lwp .AND. MOD( kt, nwrite ) == 0 ) THEN
-         WRITE(numout,*) '     diamld : write NetCDF fields'
-      ENDIF
-      
-      CALL histwrite( nidtrd,"somlttml",it,rmld          ,ndimtrd1,ndextrd1) ! Mixed-layer depth
-      CALL histwrite( nidtrd,"somltemp",it,tml           ,ndimtrd1,ndextrd1) ! Mixed-layer temperature
-      
-      
-      ! Temperature trends
-      ! ------------------
-      CALL histwrite( nidtrd,"somlttto",it,ztmltot       ,ndimtrd1,ndextrd1) ! total 
-      CALL histwrite( nidtrd,"somlttax",it,tmltrd(:,:, 1),ndimtrd1,ndextrd1) ! i- adv.
-      CALL histwrite( nidtrd,"somlttay",it,tmltrd(:,:, 2),ndimtrd1,ndextrd1) ! j- adv.
-      CALL histwrite( nidtrd,"somlttaz",it,tmltrd(:,:, 3),ndimtrd1,ndextrd1) ! vertical adv.
-      CALL histwrite( nidtrd,"somlttdh",it,tmltrd(:,:, 4),ndimtrd1,ndextrd1) ! hor. lateral diff.
-      CALL histwrite( nidtrd,"somlttfo",it,tmltrd(:,:, 5),ndimtrd1,ndextrd1) ! forcing
-
-      CALL histwrite( nidtrd,"somlbtdz",it,tmltrd(:,:, 6),ndimtrd1,ndextrd1) ! vert. diffusion 
-      CALL histwrite( nidtrd,"somlbtdt",it,ztmlres       ,ndimtrd1,ndextrd1) ! dh/dt entrainment (residual)
-      IF( l_traldf_iso ) THEN
-         CALL histwrite( nidtrd,"somlbtdv",it,tmltrd(:,:, 7),ndimtrd1,ndextrd1) ! vert. lateral diff.
-      ENDIF
-#if defined key_traldf_eiv
-      CALL histwrite( nidtrd,"somlgtax",it,tmltrd(:,:, 8),ndimtrd1,ndextrd1) ! i- adv. (eiv)
-      CALL histwrite( nidtrd,"somlgtay",it,tmltrd(:,:, 9),ndimtrd1,ndextrd1) ! j- adv. (eiv)
-      CALL histwrite( nidtrd,"somlgtaz",it,tmltrd(:,:,10),ndimtrd1,ndextrd1) ! vert. adv. (eiv)
-      z2d(:,:) = tmltrd(:,:,8) + tmltrd(:,:,9) + tmltrd(:,:,10)
-      CALL histwrite( nidtrd,"somlgtat",it,z2d           ,ndimtrd1,ndextrd1) ! total adv. (eiv)
-#endif
-      CALL histwrite( nidtrd,"somlsalt",it,sml           ,ndimtrd1,ndextrd1) ! Mixed-layer salinity
-
-      ! Salinity trends
-      ! ---------------
-      CALL histwrite( nidtrd,"somltsto",it,zsmltot       ,ndimtrd1,ndextrd1) ! total 
-      CALL histwrite( nidtrd,"somltsax",it,smltrd(:,:, 1),ndimtrd1,ndextrd1) ! i- adv.
-      CALL histwrite( nidtrd,"somltsay",it,smltrd(:,:, 2),ndimtrd1,ndextrd1) ! j- adv.
-      CALL histwrite( nidtrd,"somltsaz",it,smltrd(:,:, 3),ndimtrd1,ndextrd1) ! vert. adv.
-      CALL histwrite( nidtrd,"somltsdh",it,smltrd(:,:, 4),ndimtrd1,ndextrd1) ! hor. lateral diff.
-      CALL histwrite( nidtrd,"somltsfo",it,smltrd(:,:, 5),ndimtrd1,ndextrd1) ! forcing
-      CALL histwrite( nidtrd,"somlbsdz",it,smltrd(:,:, 6),ndimtrd1,ndextrd1) ! vert. diff.
-      CALL histwrite( nidtrd,"somlbsdt",it,zsmlres       ,ndimtrd1,ndextrd1) ! dh/dt entrainment (residual)
-      IF( l_traldf_iso ) THEN
-         CALL histwrite( nidtrd,"somlbsdv",it,smltrd(:,:, 7),ndimtrd1,ndextrd1) ! vert. lateral diff;
-      ENDIF
-#if defined key_traldf_eiv
-      CALL histwrite( nidtrd,"somlgsax",it,smltrd(:,:, 8),ndimtrd1,ndextrd1) ! i-adv. (eiv)
-      CALL histwrite( nidtrd,"somlgsay",it,smltrd(:,:, 9),ndimtrd1,ndextrd1) ! j-adv. (eiv)
-      CALL histwrite( nidtrd,"somlgsaz",it,smltrd(:,:,10),ndimtrd1,ndextrd1) ! vert. adv. (eiv)
-      z2d(:,:) = smltrd(:,:,8) + smltrd(:,:,9) + smltrd(:,:,10)
-      CALL histwrite( nidtrd,"somlgsat",it,z2d           ,ndimtrd1,ndextrd1) ! total adv. (eiv)
-#endif
-
Index: trunk/NEMO/OPA_SRC/TRD/trdmld_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdmld_oce.F90	(revision 215)
+++ trunk/NEMO/OPA_SRC/TRD/trdmld_oce.F90	(revision 215)
@@ -0,0 +1,44 @@
+MODULE trdmld_oce
+   !!======================================================================
+   !!                   ***  MODULE trdmld_oce  ***
+   !! Ocean trends :   set tracer and momentum trend variables
+   !!======================================================================
+
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE par_oce         ! ocean parameters
+
+   IMPLICIT NONE
+   PUBLIC
+
+   INTEGER, PARAMETER ::            &  !: mixed layer trends index
+      jpmldxad = 1,   &  !: zonal advection
+      jpmldyad = 2,   &  !: meridionnal advection
+      jpmldzad = 3,   &  !: vertical advection
+      jpmldldf = 4,   &  !: lateral diffusion (horiz. component+Beckman)
+      jpmldfor = 5,   &  !: forcing 
+      jpmldevd = 6,   &  !: entrainment due to vertical diffusion (TKE)
+      jpmldzdf = 7,   &  !: explicit vertical part if isopycnal diffusion
+      jpmldxei = 8,   &  !: eddy induced zonal advection
+      jpmldyei = 9,   &  !: eddy induced meridional advection
+      jpmldzei =10       !: eddy induced vertical advection
+
+#if   defined  key_trdmld   ||   defined key_esopa
+   !!----------------------------------------------------------------------
+   !!   'key_trdmld'                         mixed layer trends diagnostics
+   !!----------------------------------------------------------------------
+
+   !! Trends diagnostics parameters
+   !!---------------------------------------------------------------------
+   INTEGER, PARAMETER ::            &  !:
+# if defined key_traldf_eiv
+      jpltrd = 10,  &  !: number of mixed-layer trends arrays
+      jpktrd = jpk     !: max level for mixed-layer trends diag.
+# else
+      jpltrd = 7,   &  !: number of mixed-layer trends arrays
+      jpktrd = jpk     !: max level for mixed-layer trends diag.
+# endif
+
+#endif
+  !!======================================================================
+END MODULE trdmld_oce
Index: trunk/NEMO/OPA_SRC/TRD/trdmod.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 215)
+++ trunk/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 215)
@@ -0,0 +1,439 @@
+MODULE trdmod
+   !!======================================================================
+   !!                       ***  MODULE  trdmod  ***
+   !! Ocean diagnostics:  ocean tracers and dynamic trends
+   !!=====================================================================
+#if  defined key_trdtra || defined key_trddyn || defined key_trdmld || defined key_trdvor || defined key_esopa
+   !!----------------------------------------------------------------------
+   !!   trd_mod          : Call the trend to be computed
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce                     ! ocean dynamics and tracers variables
+   USE dom_oce                 ! ocean space and time domain variables
+   USE trdmod_oce              ! ocean variables trends
+   USE trdvor                  ! ocean vorticity trends 
+   USE trdicp                  ! ocean bassin integral constraints properties
+   USE trdmld                  ! ocean active mixed layer tracers trends 
+   USE trabbl                  ! bottom boundary layer variables
+   USE in_out_manager          ! I/O manager
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC trd_mod        ! called by all dynXX or traXX modules
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , LODYC-IPSL  (2004)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trd_mod(ptrdx, ptrdy, ktrd, ctype, kt)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mod  ***
+      !! 
+      !! ** Purpose : Dispatch all trends computation, e.g. vorticity, mld or 
+      !!              integral constrains
+      !!
+      !! ** Method :
+      !!
+      !! History :
+      !!   9.0  !  04-08  (C. Talandier) New trends organization
+      !!----------------------------------------------------------------------
+      !! * Modules used
+#if defined key_trabbl_adv
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &  ! temporary arrays
+         &         zun, zvn
+#else
+      USE oce                , zun => un,  &  ! When no bbl, zun == un
+         &                     zvn => vn      ! When no bbl, zvn == vn
+#endif
+
+      !! * Arguments
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   &
+         ptrdx,                      &   ! Temperature or U trend 
+         ptrdy                           ! Salinity    or V trend
+
+      INTEGER, INTENT( in ) ::   &
+         kt  ,                   & ! time step
+         ktrd                      ! tracer trend index
+
+      CHARACTER(len=3), INTENT( in ) ::   &
+         ctype                             ! momentum or tracers trends type
+         !                                 ! 'DYN' or 'TRA'
+
+      !! * Local save
+      REAL(wp), DIMENSION(jpi,jpj), SAVE ::   &
+         zbtr2
+
+      !! * Local declarations
+      INTEGER ::   ji, jj, jk    ! loop indices
+      REAL(wp) ::   &
+         zbtr,            &  ! temporary scalars
+         zfui, zfvj,           &  !    "         "
+         zfui1, zfvj1             !    "         "
+      REAL(wp), DIMENSION(jpi,jpj) ::   &
+         z2dx, z2dy                        ! workspace arrays
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
+         z3dx, z3dy                            ! workspace arrays
+      !!----------------------------------------------------------------------
+
+      ! Initialization of workspace arrays
+      z3dx(:,:,:) = 0.e0
+      z3dy(:,:,:) = 0.e0
+      z2dx(:,:) = 0.e0
+      z2dy(:,:) = 0.e0
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! I. Bassin averaged properties for momentum and/or tracers trends
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      IF( ( mod(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend) )   THEN
+
+         ! Active tracers trends 
+         IF( lk_trdtra .AND. ctype == 'TRA' )   THEN
+
+            IF( ktrd == jpttdnsr )   THEN
+               ! 2D array tracers surface forcing
+               z2dx(:,:) = ptrdx(:,:,1)
+               z2dy(:,:) = ptrdy(:,:,1)
+
+               CALL trd(z2dx, z2dy, ktrd, ctype)
+            ELSE
+               ! 3D array
+               CALL trd(ptrdx, ptrdy, ktrd, ctype)
+            ENDIF
+
+         ENDIF
+
+         ! Momentum trends 
+         IF( lk_trddyn .AND. ctype == 'DYN' )   THEN
+
+            IF( ktrd == jpdtdswf .OR. ktrd == jpdtdbfr )   THEN
+               ! momentum surface forcing/bottom friction  2D array
+               z2dx(:,:) = ptrdx(:,:,1)
+               z2dy(:,:) = ptrdy(:,:,1)
+
+               CALL trd(z2dx, z2dy, ktrd, ctype)
+            ELSE
+               ! 3D array
+               CALL trd(ptrdx, ptrdy, ktrd, ctype)
+            ENDIF
+
+         ENDIF
+
+      ENDIF
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! II. Vorticity trends
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      IF( lk_trdvor .AND. ctype == 'DYN' )   THEN
+
+         SELECT CASE ( ktrd )
+
+         ! Pressure Gradient trend
+         CASE ( jpdtdhpg )      
+            CALL trd_vor_zint(ptrdx, ptrdy, jpvorprg)
+
+         ! KE Gradient trend
+         CASE ( jpdtdkeg )      
+            CALL trd_vor_zint(ptrdx, ptrdy, jpvorkeg)
+
+         ! Relative Vorticity trend
+         CASE ( jpdtdrvo )      
+            CALL trd_vor_zint(ptrdx, ptrdy, jpvorrvo)
+
+         ! Planetary Vorticity Term trend
+         CASE ( jpdtdpvo )      
+            CALL trd_vor_zint(ptrdx, ptrdy, jpvorpvo)
+
+         ! Horizontal Diffusion trend
+         CASE ( jpdtdldf )      
+            CALL trd_vor_zint(ptrdx, ptrdy, jpvorldf)
+
+         ! Vertical Advection trend
+         CASE ( jpdtdzad )      
+            CALL trd_vor_zint(ptrdx, ptrdy, jpvorzad)
+
+         ! Vertical Diffusion trend
+         CASE ( jpdtdzdf )      
+            CALL trd_vor_zint(ptrdx, ptrdy, jpvorzdf)
+
+         ! Surface Pressure Grad. trend
+         CASE ( jpdtdspg )      
+            CALL trd_vor_zint(ptrdx, ptrdy, jpvorspg)
+
+         ! Beta V trend 
+         CASE ( jpdtddat )      
+            CALL trd_vor_zint(ptrdx, ptrdy, jpvorbev)
+
+         ! Wind stress forcing term
+         CASE ( jpdtdswf )      
+            z2dx(:,:) = ptrdx(:,:,1)
+            z2dy(:,:) = ptrdy(:,:,1)
+
+            CALL trd_vor_zint(z2dx, z2dy, jpvorswf)
+
+         ! Bottom friction term
+         CASE ( jpdtdbfr )      
+            z2dx(:,:) = ptrdx(:,:,1)
+            z2dy(:,:) = ptrdy(:,:,1)
+
+            CALL trd_vor_zint(z2dx, z2dy, jpvorbfr)
+
+         END SELECT
+
+      ENDIF
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! III. Mixed layer trends
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      IF( lk_trdmld .AND. ctype == 'TRA' )   THEN
+         
+         SELECT CASE ( ktrd )
+
+         ! horizontal advection trends
+         CASE ( jpttdlad )      
+
+#if defined key_trabbl_adv
+            ! Advective bottom boundary layer 
+            ! -------------------------------
+            zun(:,:,:) = un(:,:,:) - u_bbl(:,:,:)
+            zvn(:,:,:) = vn(:,:,:) - v_bbl(:,:,:)
+#endif
+            IF( kt == nit000 )   zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
+
+            SELECT CASE ( l_adv )
+
+            CASE ( 'ce2' )
+
+               ! Split horizontal trends into i- and j- compnents for trdmld case 
+               ! ----------------------------------------------------------------
+
+               ! i- advective trend computed as Uh gradh(T)
+               DO jk = 1, jpkm1
+                  DO jj = 2, jpjm1
+                     DO ji = fs_2, fs_jpim1   ! vector opt.
+# if defined key_s_coord || defined key_partial_steps
+                        zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
+
+                        zfui = 0.5 * e2u(ji  ,jj) * fse3u(ji,  jj,jk) * zun(ji,  jj,jk)
+                        zfui1= 0.5 * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zun(ji-1,jj,jk)
+# else         
+                        zbtr = zbtr2(ji,jj)
+
+                        zfui = 0.5 * e2u(ji  ,jj) * zun(ji,  jj,jk)
+                        zfui1= 0.5 * e2u(ji-1,jj) * zun(ji-1,jj,jk)
+# endif
+                        ! save i- advective trend 
+                        z3dx(ji,jj,jk) = - zbtr * ( zfui  * ( tn(ji+1,jj,jk) - tn(ji  ,jj,jk) )    &
+                            &                     + zfui1 * ( tn(ji  ,jj,jk) - tn(ji-1,jj,jk) ) )
+                        z3dy(ji,jj,jk) = - zbtr * ( zfui  * ( sn(ji+1,jj,jk) - sn(ji  ,jj,jk) )    &
+                            &                     + zfui1 * ( sn(ji  ,jj,jk) - sn(ji-1,jj,jk) ) )
+                     END DO
+                  END DO
+               END DO
+
+               ! save the i- horizontal trends for diagnostic
+               CALL trd_mld_zint(z3dx, z3dy, jpmldxad, '3D')
+
+               ! j- advective trend computed as Uh gradh(T)
+               DO jk = 1, jpkm1
+                  DO jj = 2, jpjm1
+                     DO ji = fs_2, fs_jpim1   ! vector opt.
+# if defined key_s_coord || defined key_partial_steps
+                        zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
+
+                        zfvj = 0.5 * e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * zvn(ji,jj  ,jk)
+                        zfvj1= 0.5 * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * zvn(ji,jj-1,jk)
+# else         
+                        zbtr = zbtr2(ji,jj)
+
+                        zfvj = 0.5 * e1v(ji,jj  ) * zvn(ji,jj  ,jk)
+                        zfvj1= 0.5 * e1v(ji,jj-1) * zvn(ji,jj-1,jk)
+# endif
+                        ! save j- advective trend 
+                        z3dx(ji,jj,jk) = - zbtr * ( zfvj  * ( tn(ji,jj+1,jk) - tn(ji,jj  ,jk) )   &
+                            &                     + zfvj1 * ( tn(ji,jj  ,jk) - tn(ji,jj-1,jk) ) )
+                        z3dy(ji,jj,jk) = - zbtr * ( zfvj  * ( sn(ji,jj+1,jk) - sn(ji,jj  ,jk) )   &
+                            &                     + zfvj1 * ( sn(ji,jj  ,jk) - sn(ji,jj-1,jk) ) )
+                     END DO
+                  END DO
+               END DO
+
+               ! save the j- horizontal trend for diagnostic
+               CALL trd_mld_zint(z3dx, z3dy, jpmldyad, '3D')
+
+            CASE ( 'tvd' )
+
+               ! Recompute the horizontal advection term Div(Uh.T) term 
+               z3dx(:,:,:) = ptrdx(:,:,:) - tn(:,:,:) * hdivn(:,:,:)
+               z3dy(:,:,:) = ptrdy(:,:,:) - sn(:,:,:) * hdivn(:,:,:)
+
+               ! Deduce the i- horizontal advection in substracting the j- one.
+               ! tladj()/sladj() are computed in traadv_tvd.F90 module
+               z3dx(:,:,:) = z3dx(:,:,:) - tladj(:,:,:)
+               z3dy(:,:,:) = z3dy(:,:,:) - sladj(:,:,:)
+
+               DO jk = 1, jpkm1
+                  DO jj = 2, jpjm1
+                     DO ji = fs_2, fs_jpim1
+                        zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
+
+                        ! Compute the zonal et meridional divergence
+                        zfui = e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * zun(ji  ,jj,jk)  &
+                             - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zun(ji-1,jj,jk)
+                        zfvj = e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * zvn(ji,jj  ,jk)  &
+                             - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * zvn(ji,jj-1,jk)
+
+                        ! i- advective trend computed as U gradx(T/S)
+                        z3dx(ji,jj,jk) = z3dx(ji,jj,jk) + tn(ji,jj,jk) * zfui * zbtr
+                        z3dy(ji,jj,jk) = z3dy(ji,jj,jk) + sn(ji,jj,jk) * zfui * zbtr
+
+                        ! j- advective trend computed as V grady(T/S)
+                        tladj(ji,jj,jk) = tladj(ji,jj,jk) + tn(ji,jj,jk) * zfvj * zbtr
+                        sladj(ji,jj,jk) = sladj(ji,jj,jk) + sn(ji,jj,jk) * zfvj * zbtr
+
+                     END DO
+                  END DO
+               END DO
+
+               ! save the i- horizontal trend for diagnostic
+               CALL trd_mld_zint(z3dx, z3dy, jpmldxad, '3D')
+
+               ! save the j- horizontal trend for diagnostic
+               CALL trd_mld_zint(tladj, sladi, jpmldyad, '3D')
+
+            CASE ( 'mus', 'mu2' )
+
+               !  Split horizontal trends in i- and j- direction for trdmld case 
+               ! ----------------------------------------------------------------
+
+               ! i- advective trend computed as U gradx(T/S)
+               DO jk = 1, jpkm1
+                  DO jj = 2, jpjm1      
+                     DO ji = fs_2, fs_jpim1   ! vector opt.
+# if defined key_s_coord || defined key_partial_steps
+                        zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
+                        zfui =  e2u(ji  ,jj) * fse3u(ji,  jj,jk) * zun(ji,  jj,jk)   &
+                           & -  e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zun(ji-1,jj,jk)
+# else      
+                        zbtr = zbtr2(ji,jj)
+                        zfui = e2u(ji  ,jj) * zun(ji,  jj,jk)   &
+                           & - e2u(ji-1,jj) * zun(ji-1,jj,jk)
+# endif
+                        ! save i- advective trend 
+                        z3dx(ji,jj,jk) = - zbtr * ( tladi(ji,jj,jk) - tladi(ji-1,jj,jk) )   &
+                            &                      + tn(ji,jj,jk) * zfui * zbtr
+                        z3dy(ji,jj,jk) = - zbtr * ( sladi(ji,jj,jk) - sladi(ji-1,jj,jk) )  &
+                            &                      + sn(ji,jj,jk) * zfui * zbtr
+                     END DO
+                  END DO
+               END DO        
+
+               ! save the i- horizontal trends for diagnostic
+               CALL trd_mld_zint(z3dx, z3dy, jpmldxad, '3D')
+
+               ! j- advective trend computed as V grady(T/S)
+               DO jk = 1, jpkm1
+                  DO jj = 2, jpjm1      
+                     DO ji = fs_2, fs_jpim1   ! vector opt.
+# if defined key_s_coord || defined key_partial_steps
+                        zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
+                        zfvj =  e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * zvn(ji,jj  ,jk)   &
+                           & -  e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * zvn(ji,jj-1,jk)
+# else      
+                        zbtr = zbtr2(ji,jj)
+                        zfvj = e1v(ji,jj  ) * zvn(ji,jj  ,jk)   &
+                           & - e1v(ji,jj-1) * zvn(ji,jj-1,jk)
+# endif
+                        ! save j- advective trend 
+                        z3dx(ji,jj,jk) =  - zbtr * ( tladj(ji,jj,jk) - tladj(ji,jj-1,jk) )   &
+                            &                       + tn(ji,jj,jk) * zfvj * zbtr
+                        z3dy(ji,jj,jk) =  - zbtr * ( sladj(ji,jj,jk) - sladj(ji,jj-1,jk) )   &
+                            &                       + sn(ji,jj,jk) * zfvj * zbtr
+                     END DO
+                  END DO
+               END DO        
+
+               ! save the j- horizontal trends for diagnostic
+               CALL trd_mld_zint(z3dx, z3dy, jpmldyad, '3D')
+
+            END SELECT
+
+         ! vertical advection trends
+         CASE ( jpttdzad )      
+            CALL trd_mld_zint(ptrdx, ptrdy, jpmldzad, '3D')
+
+         ! lateral diffusion trends
+         CASE ( jpttdldf )      
+            CALL trd_mld_zint(ptrdx, ptrdy, jpmldldf, '3D')
+# if defined key_traldf_eiv
+            ! Save the i- and j- eddy induce velocity trends
+            CALL trd_mld_zint(tladi, sladi, jpmldxei, '3D')
+            CALL trd_mld_zint(tladj, sladj, jpmldyei, '3D')
+# endif
+            IF( lk_trabbl_dif )   THEN
+               z3dx(:,:,:) = 0.e0
+               z3dy(:,:,:) = 0.e0
+               z3dx(:,:,1) = tldfbbl(:,:)
+               z3dy(:,:,1) = sldfbbl(:,:)
+               CALL trd_mld_zint(z3dx, z3dy, jpmldldf, '2D')
+            ENDIF
+
+         ! vertical diffusion trends
+         CASE ( jpttdzdf )      
+            CALL trd_mld_zint(ptrdx, ptrdy, jpmldzdf, '3D')
+
+         ! vertical diffusion trends
+         CASE ( jpttddoe )      
+            CALL trd_mld_zint(ptrdx, ptrdy, jpmldzei, '3D')
+
+         ! penetrative solar radiation trends
+         CASE ( jpttdqsr )      
+            CALL trd_mld_zint(ptrdx, ptrdy, jpmldfor, '3D')
+
+         ! non penetrative solar radiation trends
+         CASE ( jpttdnsr )
+            ptrdx(:,:,2:jpk) = 0.e0
+            ptrdy(:,:,2:jpk) = 0.e0
+            CALL trd_mld_zint(ptrdx, ptrdy, jpmldfor, '2D')
+
+         END SELECT   
+
+      ENDIF
+
+
+   END SUBROUTINE trd_mod
+
+#   else
+   !!----------------------------------------------------------------------
+   !!   Default case :                                         Empty module
+   !!----------------------------------------------------------------------
+   USE trdmod_oce      ! ocean variables trends
+
+CONTAINS
+   SUBROUTINE trd_mod(ptrd3dx, ptrd3dy, ktrd , ctype, kt)       ! Empty routine
+      REAL, DIMENSION(:,:,:), INTENT( in ) ::   &
+          ptrd3dx,                     &   ! Temperature or U trend 
+          ptrd3dy                          ! Salinity    or V trend
+      INTEGER, INTENT( in ) ::   ktrd      ! momentum or tracer trend index
+      INTEGER, INTENT( in ) ::   kt        ! Time step
+      CHARACTER(len=3), INTENT( in ) ::   &
+         ctype                             ! momentum or tracers trends type
+      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', ptrd3dx(1,1,1)
+      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ptrd3dy(1,1,1)
+      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ktrd
+      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ctype
+      WRITE(*,*) ' "   ": You should not have seen this print! error ?', kt
+   END SUBROUTINE trd_mod
+#   endif
+
+   !!======================================================================
+END MODULE trdmod
Index: trunk/NEMO/OPA_SRC/TRD/trdmod_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdmod_oce.F90	(revision 215)
+++ trunk/NEMO/OPA_SRC/TRD/trdmod_oce.F90	(revision 215)
@@ -0,0 +1,17 @@
+MODULE trdmod_oce
+   !!======================================================================
+   !!                   ***  MODULE trdmod_oce  ***
+   !! Ocean trends :   set tracer and momentum trend variables
+   !!======================================================================
+   !! * Modules used
+   USE trdicp_oce              ! ocean momentum/tracers bassin properties trends variables
+   USE trdmld_oce              ! ocean active mixed layer tracers trends variables
+   USE trdvor_oce              ! ocean vorticity trends variables
+
+   !! Control parameters
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC ::   l_trdtra = .FALSE.    !: tracers  trend flag
+   LOGICAL, PUBLIC ::   l_trddyn = .FALSE.    !: momentum trend flag
+
+  !!======================================================================
+END MODULE trdmod_oce
Index: trunk/NEMO/OPA_SRC/TRD/trdtra.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 129)
+++ 	(revision )
@@ -1,379 +1,0 @@
-MODULE trdtra
-   !!======================================================================
-   !!                       ***  MODULE  trdtra  ***
-   !! Ocean diagnostics:  ocean tracer trends
-   !!=====================================================================
-#if  defined key_trdtra   ||   defined key_esopa
-   !!----------------------------------------------------------------------
-   !!   'key_trdtra'  or                           tracer trend diagnostics
-   !!   'key_trdmld'                          mixed layer trend diagnostics
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   trd_tra      : verify the basin averaged properties for tracers
-   !!   trd_tra_init : ???
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables
-   USE trdtra_oce      ! ocean active tracer trend variables
-   USE trddyn_oce      ! ocean dynamics trend variables
-   USE ldftra_oce      ! ocean active tracers: lateral physics
-   USE ldfdyn_oce      ! ocean dynamics: lateral physics
-   USE zdf_oce         ! ocean vertical physics
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distibuted memory computing library
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! * Routine accessibility
-   PUBLIC trd_tra      ! called by step.F90
-   PUBLIC trd_tra_init ! called by opa.F90
-
-   !! * Shared module variables
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trdtra = .TRUE.    !: momentum trend flag
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LODYC-IPSL  (2003)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE trd_tra( kt )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_tra  ***
-      !! 
-      !! ** Purpose :   verify the basin averaged properties of the tracers
-      !!   equations at every time step frequency ntrd.
-      !!
-      !!   Method :
-      !!
-      !! History :
-      !!        !  91-12 (G. Madec)
-      !!        !  92-06 (M. Imbard) add time step frequency
-      !!        !  96-01 (G. Madec) terrain following coordinates
-      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-
-      !! * Local declarations
-      INTEGER :: ji, jj, jk, jn
-      REAL(wp)                ::   zbt 
-      REAL(wp), DIMENSION(10) ::   ztmo, zsmo, zt2, zs2
-
-      NAMELIST/namtrd/ ntrd, nctls
-      !!----------------------------------------------------------------------
-
-      ! 0. Initialization
-      ! -----------------
-
-      IF( kt == nit000 ) THEN
-
-         ! namelist namtrd : trend diagnostic
-         REWIND( numnam )
-         READ  ( numnam, namtrd )
-
-         IF(lwp) THEN
-            WRITE(numout,*) 'trd_tra : read namelist namtrd'
-            WRITE(numout,*) '~~~~~~~'
-            WRITE(numout,*) ' time step frequency trend       ntrd  = ', ntrd
-            WRITE(numout,*) ' '
-         ENDIF
-
-         ! Total volume at t-points:
-         tvolt = 0.
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  tvolt = tvolt + e1t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk)* tmask(ji,jj,jk) * tmask_i(ji,jj)
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
-
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) '***trd_tra'
-            WRITE(numout,*) ' frequency ntrd = ',ntrd
-            WRITE(numout,*) ' '
-            WRITE(numout,*) ' tvolt = ',tvolt
-            WRITE(numout,*) ' '
-         ENDIF
-
-      ENDIF
-
-
-      ! 1. Advective trends and forcing trend
-      ! -------------------------------------
-
-      IF( MOD(kt,ntrd) == 0 .OR. kt  == nit000 .OR. kt == nitend ) THEN
-         
-         ! 1.1 Mask the forcing trend and substract it from the vertical diffusion trend
-         flxtrd(:,:,1) = flxtrd(:,:,1) * tmask_i(:,:)
-         flxtrd(:,:,2) = flxtrd(:,:,2) * tmask_i(:,:)
-         
-         ! 1.2 Mask the trends
-         DO jn = 1, 6
-            DO jk = 1, jpk
-               ttrd(:,:,jk,jn) = ttrd(:,:,jk,jn)* tmask(:,:,jk) * tmask_i(:,:)
-               strd(:,:,jk,jn) = strd(:,:,jk,jn)* tmask(:,:,jk) * tmask_i(:,:)
-            END DO
-         END DO
-         
-         DO jk = 1, jpk
-            ttrd(:,:,jk,7) = ttrd(:,:,jk,7) * tmask(:,:,jk) * tmask_i(:,:)
-         END DO
-
-
-         ! 2. Basin averaged tracer trend
-         ! ------------------------------
-         
-         DO jn = 1, 6
-            ztmo(jn) = 0.
-            zsmo(jn) = 0.
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) 
-                     ztmo(jn) =  ztmo(jn) + ttrd(ji,jj,jk,jn) * zbt
-                     zsmo(jn) =  zsmo(jn) + strd(ji,jj,jk,jn) * zbt
-                  END DO
-               END DO
-            END DO
-         END DO
-         
-         ztmo(7) = 0.
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-                  ztmo(7) =  ztmo(7) + ttrd(ji,jj,jk,7) * zbt
-               END DO
-            END DO
-         END DO
-         
-         ztmo(8) = 0.
-         zsmo(8) = 0.
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
-               ztmo(8) =  ztmo(8) + flxtrd(ji,jj,1) * zbt
-               zsmo(8) =  zsmo(8) + flxtrd(ji,jj,2) * zbt
-            END DO
-         END DO
-         
-         
-         ! 3. Basin averaged tracer square trend
-         ! -------------------------------------
-         ! c a u t i o n: field before, because after the array swap
-         
-         DO jn = 1, 6
-            zt2(jn) = 0.e0
-            zs2(jn) = 0.e0
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-                     zt2(jn) = zt2(jn) + ttrd(ji,jj,jk,jn) * zbt * tb(ji,jj,jk)
-                     zs2(jn) = zs2(jn) + strd(ji,jj,jk,jn) * zbt * sb(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-         END DO
-         
-         zt2(7) = 0.e0
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-                  zt2(7) = zt2(7) + ttrd(ji,jj,jk,7) * zbt * tb(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         
-         zt2(8) = 0.e0
-         zs2(8) = 0.e0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
-               zt2(8) = zt2(8) + flxtrd(ji,jj,1) * zbt * tb(ji,jj,1)
-               zs2(8) = zs2(8) + flxtrd(ji,jj,2) * zbt * sb(ji,jj,1)
-            END DO
-         END DO
-         
-         IF( lk_mpp ) THEN
-            CALL mpp_sum( ztmo, 10 )   ! sums over the global domain
-            CALL mpp_sum( zsmo, 10 )
-            CALL mpp_sum( zt2 , 10 )
-            CALL mpp_sum( zs2 , 10 )
-         ENDIF
-         
-         ! 4. Print
-         ! --------
-         
-         IF(lwp) THEN
-            WRITE (numout,*)
-            WRITE (numout,*)
-            WRITE (numout,9400) kt
-            WRITE (numout,9401) ztmo(1) / tvolt, zsmo(1) / tvolt
-            WRITE (numout,9402) ztmo(2) / tvolt, zsmo(2) / tvolt
-            WRITE (numout,9403) ztmo(3) / tvolt, zsmo(3) / tvolt
-            WRITE (numout,9404) ztmo(4) / tvolt, zsmo(4) / tvolt
-            WRITE (numout,9405) ztmo(5) / tvolt, zsmo(5) / tvolt
-            WRITE (numout,9406) ztmo(6) / tvolt, zsmo(6) / tvolt
-            WRITE (numout,9407) ztmo(7) / tvolt
-            WRITE (numout,9408) ztmo(8) / tvolt, zsmo(8) / tvolt
-            WRITE (numout,9409)
-            WRITE (numout,9410) (  ztmo(1) + ztmo(2) + ztmo(3) + ztmo(4)              &
-            &                    + ztmo(5) + ztmo(6) + ztmo(7) + ztmo(8) ) / tvolt,   &
-            &                   (  zsmo(1) + zsmo(2) + zsmo(3) + zsmo(4)              &
-            &                    + zsmo(5) + zsmo(6)           + zsmo(8) ) / tvolt
-         ENDIF
-
-9400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
-              '              salinity',/' ============================')
-9401     FORMAT(' horizontal advection        ',e20.13,'     ',e20.13)
-9402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
-9403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
-9404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
-9405     FORMAT(' STATIC instability mixing   ',e20.13,'     ',e20.13)
-9406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
-9407     FORMAT(' penetrative qsr             ',e20.13,'     ',e20.13)
-9408     FORMAT(' forcing term                ',e20.13,'     ',e20.13)
-9409     FORMAT(' -------------------------------------------------------------------------')
-9410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
-
-
-         IF(lwp) THEN
-            WRITE (numout,*)
-            WRITE (numout,*)
-            WRITE (numout,9420) kt
-            WRITE (numout,9421) zt2(1) / tvolt, zs2(1) / tvolt
-            WRITE (numout,9422) zt2(2) / tvolt, zs2(2) / tvolt
-            WRITE (numout,9423) zt2(3) / tvolt, zs2(3) / tvolt
-            WRITE (numout,9424) zt2(4) / tvolt, zs2(4) / tvolt
-            WRITE (numout,9425) zt2(5) / tvolt, zs2(5) / tvolt
-            WRITE (numout,9426) zt2(6) / tvolt, zs2(6) / tvolt
-            WRITE (numout,9427) zt2(7) / tvolt
-            WRITE (numout,9428) zt2(8) / tvolt, zs2(8) / tvolt
-            WRITE (numout,9429)
-            WRITE (numout,9430) (  zt2(1) + zt2(2) + zt2(3) + zt2(4)              &
-            &                    + zt2(5) + zt2(6) + zt2(7) + zt2(8) ) / tvolt,   &
-            &                   (  zs2(1) + zs2(2) + zs2(3) + zs2(4)              &
-            &                    + zs2(5) + zs2(6)          + zs2(8) ) / tvolt
-         ENDIF
-
-9420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
-            '               salinity', /, ' ===============================')
-9421     FORMAT(' horizontal advection      * t   ', e20.13, '     ', e20.13)
-9422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
-9423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
-9424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
-9425     FORMAT(' STATIC instability mixing * t   ', e20.13, '     ', e20.13)
-9426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
-9427     FORMAT(' penetrative qsr           * t   ', e20.13, '     ', e20.13)
-9428     FORMAT(' forcing term              * t   ', e20.13, '     ', e20.13)
-9429     FORMAT(' -----------------------------------------------------------------------------')
-9430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
-
-
-         IF(lwp) THEN
-            WRITE (numout,*)
-            WRITE (numout,*)
-            WRITE (numout,9440) kt
-            WRITE (numout,9441) ( ztmo(1)+ztmo(2) )/tvolt, ( zsmo(1)+zsmo(2) )/tvolt
-            WRITE (numout,9442)   ztmo(3)/tvolt,  zsmo(3)/tvolt
-            WRITE (numout,9443)   ztmo(4)/tvolt,  zsmo(4)/tvolt
-            WRITE (numout,9444)   ztmo(5)/tvolt,  zsmo(5)/tvolt
-            WRITE (numout,9445) ( zt2(1)+zt2(2) )/tvolt, ( zs2(1)+zs2(2) )/tvolt
-            WRITE (numout,9446)   zt2(3)/tvolt,   zs2(3)/tvolt
-            WRITE (numout,9447)   zt2(4)/tvolt,   zs2(4)/tvolt
-            WRITE (numout,9448)   zt2(5)/tvolt,   zs2(5)/tvolt
-         ENDIF
-
-9440     FORMAT(' tracer consistency at it= ',i6,   &
-            ' :         temperature','                salinity',/,   &
-            ' ==================================')
-9441     FORMAT(' 0 = horizontal+vertical advection      ',e20.13,'       ',e20.13)
-9442     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
-9443     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
-9444     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
-9445     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
-9446     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
-9447     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
-9448     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
-      ENDIF
-
-   END SUBROUTINE trd_tra
-
-
-   SUBROUTINE trd_tra_init
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_tra_init  ***
-      !! 
-      !! ** Purpose :   
-      !!
-      !! ** Method  :
-      !!
-      !! History :
-      !!   9.0  !  03-09  (G. Madec)  Original code
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER :: ji, jj, jk
-
-      NAMELIST/namtrd/ ntrd, nctls
-      !!----------------------------------------------------------------------
-
-      ! set to zero the tracers trends
-      ttrd  (:,:,:,:) = 0.e0
-      strd  (:,:,:,:) = 0.e0
-      ttrdh (:,:,:,:) = 0.e0
-      strdh (:,:,:,:) = 0.e0
-      flxtrd(:,:,  :) = 0.e0
-
-      ! namelist namtrd : trend diagnostic
-      REWIND( numnam )
-      READ  ( numnam, namtrd )
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'trd_tra : read namelist namtrd'
-         WRITE(numout,*) '~~~~~~~'
-         WRITE(numout,*) '          time step frequency trend       ntrd  = ', ntrd
-      ENDIF
-
-      ! Total volume at t-points:
-      tvolt = 0.
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               tvolt = tvolt + e1t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk)* tmask(ji,jj,jk) * tmask_i(ji,jj)
-            END DO
-         END DO
-      END DO
-      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
-
-      IF(lwp) WRITE(numout,*) '          total ocean volume at T-point   tvolt = ',tvolt
-
-   END SUBROUTINE trd_tra_init
-
-#   else
-   !!----------------------------------------------------------------------
-   !!   Default case :                                         Empty module
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trdtra = .FALSE.   !: momentum trend flag
-CONTAINS
-   SUBROUTINE trd_tra( kt )        ! Empty routine
-      WRITE(*,*) 'trd_tra: You should not have seen this print! error?', kt
-   END SUBROUTINE trd_tra
-   SUBROUTINE trd_tra_init         ! Empty routine
-   END SUBROUTINE trd_tra_init
-#endif
-
-   !!======================================================================
-END MODULE trdtra
Index: trunk/NEMO/OPA_SRC/TRD/trdtra_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdtra_oce.F90	(revision 129)
+++ 	(revision )
@@ -1,70 +1,0 @@
-MODULE trdtra_oce
-   !!======================================================================
-   !!                   ***  MODULE trdtra_oce  ***
-   !! Ocean trends :   set tracer trend variables
-   !!======================================================================
-
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE par_oce      ! ocean parameters
-
-   IMPLICIT NONE
-   PUBLIC
-
-   INTEGER  ::      & !!: namdia :  diagnostics on dynamics and/or tracer trends
-      ntrd  = 10 ,  &  !: time step frequency dynamics and tracers trends
-      nctls =  0       !: control surface type for trends vertical integration
-   !!----------------------------------------------------------------------
-
-#if   defined key_trdtra   ||   defined key_trdmld   ||   defined key_esopa
-   !!----------------------------------------------------------------------
-   !!   'key_trdtra'   or                         tracer trends diagnostics
-   !!   'key_trdmld'                         mixed layer trends diagnostics
-   !!----------------------------------------------------------------------
-
-   !! Trends diagnostics parameters
-   !!---------------------------------------------------------------------
-   INTEGER, PARAMETER ::            &  !:
-# if defined key_traldf_eiv
-      jptrdh = 4,   &  !: number of 3D horiz trends arrays
-      jpltrd = 10,  &  !: number of mixed-layer trends arrays
-      jpktrd = jpk     !: max level for mixed-layer trends diag.
-# else
-      jptrdh = 2,   &  !: number of 3D horiz trends arrays
-      jpltrd = 7,   &  !: number of mixed-layer trends arrays
-      jpktrd = jpk     !: max level for mixed-layer trends diag.
-# endif
-
-   !! Trends diagnostics variables
-   !!---------------------------------------------------------------------
-   REAL(wp) ::   &  !:
-      tvolt         !: volume of the whole ocean computed at t-points
-   REAL(wp), DIMENSION(jpi,jpj,jpk,7) ::   &  !:
-      ttrd             !: trends of the temperature tracer equations
-      !                !  ttrd(,,,1) : horizontal advection
-      !                !  ttrd(,,,2) : vertical advection
-      !                !  ttrd(,,,3) : horizontal diffusion
-      !                !  ttrd(,,,4) : vertical diffusion
-      !                !  ttrd(,,,5) : static instability
-      !                !  ttrd(,,,6) : damping OR vertical EIV
-      !                !  ttrd(,,,7) : penetrative solar radiation (T only)
-   REAL(wp), DIMENSION(jpi,jpj,jpk,6) ::   &  !:
-      strd             !: trends of the salinity tracer equations
-      !                !  same as ttrd()
-   REAL(wp), DIMENSION(jpi,jpj,jpk,jptrdh) ::   &  !:
-      ttrdh, strdh     !: ttrdh(,,,1) : zonal advection
-      !                !  ttrdh(,,,2) : meridional advection
-      !                !  ttrdh(,,,3) : zonal EIV
-      !                !  ttrdh(,,,4) : meridional EIV
-   REAL(wp), DIMENSION(jpi,jpj,2) ::   &  !:
-      flxtrd,       &  !: tracer forcing trends
-      bbltrd           !: tracer bottom boundary layer trends
-   
-#else
-  !!----------------------------------------------------------------------
-  !!   Default option :                        NO tracer trend diagnostics
-  !!----------------------------------------------------------------------
-#endif
-
-  !!======================================================================
-END MODULE trdtra_oce
Index: trunk/NEMO/OPA_SRC/TRD/trdvor_ncinit.h90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdvor_ncinit.h90	(revision 129)
+++ 	(revision )
@@ -1,61 +1,0 @@
-      !-----------------------------------------
-      ! II.1 Define frequency of output and means
-      ! -----------------------------------------
-#if defined key_diainstant
-      zsto = nwrite*rdt
-      clop ="inst(only(x))"
-#else
-      zsto = rdt
-      clop ="ave(x)"
-#endif
-      zout = ntrd*rdt
-
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE (numout,*) ' Diavor: netCDF initialization'
-
-         ! II.2 Compute julian date from starting date of the run
-         ! ------------------------
-         CALL ymds2ju( nyear, nmonth, nday, 0.e0, zjulian )
-         IF (lwp) WRITE(numout,*)' '  
-         IF (lwp) WRITE(numout,*)' Date 0 used :',nit000         &
-              ,' YEAR ', nyear,' MONTH ', nmonth,' DAY ', nday   &
-              ,'Julian day : ', zjulian
-
-         ! II.3 Define the T grid trend file (nidvor)
-         ! ---------------------------------
-         CALL dia_nam( clhstnam, ntrd, 'vort' )                  ! filename
-         IF(lwp) WRITE(numout,*) ' Name of NETCDF file ', clhstnam
-         CALL histbeg( clhstnam, jpi, glamf, jpj, gphif,1, jpi,   &  ! Horizontal grid : glamt and gphit
-            &          1, jpj, 0, zjulian, rdt, nh_t, nidvor)
-         CALL wheneq( jpi*jpj, fmask, 1, 1., ndexvor1, ndimvor1 )    ! surface
-
-         ! Declare output fields as netCDF variables
-         CALL histdef( nidvor, "sovortPh", cvort//"grad Ph" , "s-2",        & ! grad Ph
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "sovortEk", cvort//"Energy", "s-2",          & ! Energy
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "sovozeta", cvort//"rel vorticity", "s-2",   & ! rel vorticity
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "sovortif", cvort//"coriolis", "s-2",        & ! coriolis
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "sovodifl", cvort//"lat diff ", "s-2",       & ! lat diff
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "sovoadvv", cvort//"vert adv", "s-2",        & ! vert adv
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "sovodifv", cvort//"vert diff" , "s-2",      & ! vert diff
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "sovortPs", cvort//"grad Ps", "s-2",         & ! grad Ps
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "sovortbv", cvort//"Beta V", "s-2",          & ! beta.V
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "sovowind", cvort//"wind stress", "s-2",     & ! wind stress
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "sovobfri", cvort//"bottom friction", "s-2", & ! bottom friction
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "1st_mbre", cvort//"1st mbre", "s-2",        & ! First membre
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histdef( nidvor, "sovorgap", cvort//"gap", "s-2",             & ! gap between 1st and 2 nd mbre
-            &          jpi,jpj,nh_t,1,1,1,-99,32,clop,zsto,zout)
-         CALL histend( nidvor )
-
-      ENDIF
Index: trunk/NEMO/OPA_SRC/TRD/trdvor_ncwrite.h90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdvor_ncwrite.h90	(revision 129)
+++ 	(revision )
@@ -1,19 +1,0 @@
-     ! define time axis
-      it= kt-nit000+1
-      IF( lwp .AND. MOD( kt, ntrd ) == 0 ) THEN
-         WRITE(numout,*) '     diavor : write NetCDF fields'
-      ENDIF
- 
-      CALL histwrite( nidvor,"sovortPh",it,vortrd(:,:,1),ndimvor1,ndexvor1)  ! grad Ph
-      CALL histwrite( nidvor,"sovortEk",it,vortrd(:,:,2),ndimvor1,ndexvor1)  ! Energy
-      CALL histwrite( nidvor,"sovozeta",it,vortrd(:,:,3),ndimvor1,ndexvor1)  ! rel vorticity
-      CALL histwrite( nidvor,"sovortif",it,vortrd(:,:,4),ndimvor1,ndexvor1)  ! coriolis
-      CALL histwrite( nidvor,"sovodifl",it,vortrd(:,:,5),ndimvor1,ndexvor1)  ! lat diff
-      CALL histwrite( nidvor,"sovoadvv",it,vortrd(:,:,6),ndimvor1,ndexvor1)  ! vert adv
-      CALL histwrite( nidvor,"sovodifv",it,vortrd(:,:,7),ndimvor1,ndexvor1)  ! vert diff
-      CALL histwrite( nidvor,"sovortPs",it,vortrd(:,:,8),ndimvor1,ndexvor1)  ! grad Ps
-      CALL histwrite( nidvor,"sovortbv",it,vortrd(:,:,9),ndimvor1,ndexvor1)  ! beta.V
-      CALL histwrite( nidvor,"sovowind",it,vortrd(:,:,10),ndimvor1,ndexvor1) ! wind stress
-      CALL histwrite( nidvor,"sovobfri",it,vortrd(:,:,11),ndimvor1,ndexvor1) ! bottom friction
-      CALL histwrite( nidvor,"1st_mbre",it,vor_avrtot    ,ndimvor1,ndexvor1) ! First membre
-      CALL histwrite( nidvor,"sovorgap",it,vor_avrres    ,ndimvor1,ndexvor1) ! gap between 1st and 2 nd mbre
Index: trunk/NEMO/OPA_SRC/TRD/trdvor_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdvor_oce.F90	(revision 215)
+++ trunk/NEMO/OPA_SRC/TRD/trdvor_oce.F90	(revision 215)
@@ -0,0 +1,30 @@
+MODULE trdvor_oce
+   !!======================================================================
+   !!                   ***  MODULE trdvor_oce  ***
+   !! Ocean trends :   set vorticity trend variables
+   !!======================================================================
+
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE par_oce      ! ocean parameters
+
+   IMPLICIT NONE
+   PUBLIC
+
+   INTEGER,PARAMETER :: jplvor = 11     ! Number of vorticity trend terms
+
+   INTEGER, PARAMETER ::            &  !: vorticity trends index
+      jpvorprg = 1,   &  !: Pressure Gradient Trend
+      jpvorkeg = 2,   &  !: KE Gradient Trend
+      jpvorrvo = 3,   &  !: Relative Vorticity Trend
+      jpvorpvo = 4,   &  !: Planetary Vorticity Term Trend
+      jpvorldf = 5,   &  !: Horizontal Diffusion Trend
+      jpvorzad = 6,   &  !: Vertical Advection Trend
+      jpvorzdf = 7,   &  !: Vertical Diffusion Trend
+      jpvorspg = 8,   &  !: Surface Pressure Grad. Trend
+      jpvorbev = 9,   &  !: Beta V
+      jpvorswf =10,   &  !: wind stress forcing term
+      jpvorbfr =11       !: bottom friction term
+
+  !!======================================================================
+END MODULE trdvor_oce
