Index: branches/dev_001_GM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 790)
@@ -1,16 +1,13 @@
 MODULE dom_oce
+   !!======================================================================
+   !!                       ***  MODULE dom_oce  ***
+   !! ** Purpose :   Define in memory all the ocean space domain variables
+   !!======================================================================
+   !! History :  1.0  !  05-10  (A. Beckmann, G. Madec)  reactivate s-coordinate 
    !!----------------------------------------------------------------------
-   !!                       ***  MODULE dom_oce  ***
-   !!       
-   !! ** Purpose :   Define in memory all the ocean space domain variables
+   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2007) 
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-   !! History :
-   !!   9.0  !  05-10  (A. Beckmann, G. Madec)  reactivate s-coordinate 
-   !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2006) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_oce      ! ocean parameters
 
@@ -22,55 +19,49 @@
    !! space domain parameters
    !! -----------------------
-   LOGICAL, PUBLIC ::   &   !:
-      lzoom      =  .FALSE. ,   &  !: zoom flag
-      lzoom_e    =  .FALSE. ,   &  !: East  zoom type flag
-      lzoom_w    =  .FALSE. ,   &  !: West  zoom type flag
-      lzoom_s    =  .FALSE. ,   &  !: South zoom type flag
-      lzoom_n    =  .FALSE. ,   &  !: North zoom type flag
-      lzoom_arct =  .FALSE. ,   &  !: ORCA    arctic zoom flag
-      lzoom_anta =  .FALSE.        !: ORCA antarctic zoom flag
+   LOGICAL, PUBLIC ::   lzoom      =  .FALSE.      !: zoom flag
+   LOGICAL, PUBLIC ::   lzoom_e    =  .FALSE.      !: East  zoom type flag
+   LOGICAL, PUBLIC ::   lzoom_w    =  .FALSE.      !: West  zoom type flag
+   LOGICAL, PUBLIC ::   lzoom_s    =  .FALSE.      !: South zoom type flag
+   LOGICAL, PUBLIC ::   lzoom_n    =  .FALSE.      !: North zoom type flag
+   LOGICAL, PUBLIC ::   lzoom_arct =  .FALSE.      !: ORCA    arctic zoom flag
+   LOGICAL, PUBLIC ::   lzoom_anta =  .FALSE.      !: ORCA antarctic zoom flag
 
-   INTEGER, PUBLIC ::           & !!: namdom : space domain (bathymetry, mesh)
-      ntopo   =  0 ,            &  !: = 0/1 ,compute/read the bathymetry file
-      ngrid   =  0 ,            &  !: = 0/1, compute/read the horizontal mesh file
-      nmsh    =  0                 !: = 1 create a mesh-mask file
+   !                                              !!! namdom : space domain (bathymetry, mesh)
+   INTEGER, PUBLIC ::   ntopo   =  0               !: = 0/1 ,compute/read the bathymetry file
+   INTEGER, PUBLIC ::   ngrid   =  0               !: = 0/1, compute/read the horizontal mesh file
+   INTEGER, PUBLIC ::   nmsh    =  0               !: = 1 create a mesh-mask file
 
-   INTEGER, PUBLIC ::         &   !:
-      ! domain parameters linked to mpp
-      nperio,          &  !: type of lateral boundary condition
-      nimpp, njmpp,    &  !: i- & j-indexes for mpp-subdomain left bottom
-      nreci, nrecj,    &  !: overlap region in i and j
-      nproc,           &  !: number for local processor
-      narea,           &  !: number for local area
-      nbondi, nbondj,  &  !: mark of i- and j-direction local boundaries
-      npolj,           &  !: north fold mark (0, 3 or 4)
-      nlci, nlcj,      &  !: i- & j-dimensions of the local subdomain
-      nldi, nlei,      &  !: first and last indoor i- and j-indexes
-      nldj, nlej,      &  !:
-      noea, nowe,      &  !: index of the local neighboring processors in
-      noso, nono,      &  !: east, west, south and north directions
-      npne, npnw,      &  !: index of north east and north west processor
-      npse, npsw,      &  !: index of south east and south west processor
-      nbne, nbnw,      &  !: logical of north east & north west processor
-      nbse, nbsw,      &  !: logical of south east & south west processor
-      nidom               !: ???
+   !                                              !!! domain parameters linked to mpp
+   INTEGER, PUBLIC ::   nperio                     !: type of lateral boundary condition
+   INTEGER, PUBLIC ::   nimpp , njmpp              !: i- & j-indexes for mpp-subdomain left bottom
+   INTEGER, PUBLIC ::   nreci , nrecj              !: overlap region in i and j
+   INTEGER, PUBLIC ::   nproc                      !: number for local processor
+   INTEGER, PUBLIC ::   narea                      !: number for local area
+   INTEGER, PUBLIC ::   nbondi, nbondj             !: mark of i- and j-direction local boundaries
+   INTEGER, PUBLIC ::   npolj                      !: north fold mark (0, 3 or 4)
+   INTEGER, PUBLIC ::   nlci  , nlcj               !: i- & j-dimensions of the local subdomain
+   INTEGER, PUBLIC ::   nldi  , nlei               !: first and last indoor i- and j-indexes
+   INTEGER, PUBLIC ::   nldj  , nlej               !:
+   INTEGER, PUBLIC ::   noea  , nowe               !: index of the local neighboring processors in
+   INTEGER, PUBLIC ::   noso  , nono               !: east, west, south and north directions
+   INTEGER, PUBLIC ::   npne  , npnw               !: index of north east and north west processor
+   INTEGER, PUBLIC ::   npse  , npsw               !: index of south east and south west processor
+   INTEGER, PUBLIC ::   nbne  , nbnw               !: logical of north east & north west processor
+   INTEGER, PUBLIC ::   nbse  , nbsw               !: logical of south east & south west processor
+   INTEGER, PUBLIC ::   nidom                      !: ???
 
-   INTEGER, PUBLIC, DIMENSION(jpi) ::   &   !:
-      mig                 !: local  ==> global  domain i-indice
-   INTEGER, PUBLIC, DIMENSION(jpj) ::   &   !:
-      mjg                 !: local  ==> global  domain j-indice
-   INTEGER, PUBLIC, DIMENSION( jpidta ) ::   &  !:  !!bug ==> other solution?
-      mi0, mi1            !: global ==> local domain i-indice
-      !                   !  (mi0=1 and mi1=0 if the global indice is not in the local domain)
-   INTEGER, PUBLIC, DIMENSION( jpjdta ) ::   &  !:
-      mj0, mj1            !: global ==> local domain j-indice
-      !                   ! (mi0=1 and mi1=0 if the global indice is not in the local domain)
+   INTEGER, PUBLIC, DIMENSION(jpi)    ::   mig        !: local  ==> global  domain i-indice
+   INTEGER, PUBLIC, DIMENSION(jpj)    ::   mjg        !: local  ==> global  domain j-indice
+   INTEGER, PUBLIC, DIMENSION(jpidta) ::   mi0, mi1   !: global ==> local domain i-indice  (if the global indice is
+   !                                                  !  not in the local domain then mi0=1 and mi1=0)
+   !                                                  !!gm bug ==> other solution?
+   INTEGER, PUBLIC, DIMENSION(jpjdta) ::   mj0, mj1   !: global ==> local domain j-indice  (if the global indice is
+   !                                                  ! not in the local domain then mj0=1 and mj1=0)
 
-   INTEGER, PUBLIC, DIMENSION(jpnij) ::   &  !:
-      nimppt, njmppt,  &  !: i-, j-indexes for each processor
-      ibonit, ibonjt,  &  !: i-, j- processor neighbour existence
-      nlcit , nlcjt,   &  !: dimensions of every subdomain
-      nldit , nldjt,   &  !: first, last indoor index for each i-domain
-      nleit , nlejt       !: first, last indoor index for each j-domain
+   INTEGER, PUBLIC, DIMENSION(jpnij)  ::   nimppt, njmppt   !: i-, j-indexes for each processor
+   INTEGER, PUBLIC, DIMENSION(jpnij)  ::   ibonit, ibonjt   !: i-, j- processor neighbour existence
+   INTEGER, PUBLIC, DIMENSION(jpnij)  ::   nlcit , nlcjt    !: dimensions of every subdomain
+   INTEGER, PUBLIC, DIMENSION(jpnij)  ::   nldit , nldjt    !: first, last indoor index for each i-domain
+   INTEGER, PUBLIC, DIMENSION(jpnij)  ::   nleit , nlejt    !: first, last indoor index for each j-domain
 
    !!----------------------------------------------------------------------
@@ -78,24 +69,21 @@
    !! ---------------------------------------------------------------------
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
-      glamt, glamu,    &  !: longitude of t-, u-, v- and f-points (degre)
-      glamv, glamf,    &  !:
-      gphit, gphiu,    &  !: latitude  of t-, u-, v- and f-points (degre)
-      gphiv, gphif,    &  !:
-      e1t, e2t,        &  !: horizontal scale factors at t-point (m)
-      e1u, e2u,        &  !: horizontal scale factors at u-point (m)
-      e1v, e2v,        &  !: horizontal scale factors at v-point (m)
-      e1f, e2f,        &  !: horizontal scale factors at f-point (m)
-      ff                  !: coriolis factor (2.*omega*sin(yphi) ) (s-1)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   glamt, glamu   !: longitude of t-, u-, v- and f-points (degre)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   glamv, glamf   !:
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gphit, gphiu   !: latitude  of t-, u-, v- and f-points (degre)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gphiv, gphif   !:
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   e1t  , e2t     !: horizontal scale factors at t-point (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   e1u  , e2u     !: horizontal scale factors at u-point (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   e1v  , e2v     !: horizontal scale factors at v-point (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   e1f  , e2f     !: horizontal scale factors at f-point (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ff             !: coriolis factor (2.*omega*sin(yphi) ) (s-1)
 
    !!----------------------------------------------------------------------
    !! vertical coordinate and scale factors
    !! --------------------------------------
-
-   LOGICAL, PUBLIC ::           & !!: namzgr : vertical coordinate
-      ln_zco     =  .TRUE.  ,   &  !: z-coordinate - full step
-      ln_zps     =  .FALSE. ,   &  !: z-coordinate - partial step
-      ln_sco     =  .FALSE.        !: s-coordinate or hybrid z-s coordinate
-
+   !                                           !!* namzgr : vertical coordinate *
+   LOGICAL, PUBLIC ::   ln_zco     =  .TRUE.    !: z-coordinate - full step
+   LOGICAL, PUBLIC ::   ln_zps     =  .FALSE.   !: z-coordinate - partial step
+   LOGICAL, PUBLIC ::   ln_sco     =  .FALSE.   !: s-coordinate or hybrid z-s coordinate
 #if defined key_zco
    LOGICAL, PUBLIC, PARAMETER ::   lk_zco = .TRUE.    !: z-coordinate flag (1D arrays)
@@ -105,11 +93,10 @@
    !! All coordinates
    !! ---------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
-      gdep3w        ,    &  !: depth of T-points (sum of e3w) (m)
-      gdept , gdepw ,    &  !: analytical depth at T-W  points (m)
-      e3v   , e3f   ,    &  !: analytical vertical scale factors at  V--F
-      e3t   , e3u   ,    &  !:                                       T--U  points (m)
-      e3vw          ,    &  !: analytical vertical scale factors at  VW--
-      e3w   , e3uw          !:                                        W--UW  points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   gdep3w          !: depth of T-points (sum of e3w) (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   gdept , gdepw   !: depth at T-W  points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3v   , e3f     !: vertical scale factors at  V--F
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3t   , e3u     !:                            T--U  points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3vw            !: vertical scale factors at VW--
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3w   , e3uw    !:                            W--UW  points (m)
 #endif
 #if defined key_vvl
@@ -118,46 +105,41 @@
    !! All coordinates
    !! ---------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
-      gdep3w_1        ,  &  !: depth of T-points (sum of e3w) (m)
-      gdept_1, gdepw_1,  &  !: analytical depth at T-W  points (m)
-      e3v_1  , e3f_1  ,  &  !: analytical vertical scale factors at  V--F
-      e3t_1  , e3u_1  ,  &  !:                                       T--U  points (m)
-      e3vw_1          ,  &  !: analytical vertical scale factors at  VW--
-      e3w_1  , e3uw_1       !:                                       W--UW  points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   gdep3w_1          !: depth of T-points (sum of e3w) (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   gdept_1, gdepw_1  !: depth at T-W  points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3v_1  , e3f_1    !: vertical scale factors at  V--F
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3t_1  , e3u_1    !:                            T--U  points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3vw_1            !: vertical scale factors at VW--
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e3w_1  , e3uw_1   !:                            W--UW  points (m)
 #else
    LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag
 #endif
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      hur, hvr,          &  !: inverse of u and v-points ocean depth (1/m)
-      hu , hv               !: depth at u- and v-points (meters)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hur, hvr              !: inverse of u and v-points ocean depth (1/m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hu , hv               !: depth at u- and v-points (meters)
 
    !! z-coordinate with full steps (also used in the other cases as reference z-coordinate)
    !! =-----------------====------
-   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &  !:
-      gdept_0, gdepw_0,       &  !: reference depth of t- and w-points (m)
-      e3t_0  , e3w_0             !: reference vertical scale factors at T- and W-pts (m)
+   REAL(wp), PUBLIC, DIMENSION(jpk) ::   gdept_0, gdepw_0          !: reference depth of t- and w-points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpk) ::   e3t_0  , e3w_0            !: reference vertical scale factors at T- & W-pts (m)
 
    !! z-coordinate with partial steps
    !! =-----------------=======------
-   REAL(wp), PUBLIC ::      & !!: * namelist namdom *
-      e3zps_min = 5.0_wp,   &  !: miminum thickness for partial steps (meters)
-      e3zps_rat = 0.1_wp       !: minimum thickness ration for partial steps
+   !                          !!: * namelist namdom *
+   REAL(wp), PUBLIC ::   e3zps_min = 5.0_wp       !: miminum thickness for partial steps (meters)
+   REAL(wp), PUBLIC ::   e3zps_rat = 0.1_wp       !: minimum thickness ration for partial steps
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
-      hdept, hdepw, e3tp, e3wp   !: ???
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hdept, hdepw, e3tp, e3wp   !: ???
 
    !! s-coordinate and hybrid z-s-coordinate
    !! =----------------======---------------
-   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &   !:
-      gsigt, gsigw ,   &  !: model level depth coefficient at t-, w-levels (analytic)
-      gsi3w        ,   &  !: model level depth coefficient at w-level (sum of gsigw)
-      esigt, esigw        !: vertical scale factor coef. at t-, w-levels
+   REAL(wp), PUBLIC, DIMENSION(jpk) ::   gsigt, gsigw   !: model level depth coefficient at t-, w-levels (analytic)
+   REAL(wp), PUBLIC, DIMENSION(jpk) ::   gsi3w          !: model level depth coefficient at w-level (sum of gsigw)
+   REAL(wp), PUBLIC, DIMENSION(jpk) ::   esigt, esigw   !: vertical scale factor coef. at t-, w-levels
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      hbatv , hbatf ,   &  !: ocean depth at the vertical of  V--F
-      hbatt , hbatu ,   &  !:                                 T--U  points (m)
-      scosrf, scobot,   &  !:  ocean surface and bottom topographies (if deviating from coordinate surfaces in HYBRID)
-      hifv  , hiff  ,   &  !: interface depth between stretching    at  V--F
-      hift  , hifu         !: and quasi-uniform spacing                 T--U  points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hbatv , hbatf   !: ocean depth at the vertical of  V--F
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hbatt , hbatu   !:                                 T--U  points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   scosrf, scobot  !:  ocean surface and bottom topographies
+   !                                                         !   (if deviating from coordinate surfaces in HYBRID)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hifv  , hiff    !: interface depth between stretching at  V--F
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hift  , hifu    !: and quasi-uniform spacing              T--U  points (m)
 
    !!----------------------------------------------------------------------
@@ -165,25 +147,16 @@
    !! -----------------
 
-   INTEGER , PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      mbathy     !: number of ocean level (=0, 1, ... , jpk-1)
+   INTEGER , PUBLIC, DIMENSION(jpi,jpj)     ::   mbathy         !: number of ocean level (=0, 1, ... , jpk-1)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   bathy          !: ocean depth (meters)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   tmask_i        !: interior domain T-point mask
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   bmask          !: land/ocean mask of barotropic stream function
+   REAL(wp), PUBLIC, DIMENSION(jpiglo)      ::   tpol, fpol     !: north fold mask (nperio= 3 or 4)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   vmask, fmask   !: land/ocean mask at V--F
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   tmask, umask   !:             and at V--U  points
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      bathy  ,         &  !: ocean depth (meters)
-      tmask_i,         &  !: interior domain T-point mask
-      bmask               !: land/ocean mask of barotropic stream function
-
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
-      tmask, umask,    &  !: land/ocean mask at T-, U-, V- and F-points
-      vmask, fmask        !:
-
-   REAL(wp), PUBLIC, DIMENSION(jpiglo) ::   &   !:
-      tpol, fpol          !: north fold mask (nperio= 3 or 4)
 
 #if defined key_noslip_accurate
-   INTEGER, PUBLIC, DIMENSION(4,jpk) ::   &   !:
-      npcoa               !: ???
-   INTEGER, PUBLIC, DIMENSION(2*(jpi+jpj),4,jpk) ::   &   !:
-      nicoa,           &  !: ???
-      njcoa               !: ???
+   INTEGER, PUBLIC, DIMENSION(            4,jpk) ::   npcoa          !: ???
+   INTEGER, PUBLIC, DIMENSION(2*(jpi+jpj),4,jpk) ::   nicoa, njcoa   !: ???
 #endif
 
@@ -191,28 +164,25 @@
    !! time domain
    !!----------------------------------------------------------------------
-   INTEGER, PUBLIC ::    & !!: * Namelist * ???
-      nacc   = 0 ,       &  !: = 0/1 use of the acceleration of convergence technique
-      neuler                !: restart euler forward option (0=Euler)
+   !                                        !!! * Namelist * ???
+   INTEGER, PUBLIC  ::   nacc   = 0          !: = 0/1 use of the acceleration of convergence technique
+   INTEGER, PUBLIC  ::   neuler              !: restart euler forward option (0=Euler)
 
+   !                                        !!! * Namelist ??? *
+   REAL(wp), PUBLIC ::   rdt    = 3600._wp   !: time step for the dynamics (and tracer if nacc=0)
+   REAL(wp), PUBLIC ::   rdtmin = 3600._wp   !: minimum time step on tracers
+   REAL(wp), PUBLIC ::   rdtmax = 3600._wp   !: maximum time step on tracers
+   REAL(wp), PUBLIC ::   rdth   =  800._wp   !: depth variation of tracer step
+   REAL(wp), PUBLIC ::   rdtbt  =   90._wp   !: barotropic time step for the dynamics (lk_dynspg_ts=T)
+   REAL(wp), PUBLIC ::   atfp   =   0.1_wp   !: asselin time filter parameter
+   REAL(wp), PUBLIC ::   atfp1               !: asselin time filter coeff. (atfp1= 1-2*atfp)
 
-   REAL(wp), PUBLIC ::       & !!: * Namelist ??? *
-      rdt    = 3600._wp ,    &  !: time step for the dynamics (and tracer if nacc=0)
-      rdtmin = 3600._wp ,    &  !: minimum time step on tracers
-      rdtmax = 3600._wp ,    &  !: maximum time step on tracers
-      rdth   =  800._wp ,    &  !: depth variation of tracer step
-      rdtbt  =   90._wp ,    &  !: barotropic time step for the dynamics (lk_dynspg_ts=T)
-      atfp   = 0.1_wp   ,    &  !: asselin time filter parameter
-      atfp1                     !: asselin time filter coeff. (atfp1= 1-2*atfp)
-
-   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &  !:
-      rdttra                    !: vertical profile of tracer time step
+   REAL(wp), PUBLIC, DIMENSION(jpk) ::   rdttra   !: vertical profile of tracer time step
 
    !!----------------------------------------------------------------------
    !! cross land advection
    !!----------------------------------------------------------------------
-
-   INTEGER, PUBLIC ::       & !!: namelist ???
-      n_cla                    !: flag (0/1) for cross land advection to
-      !                        ! parameterize exchanges through straits
+   !                                        !!! * Namelist ??? *
+   INTEGER, PUBLIC ::   n_cla                !: flag (0/1) for cross land advection to
+      !                                      ! parameterize exchanges through straits
 
 #if defined key_agrif
@@ -221,8 +191,8 @@
    !!----------------------------------------------------------------------
    LOGICAL, PUBLIC :: spongedoneT = .FALSE.  !: ???
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &
-      zspe1ur, zspe2vr ,zspbtr2    !: ???
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   zspe1ur, zspe2vr ,zspbtr2    !: ???
    !!----------------------------------------------------------------------
 #endif
 
+   !!======================================================================
 END MODULE dom_oce
Index: branches/dev_001_GM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 790)
@@ -1,7 +1,11 @@
 MODULE domzgr
-   !!==============================================================================
+   !!======================================================================
    !!                       ***  MODULE domzgr   ***
    !! Ocean initialization : domain initialization
-   !!==============================================================================
+   !!======================================================================
+   !! History :  NEMO !  03-08  (G. Madec)  original code
+   !!            1.0  !  05-10  (A. Beckmann)  modifications for hybrid s-ccordinates
+   !!             -   !  06-04  (R. Benshila, G. Madec)  add zgr_zco routine
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -14,6 +18,5 @@
    !!       zgr_zps      : z-coordinate with partial steps
    !!       zgr_sco      : s-coordinate
-   !!---------------------------------------------------------------------
-   !! * Modules used
+   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -28,16 +31,12 @@
    PRIVATE
 
-   !! * Routine accessibility
    PUBLIC dom_zgr        ! called by dom_init.F90
 
-   !! * Module variables
-      REAL(wp) ::          & !!: Namelist nam_zgr_sco
-      sbot_min =  300.  ,  &  !: minimum depth of s-bottom surface (>0) (m)
-      sbot_max = 5250.  ,  &  !: maximum depth of s-bottom surface (= ocean depth) (>0) (m)
-      theta    =    6.0 ,  &  !: surface control parameter (0<=theta<=20)
-      thetb    =    0.75,  &  !: bottom control parameter  (0<=thetb<= 1)
-      r_max    =    0.15      !: maximum cut-off r-value allowed (0<r_max<1)
-
-
+   !                                 !!! * Namelist nam_zgr_sco *
+   REAL(wp) ::   sbot_min =  300.     !: minimum depth of s-bottom surface (>0) (m)
+   REAL(wp) ::   sbot_max = 5250.     !: maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   REAL(wp) ::   theta    =    6.0    !: surface control parameter (0<=theta<=20)
+   REAL(wp) ::   thetb    =    0.75   !: bottom control parameter  (0<=thetb<= 1)
+   REAL(wp) ::   r_max    =    0.15   !: maximum cut-off r-value allowed (0<r_max<1)
  
    !! * Substitutions
@@ -45,5 +44,7 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -58,10 +59,4 @@
       !!
       !! ** Method  reference vertical coordinate
-      !!
-      !! ** Action : 
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  original code
-      !!   9.0  !  05-10  (A. Beckmann)  modifications for hybrid s-ccordinates
       !!----------------------------------------------------------------------
       INTEGER ::   ioptio = 0      ! temporary integer
@@ -70,13 +65,8 @@
       !!----------------------------------------------------------------------
 
-      ! Read Namelist nam_zgr : vertical coordinate'
-      ! ---------------------
-      REWIND ( numnam )
+      REWIND ( numnam )                ! Read Namelist nam_zgr : vertical coordinate'
       READ   ( numnam, nam_zgr )
 
-      ! Parameter control and print
-      ! ---------------------------
-      ! Control print
-      IF(lwp) THEN
+      IF(lwp) THEN                     ! Control print
          WRITE(numout,*)
          WRITE(numout,*) 'dom_zgr : vertical coordinate'
@@ -88,5 +78,5 @@
       ENDIF
 
-      ! Check Vertical coordinate options
+      !                                ! Check Vertical coordinate options
       ioptio = 0
       IF( ln_zco ) ioptio = ioptio + 1
@@ -131,5 +121,5 @@
       ENDIF
 !!!bug gm
-
+      !
    END SUBROUTINE dom_zgr
 
@@ -154,16 +144,10 @@
       !!              -  e3t_0, e3w_0    : scale factors at T- and W-levels (m)
       !!
-      !! Reference :
-      !!      Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER  ::   jk                     ! dummy loop indices
-      REAL(wp) ::   zt, zw                 ! temporary scalars
-      REAL(wp) ::   &
-      zsur , za0, za1, zkth, zacr,      &  ! Values set from parameters in
-      zdzmin, zhmax                        ! par_CONFIG_Rxx.h90
+      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
+      !!----------------------------------------------------------------------
+      INTEGER  ::   jk                          ! dummy loop indices
+      REAL(wp) ::   zt, zw                      ! temporary scalars
+      REAL(wp) ::   zsur , za0, za1             ! value set from parameters 
+      REAL(wp) ::   zkth, zacr, zdzmin, zhmax   !    "             "
       !!----------------------------------------------------------------------
 
@@ -176,5 +160,5 @@
       !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
       !
-       IF(  ppa1  == pp_to_be_computed  .AND.  &
+      IF(   ppa1  == pp_to_be_computed  .AND.  &
          &  ppa0  == pp_to_be_computed  .AND.  &
          &  ppsur == pp_to_be_computed           ) THEN
@@ -183,16 +167,11 @@
              &                         *  (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
              &                             - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
-
          za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
          zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
-
-       ELSE
-         za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
-       ENDIF
-
-
-      ! Parameter print
-      ! ---------------
-      IF(lwp) THEN
+      ELSE
+         za1 = ppa1   ;   za0 = ppa0   ;   zsur = ppsur
+      ENDIF
+
+      IF(lwp) THEN                    ! control print
          WRITE(numout,*)
          WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
@@ -220,6 +199,5 @@
       ! Reference z-coordinate (depth - scale factor at T- and W-points)
       ! ======================
-      IF( ppkth == 0.e0 ) THEN            !  uniform vertical grid       
-
+      IF( ppkth == 0.e0 ) THEN            ! uniform vertical grid       
          za1 = zhmax / FLOAT(jpk-1) 
          DO jk = 1, jpk
@@ -231,7 +209,5 @@
             e3t_0  (jk) =  za1
          END DO
-
-      ELSE
-
+      ELSE                                ! stretched vertical grid
          DO jk = 1, jpk
             zw = FLOAT( jk )
@@ -242,12 +218,8 @@
             e3t_0  (jk) =          za0      + za1        * TANH(      (zt-zkth)/zacr   )
          END DO
-         gdepw_0(1) = 0.e0      ! force first w-level to be exactly at zero
-
-      ENDIF
-
-      ! Control and  print
-      ! ==================
-
-      IF(lwp) THEN
+         gdepw_0(1) = 0.e0                   ! force first w-level to be exactly at zero
+      ENDIF
+
+      IF(lwp) THEN                        ! Control and  print
          WRITE(numout,*)
          WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
@@ -255,10 +227,10 @@
          WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
       ENDIF
-
+      !
       DO jk = 1, jpk
          IF( e3w_0(jk)  <= 0. .OR. e3t_0(jk)  <= 0. ) CALL ctl_stop( ' e3w or e3t =< 0 ' )
          IF( gdepw_0(jk) < 0. .OR. gdept_0(jk) < 0. ) CALL ctl_stop( ' gdepw or gdept < 0 ' )
       END DO
-
+      !
    END SUBROUTINE zgr_z
 
@@ -297,24 +269,14 @@
       !! ** Action  : - mbathy: level bathymetry (in level index)
       !!              - bathy : meter bathymetry (in meters)
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  Original code
-      !!   9.0  !  05-10  (A. Beckmann)  modifications for s-ccordinates
-      !!----------------------------------------------------------------------
-      !! * Modules used
+      !!----------------------------------------------------------------------
       USE iom
 
-      !! * Local declarations
-      INTEGER ::   ji, jj, jl, jk       ! dummy loop indices
-      INTEGER ::   inum                 ! temporary logical unit
-      INTEGER  ::   &
-         ii_bump, ij_bump, ih           ! bump center position
-      INTEGER , DIMENSION(jpidta,jpjdta) ::   &
-         idta                           ! global domain integer data
-      REAL(wp) ::   &
-         r_bump, h_bump, h_oce,      &  ! bump characteristics 
-         zi, zj, zh                     ! temporary scalars
-      REAL(wp), DIMENSION(jpidta,jpjdta) ::   &
-         zdta                           ! global domain scalar data
+      INTEGER  ::   ji, jj, jl, jk          ! dummy loop indices
+      INTEGER  ::   inum                    ! temporary logical unit
+      INTEGER  ::   ii_bump, ij_bump, ih    ! bump center position
+      REAL(wp) ::   r_bump, h_bump, h_oce   ! bump characteristics 
+      REAL(wp) ::   zi, zj, zh                     ! temporary scalars
+      INTEGER , DIMENSION(jpidta,jpjdta) ::   idta   ! global domain integer data
+      REAL(wp), DIMENSION(jpidta,jpjdta) ::   zdta   ! global domain scalar data
       !!----------------------------------------------------------------------
 
@@ -331,16 +293,13 @@
 
          IF( ntopo == 0 ) THEN                        ! flat basin
-
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
-
             idta(:,:) = jpkm1                            ! flat basin 
             zdta(:,:) = gdepw_0(jpk)
             h_oce     = gdepw_0(jpk)
-
+            !
          ELSE                                         ! bump
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
-
             ii_bump = jpidta / 2           ! i-index of the bump center
             ij_bump = jpjdta / 2           ! j-index of the bump center
@@ -361,5 +320,4 @@
                END DO
             END DO
-
             ! idta :
             IF( ln_sco ) THEN      ! s-coordinate (zsc       ): idta()=jpk
@@ -478,9 +436,5 @@
       !!
       !! ** Action  : - update mbathy: level bathymetry (in level index)
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  Original code
-      !!----------------------------------------------------------------------
-      !! * Local variables
+      !!----------------------------------------------------------------------
       INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
       !!----------------------------------------------------------------------
@@ -523,5 +477,5 @@
          !
       ENDIF
-
+      !
    END SUBROUTINE zgr_bat_zoom
 
@@ -548,15 +502,8 @@
       !! ** Action  : - update mbathy: level bathymetry (in level index)
       !!              - update bathy : meter bathymetry (in meters)
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  Original code
-      !!   9.0  !  05-10  (A. Beckmann)  modifications for s-ccordinates
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::   ji, jj, jl           ! dummy loop indices
-      INTEGER ::   &
-         icompt, ibtest, ikmax          ! temporary integers
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zbathy                         ! temporary workspace
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jl                    ! dummy loop indices
+      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
+      REAL(wp), DIMENSION(jpi,jpj) ::   zbathy   ! temporary workspace
       !!----------------------------------------------------------------------
 
@@ -580,5 +527,4 @@
 
          DO jl = 1, 2
-
             IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
                mbathy( 1 ,:) = mbathy(jpim1,:)
@@ -597,11 +543,10 @@
                END DO
             END DO
-
-         END DO
-         IF( icompt == 0 ) THEN
-            IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
-         ELSE
-            IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
+         END DO
+         !
+         IF( icompt == 0 ) THEN   ;   IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
+         ELSE                     ;   IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
          ENDIF
+         !
          IF( lk_mpp ) THEN
             zbathy(:,:) = FLOAT( mbathy(:,:) )
@@ -610,5 +555,5 @@
          ENDIF
 
-         ! 3.2 East-west cyclic boundary conditions
+         ! East-west cyclic boundary conditions
 
          IF( nperio == 0 ) THEN
@@ -661,7 +606,7 @@
                mbathy(:,:) = INT( zbathy(:,:) )
             ENDIF
-
+            !
          ENDIF
-
+         !
       ENDIF
 
@@ -674,5 +619,5 @@
          END DO
       END DO
-      !!! test a faire:   ikmax = MAX( mbathy(:,:) )   ???
+!!gm      !!! test a faire:   ikmax = MAX( mbathy(:,:) )   ???
 
       IF( ikmax > jpkm1 ) THEN
@@ -694,5 +639,5 @@
          WRITE(numout,*)
       ENDIF
-
+      !
    END SUBROUTINE zgr_bat_ctl
 
@@ -705,11 +650,8 @@
       !!
       !! ** Method  :   set 3D coord. arrays to reference 1D array if lk_zco=T
-      !!
-      !! History :
-      !!        !  06-04  (R. Benshila, G. Madec)  Original code 
       !!----------------------------------------------------------------------
       INTEGER  ::   jk
       !!----------------------------------------------------------------------
-
+      !
       IF( .NOT.lk_zco ) THEN
          DO jk = 1, jpk
@@ -726,7 +668,6 @@
          END DO
       ENDIF
-
+      !
    END SUBROUTINE zgr_zco
-
 
 
@@ -769,25 +710,14 @@
       !!      schemes.
       !!
-      !! Reference :
-      !!      Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
-      !!
-      !! History :
-      !!        !  95-12  (G. Madec)  Original code : s vertical coordinate
-      !!        !  97-07  (G. Madec)  lbc_lnk call
-      !!        !  97-04  (J.-O. Beismann) 
-      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
-      !!   9.0  !  05-10  (A. Beckmann)  new stretching function
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
+      !!----------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl
-      REAL(wp) ::   fssig, fsdsig, pfkk
+      REAL(wp) ::   fssigt, fssigw, pfkk
 
       INTEGER  ::   iip1, ijp1, iim1, ijm1
-      REAL(wp) ::   &
-         fssigt, fssigw, zcoeft, zcoefw,   &
-         zrmax, ztaper
-
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zenv, ztmp, zmsk, zri, zrj, zhbat
+!     REAL(wp) ::   zcoeft, zcoefw
+      REAL(wp) ::   zrmax, ztaper
+
+      REAL(wp), DIMENSION(jpi,jpj) ::   zenv, ztmp, zmsk, zri, zrj, zhbat
 
       NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max
@@ -1120,5 +1050,6 @@
       END DO
 
-      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
+      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
+         &                                        ' MAX '               , MAXVAL( mbathy(:,:) )
 
 
@@ -1129,7 +1060,5 @@
       IF( lzoom )   CALL zgr_bat_zoom
 
-      ! 2.4 Control print
-
-      IF(lwp) THEN
+      IF(lwp) THEN                           ! Control print
          WRITE(numout,*) 
          WRITE(numout,*) ' domzgr: vertical coefficients for model level'
@@ -1202,5 +1131,5 @@
       END DO
 !!bug gm    #endif
-
+      !
    END SUBROUTINE zgr_sco
 
Index: branches/dev_001_GM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90	(revision 790)
@@ -99,8 +99,4 @@
                ua(ji,jj,jk) = ua(ji,jj,jk) + zua
                va(ji,jj,jk) = va(ji,jj,jk) + zva
-#if defined key_trddyn
-               utrd(ji,jj,jk,1) = zua                    ! save the horizontal advective trend of momentum
-               vtrd(ji,jj,jk,1) = zva 
-#endif
             END DO
          END DO
@@ -109,4 +105,5 @@
       !                                                ! ===============
 
+!!gm  momentum trend diagnostics is missing
 
       ! II. Vertical advection
@@ -158,4 +155,5 @@
          END DO
       END DO
+!!gm  momentum trend diagnostics is missing
 
    END SUBROUTINE dyn_adv_cen2
Index: branches/dev_001_GM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90	(revision 790)
@@ -191,8 +191,4 @@
                ua(ji,jj,jk) = ua(ji,jj,jk) + zua
                va(ji,jj,jk) = va(ji,jj,jk) + zva
-#if defined key_trddyn
-               utrd(ji,jj,jk,1) = zua                    ! save the horizontal advective trend of momentum
-               vtrd(ji,jj,jk,1) = zva 
-#endif
             END DO
          END DO
@@ -201,4 +197,5 @@
       !                                                ! ===============
 
+!!gm  momentum trend diagnostics is missing
 
       ! II. Vertical advection
Index: branches/dev_001_GM/NEMO/OPA_SRC/SBC/ocesbc.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/SBC/ocesbc.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/SBC/ocesbc.F90	(revision 790)
@@ -4,4 +4,9 @@
    !!                     Ocean surface boundary conditions
    !!======================================================================
+   !! History :   -   !  00-10  (O. Marti)  Original code
+   !!            1.0  !  02-12  (G. Madec)  F90: Free form and module
+   !!            2.0  !  05-11  (V. Garnier) Surface pressure gradient organization
+   !!            2.4  !  08-01  (G. Madec)  suppress dmp (internal damping contribution)
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -9,5 +14,4 @@
    !!   oce_sbc_dmp : ???
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce            ! dynamics and tracers variables
    USE dom_oce        ! ocean space domain variables
@@ -20,5 +24,4 @@
    USE flxmod         ! thermohaline fluxes
    USE flxrnf         ! runoffs forcing
-   USE tradmp         ! damping salinity trend
    USE dtatem         ! ocean temperature data
    USE dtasal         ! ocean salinity data
@@ -49,8 +52,4 @@
       qrp ,         &  !: heat flux damping (w/m2)
       erp              !: evaporation damping (kg/m2/s = mm/s)
-#if ! defined key_dynspg_rl
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
-      dmp              !: internal dampind term
-#endif
 
 #  include "domzgr_substitute.h90"
@@ -80,8 +79,4 @@
       !!        Louvain la Neuve Sea Ice Model in coupled mode
       !!
-      !! History :
-      !!   1.0  !  00-10  (O. Marti)  Original code
-      !!   2.0  !  02-12  (G. Madec)  F90: Free form and module
-      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -105,7 +100,4 @@
          emps  (:,:) = 0.e0
          erp   (:,:) = 0.e0
-#if ! defined key_dynspg_rl 
-         dmp   (:,:) = 0.e0
-#endif
       ENDIF
 
@@ -143,5 +135,5 @@
 
          ! total volume flux (use on sea-surface height)
-         emp (:,:) = fmass(:,:)  -  dmp(:,:) + runoff(:,:) + erp(:,:)
+         emp (:,:) = fmass(:,:)  + runoff(:,:) + erp(:,:)
 #else
          ! Rigid-lid (emp=emps=E-P-R+Erp) 
@@ -216,7 +208,4 @@
          emps   (:,:) = 0.e0
          erp    (:,:) = 0.e0
-#if ! defined key_dynspg_rl 
-         dmp    (:,:) = 0.e0
-#endif
       ENDIF
 #if defined key_flx_core
@@ -243,5 +232,5 @@
 
          ! total volume flux (use on sea-surface height)
-         emp (:,:) = fmass(:,:) -   dmp(:,:) + runoff(:,:) + erp(:,:) + empold      
+         emp (:,:) = fmass(:,:)              + runoff(:,:) + erp(:,:) + empold      
 #else
          ! Rigid-lid (emp=emps=E-P-R+Erp)
@@ -446,7 +435,4 @@
          emps   (:,:) = 0.e0
          erp    (:,:) = 0.e0
-#if ! defined key_dynspg_rl 
-         dmp    (:,:) = 0.e0
-#endif
       ENDIF
 
@@ -537,9 +523,8 @@
 
          ! Contribution to sea level:
-         ! net upward water flux emp() = e-p + runoff() + erp() + dmp + empold
+         ! net upward water flux emp() = e-p + runoff() + erp() + empold
          emp (:,:) = zemp(:,:)     &   ! e-p data
             &      + runoff(:,:)   &   ! runoff data
             &      + erp(:,:)      &   ! restoring term to SSS data
-            &      + dmp(:,:)      &   ! freshwater flux associated with internal damping
             &      + empold            ! domain averaged annual mean correction
 
@@ -809,8 +794,4 @@
       REAL(wp) ::   zerpplus(jpi,jpj), zerpminus(jpi,jpj)
       REAL(wp) ::   zplus, zminus, zadefi
-# if defined key_tradmp
-      INTEGER jk
-      REAL(wp), DIMENSION(jpi,jpj) ::   zstrdmp
-# endif
 #endif
       !!----------------------------------------------------------------------
@@ -838,17 +819,4 @@
       ! Free-surface
          
-      ! Internal damping
-# if defined key_tradmp
-      ! Vertical mean of dampind trend (computed in tradmp module)
-      zstrdmp(:,:) = 0.e0
-      DO jk = 1, jpk
-         zstrdmp(:,:) = zstrdmp(:,:) + strdmp(:,:,jk) * fse3t(:,:,jk)
-      END DO
-      ! volume flux associated to internal damping to climatology
-      dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + 1.e-20 )
-# else
-      dmp(:,:) = 0.e0            ! No internal damping
-# endif
-      
       !   evaporation damping term ( Surface restoring )
       zerpplus (:,:) = 0.e0
Index: branches/dev_001_GM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 790)
@@ -152,7 +152,7 @@
       ! -----------------------------------------------------------------
       ! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
-#  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = 1, jpij   ! vector opt. (forced unrolling)
+#  if defined key_vectopt_loop
+      DO jj = 1, 1         ! vector opt.
+         DO ji = 1, jpij   ! forced the loop collapse
 #  else
       DO jj = 1, jpj
@@ -165,13 +165,11 @@
             zsbb(ji,jj) = sb(ji,jj,ik) * tmask(ji,jj,1)
             zdep(ji,jj) = fsdept(ji,jj,ik)                ! depth of the ocean bottom T-level
-#  if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#  endif
+         END DO
       END DO
 
       IF( ln_zps ) THEN      ! partial steps correction 
-# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+# if defined key_vectopt_loop
+      DO jj = 1, 1             ! vector opt.
+         DO ji = 1, jpij-jpi   ! forced the loop collapse
 # else
          DO jj = 1, jpjm1
@@ -186,12 +184,10 @@
                zahu(ji,jj) = atrbbl * e2u(ji,jj) * ze3u / e1u(ji,jj) * umask(ji,jj,1)
                zahv(ji,jj) = atrbbl * e1v(ji,jj) * ze3v / e2v(ji,jj) * vmask(ji,jj,1)
-# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
             END DO
-# endif
          END DO
       ELSE                    ! z-coordinate - full steps or s-coordinate
 #   if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+      DO jj = 1, 1             ! vector opt.
+         DO ji = 1, jpij-jpi   ! forced the loop collapse
 #   else
          DO jj = 1, jpjm1
@@ -202,7 +198,5 @@
                zahu(ji,jj) = atrbbl * e2u(ji,jj) * fse3u(ji,jj,iku) / e1u(ji,jj) * umask(ji,jj,1)
                zahv(ji,jj) = atrbbl * e1v(ji,jj) * fse3v(ji,jj,ikv) / e2v(ji,jj) * vmask(ji,jj,1)
-#   if ! defined key_vectopt_loop   ||   defined key_mpp_omp
             END DO
-#   endif
          END DO
       ENDIF
@@ -218,6 +212,6 @@
 
 #  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+      DO jj = 1, 1             ! vector opt.
+         DO ji = 1, jpij-jpi   ! forced the loop collapse
 #  else
       DO jj = 1, jpjm1
@@ -236,12 +230,10 @@
             zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
             zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
-#  if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#  endif
-      END DO
-
-#  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+         END DO
+      END DO
+
+#  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
+      DO jj = 1, 1             ! vector opt.
+         DO ji = 1, jpij-jpi   ! forced the loop collapse
 #  else
       DO jj = 1, jpjm1
@@ -260,7 +252,5 @@
             zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
             zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
-#  if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#  endif
+         END DO
       END DO
 
@@ -268,6 +258,6 @@
                                ! 
 #  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+      DO jj = 1, 1             ! vector opt.
+         DO ji = 1, jpij-jpi   ! forced the loop collapse
 #  else
       DO jj = 1, jpjm1
@@ -279,12 +269,10 @@
             zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
             zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
-#  if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#  endif
-      END DO
-
-#  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+         END DO
+      END DO
+
+#  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
+      DO jj = 1, 1             ! vector opt.
+         DO ji = 1, jpij-jpi   ! forced the loop collapse
 #  else
       DO jj = 1, jpjm1
@@ -296,7 +284,5 @@
             zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
             zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
-#  if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#  endif
+         END DO
       END DO
 
@@ -304,6 +290,6 @@
 
 #  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+      DO jj = 1, 1             ! vector opt.
+         DO ji = 1, jpij-jpi   ! forced the loop collapse
 #  else
       DO jj = 1, jpjm1
@@ -312,16 +298,14 @@
             ! local density gradient along i-bathymetric slope
             zgdrho = - ( rbeta*( zsnb(ji+1,jj) - zsnb(ji,jj) )   &
-                     -  ralpha*( ztnb(ji+1,jj) - ztnb(ji,jj) ) )
+               &     -  ralpha*( ztnb(ji+1,jj) - ztnb(ji,jj) ) )
             ! sign of local i-gradient of density multiplied by the i-slope
             zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
             zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
-#  if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#  endif
-      END DO
-
-#  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+         END DO
+      END DO
+
+#  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
+      DO jj = 1, 1             ! vector opt.
+         DO ji = 1, jpij-jpi   ! forced the loop collapse
 #  else
       DO jj = 1, jpjm1
@@ -330,11 +314,9 @@
             ! local density gradient along j-bathymetric slope
             zgdrho = - ( rbeta*( zsnb(ji,jj+1) - zsnb(ji,jj) )   &
-                     -  ralpha*( ztnb(ji,jj+1) - ztnb(ji,jj) ) )   
+               &     -  ralpha*( ztnb(ji,jj+1) - ztnb(ji,jj) ) )   
             ! sign of local j-gradient of density multiplied by the j-slope
-            zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
+            zsign = SIGN( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
             zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
-#  if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#  endif
+         END DO
       END DO
       
@@ -351,6 +333,6 @@
       ! first derivative (gradient)
 #  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+      DO jj = 1, 1             ! vector opt.
+         DO ji = 1, jpij-jpi   ! forced the loop collapse
 #  else
       DO jj = 1, jpjm1
@@ -362,7 +344,5 @@
             zky(ji,jj) = zkj(ji,jj) * ( ztbb(ji,jj+1) - ztbb(ji,jj) )
             zkw(ji,jj) = zkj(ji,jj) * ( zsbb(ji,jj+1) - zsbb(ji,jj) )
-#  if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#  endif
+         END DO
       END DO
 
@@ -401,6 +381,6 @@
       ! second derivative (divergence) and add to the general tracer trend
 #  if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+      DO jj = 1, 1                   ! vector opt.
+         DO ji = jpi+2, jpij-jpi-1   ! forced the loop collapse
 #  else
       DO jj = 2, jpjm1
@@ -415,7 +395,5 @@
             ta(ji,jj,ik) = ta(ji,jj,ik) + zta
             sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
-#  if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#  endif
+         END DO
       END DO
 
@@ -441,6 +419,6 @@
    !!   Default option :                 NO advective bottom boundary layer
    !!----------------------------------------------------------------------
-   SUBROUTINE tra_bbl_adv (kt )              ! Empty routine
-      INTEGER, INTENT(in) :: kt
+   SUBROUTINE tra_bbl_adv( kt )              ! Empty routine
+      INTEGER :: kt
       WRITE(*,*) 'tra_bbl_adv: You should not have seen this print! error?', kt
    END SUBROUTINE tra_bbl_adv
Index: branches/dev_001_GM/NEMO/OPA_SRC/TRA/trabbl_adv.h90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/TRA/trabbl_adv.h90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/TRA/trabbl_adv.h90	(revision 790)
@@ -99,7 +99,7 @@
       ! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
 
-#if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = 1, jpij   ! vector opt. (forced unrolling)
+#if defined key_vectopt_loop 
+      DO jj = 1, 1
+         DO ji = 1, jpij   ! vector opt. (forced unrolling)
 #else
       DO jj = 1, jpj
@@ -115,7 +115,5 @@
             zunb(ji,jj) = un(ji,jj,mbku(ji,jj)) 
             zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
-#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#endif
+         END DO
       END DO
 
@@ -224,12 +222,11 @@
 
       ! ... is equal to zero but where bbl will work.
-      u_bbl(:,:,:) = 0.e0
-      v_bbl(:,:,:) = 0.e0
+      u_bbl(:,:,:) = 0.e0   ;   v_bbl(:,:,:) = 0.e0
 
       IF( ln_zps ) THEN     ! partial steps correction   
       
-# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+# if defined key_vectopt_loop
+         DO jj = 1, 1
+            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
 # else   
          DO jj = 1, jpjm1
@@ -249,7 +246,5 @@
                   v_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv) * ze3v / fse3v(ji,jj,ikv)       
                ENDIF
-# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-            END DO
-# endif
+            END DO
          END DO
 
@@ -259,7 +254,7 @@
       ELSE       ! if not partial step loop over the whole domain no lbc call
 
-#if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-         jj = 1
-         DO ji = 1, jpij   ! vector opt. (forced unrolling)
+#if defined key_vectopt_loop
+         DO jj = 1, 1
+            DO ji = 1, jpij   ! vector opt. (forced unrolling)
 #else
          DO jj = 1, jpj
@@ -272,7 +267,5 @@
                   v_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv)       
                ENDIF
-#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-            END DO
-#endif
+            END DO
          END DO
        
@@ -284,7 +277,7 @@
       ! ... Second order centered tracer flux at u and v-points
 
-# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+# if defined key_vectopt_loop
+      DO jj = 1, 1
+         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
 # else
       DO jj = 1, jpjm1
@@ -309,11 +302,9 @@
             zwz(ji,jj) = ( ( zfvj + ABS( zfvj ) ) * zsbb(ji  ,jj  )   &
                &          +( zfvj - ABS( zfvj ) ) * zsbb(ji  ,jj+1) ) * 0.5
-#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#endif
-        END DO
-# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+         END DO
+      END DO
+# if defined key_vectopt_loop
+      DO jj = 1, 1
+         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
 # else
       DO jj = 2, jpjm1
@@ -331,7 +322,5 @@
             ta(ji,jj,ik) = ta(ji,jj,ik) + zta
             sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
-#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-#endif
+         END DO
       END DO
 
@@ -365,7 +354,7 @@
       IF( ln_zps ) THEN
 	  
-# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+# if defined key_vectopt_loop 
+         DO jj = 1, 1
+            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
 # else
          DO jj = 1, jpjm1
@@ -383,14 +372,12 @@
                zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * ze3u  
                zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * ze3v
-#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-            END DO
-#endif
+            END DO
          END DO  
 	
       ELSE
 
-# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+# if defined key_vectopt_loop
+         DO jj = 1, 1
+            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
 # else
          DO jj = 1, jpjm1
@@ -401,7 +388,5 @@
                zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,iku) 
                zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,ikv) 
-#if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-            END DO
-#endif
+            END DO
          END DO
 
@@ -409,7 +394,7 @@
  
 
-# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
-      jj = 1
-      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+# if defined key_vectopt_loop
+      DO jj = 1, 1
+         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
 # else
       DO jj = 2, jpjm1
@@ -425,8 +410,6 @@
                &   ) / zbt
 
-# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
-         END DO
-# endif
-        END DO
+         END DO
+      END DO
 
       ! 7. compute additional vertical velocity to be used in t boxes
Index: branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdicp.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 790)
@@ -36,5 +36,5 @@
 
    INTERFACE trd_icp
-      MODULE PROCEDURE trd_2d, trd_3d
+      MODULE PROCEDURE trd_2d, trd_3d, trd_u2d, trd_u3d
    END INTERFACE
 
@@ -44,4 +44,19 @@
    PUBLIC   trd_icp_init  ! called by opa.F90
 
+   !! Variables used for diagnostics
+   REAL(wp) ::   tvolt        !: volume of the whole ocean computed at t-points
+   REAL(wp) ::   tvolu        !: volume of the whole ocean computed at u-points
+   REAL(wp) ::   tvolv        !: volume of the whole ocean computed at v-points
+
+   !! Active Tracer trend diagnostics variables
+   REAL(wp), DIMENSION(jpt_trd,2) ::   tsmo         !: tracers trends average 
+   REAL(wp), DIMENSION(jpt_trd,2) ::   ts2          !: tracers square trends average 
+
+   !! Momentum trends diagnostics variables
+   REAL(wp), DIMENSION(jptot_dyn) ::   umo, vmo         !: momentum trends average 
+   REAL(wp), DIMENSION(jptot_dyn) ::   hke              !: momentum square trends average 
+   REAL(wp) ::   rpktrd   !: potential to kinetic energy conversion
+   REAL(wp) ::   peke     !: conversion potential energy - kinetic energy trend
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -55,5 +70,5 @@
 CONTAINS
 
-   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype, clpas )
+   SUBROUTINE trd_2d( ptrd, ktra, ktrd , ctype, clpas )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_2d  ***
@@ -62,94 +77,167 @@
       !!              momentum equations at every time step frequency ntrd.
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dx             ! Temperature or U trend 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dy             ! Salinity    or V trend
-      INTEGER                     , INTENT(in   ) ::   ktrd                ! tracer trend index
-      CHARACTER(len=3)            , INTENT(in   ) ::   ctype               ! momentum ('DYN') or tracers ('TRA') trends
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptrd      ! tracer or u trend
+      INTEGER                     , INTENT(in   )           ::   ktra      ! tracer index
+      INTEGER                     , INTENT(in   )           ::   ktrd      ! tracer trend index
+      CHARACTER(len=3)            , INTENT(in   )           ::   ctype     ! tracer type (='TRA' or 'TRC')
       CHARACTER(len=3)            , INTENT(in   ), OPTIONAL ::   clpas     ! number of passage
       !!
-      INTEGER  ::   ji, jj                                                 ! loop indices
-      CHARACTER(len=3) ::   cpas                                           ! number of passage
-      REAL(wp) ::   zmsku, zbtu, zbt                                       ! temporary scalars
-      REAL(wp) ::   zmskv, zbtv                                            !    "         "
-      !!----------------------------------------------------------------------
-
-      ! Control of optional arguments
-      cpas = 'fst'
+      CHARACTER(len=3) ::   cpas   ! number of passage
+      REAL(wp)         ::   zsum   ! temporary scalars
+      !!----------------------------------------------------------------------
+
+      cpas = 'fst'                             ! Control of optional arguments
       IF( PRESENT(clpas) )  cpas = clpas
+
+      ! 1. Mask volumic trends
+      ptrd(:,:) = e1t(:,:) * e2t(:,:) * fse3t(:,:,1) * ptrd(:,:) * tmask_i(:,:)
+      
+      ! 2. Basin averaged tracer trends
+      SELECT CASE( ctype )
+      CASE( 'TRA' )              ! Tracers
+         zsum = SUM( ptrd(:,:) )
+         IF( cpas == 'fst' ) THEN   ;   tsmo(ktrd,ktra) =                   zsum
+         ELSE                       ;   tsmo(ktrd,ktra) = tsmo(ktrd,ktra) + zsum
+         ENDIF
+      CASE( 'TRC' )              ! Passive tracers
+         !      ....  to be done
+      END SELECT
+      
+      ! 3. Basin averaged tracer square trends (i.e. trd * now tracer field)
+      SELECT CASE( ctype )
+      CASE( 'TRA' )              ! Tracers
+         IF( ktra == jp_tem )   zsum = SUM(  ptrd(:,:) * tn(:,:,1) )
+         IF( ktra == jp_sal )   zsum = SUM(  ptrd(:,:) * sn(:,:,1) )
+         IF( cpas == 'fst' ) THEN   ;   ts2(ktrd,ktra) =                  zsum
+         ELSE                       ;   ts2(ktrd,ktra) = ts2(ktrd,ktra) + zsum
+         ENDIF
+      CASE( 'TRC' )              ! Passive tracers
+         !      ....  to be done
+      END SELECT
+      !
+   END SUBROUTINE trd_2d
+
+
+   SUBROUTINE trd_3d( ptrd, ktra, ktrd , ctype, clpas )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_3d  ***
+      !! 
+      !! ** Purpose : verify the basin averaged properties of tracers and/or
+      !!              momentum equations at every time step frequency ntrd.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::   ptrd      ! tracer or u trend
+      INTEGER                         , INTENT(in   )           ::   ktra      ! tracer index
+      INTEGER                         , INTENT(in   )           ::   ktrd      ! tracer trend index
+      CHARACTER(len=3)                , INTENT(in   )           ::   ctype     ! tracer type (='TRA' or 'TRC')
+      CHARACTER(len=3)                , INTENT(in   ), OPTIONAL ::   clpas     ! number of passage
+      !!
+      INTEGER          ::   jk     ! dummy loop indices
+      CHARACTER(len=3) ::   cpas   ! number of passage
+      REAL(wp)         ::   zsum   ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsurf          ! 2D workspace array
+      !!----------------------------------------------------------------------
+
+      cpas = 'fst'                             ! Control of optional arguments
+      IF( PRESENT(clpas) )  cpas = clpas
+
+      ! 1. Mask volumic trends
+      zsurf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
+      DO jk = 1, jpk
+         ptrd(:,:,jk) = zsurf(:,:) * fse3t(:,:,jk) * ptrd(:,:,jk) * tmask(:,:,jk)
+      END DO
+
+      ! 2. Basin averaged tracer trends
+      SELECT CASE( ctype )
+      CASE( 'TRA' )              ! Tracers
+         zsum = SUM( ptrd(:,:,:) )
+         IF( cpas == 'fst' ) THEN   ;   tsmo(ktrd,ktra) =                   zsum
+         ELSE                       ;   tsmo(ktrd,ktra) = tsmo(ktrd,ktra) + zsum
+         ENDIF
+      CASE( 'TRC' )              ! Passive tracers
+         !      ....  to be done
+      END SELECT
+
+      ! 3. Basin averaged tracer square trends (i.e. trd * now tracer field)
+      SELECT CASE( ctype )
+      CASE( 'TRA' )              ! Tracers
+         IF( ktra == jp_tem )   zsum = SUM(  ptrd(:,:,:) * tn(:,:,:) )
+         IF( ktra == jp_sal )   zsum = SUM(  ptrd(:,:,:) * sn(:,:,:) )
+         IF( cpas == 'fst' ) THEN   ;   ts2(ktrd,ktra) =                  zsum
+         ELSE                       ;   ts2(ktrd,ktra) = ts2(ktrd,ktra) + zsum
+         ENDIF
+      CASE( 'TRC' )              ! Passive tracers
+         !      ....  to be done
+      END SELECT
+      !
+   END SUBROUTINE trd_3d
+
+
+   SUBROUTINE trd_u2d( ptrdu, ptrdv, ktrd , ctype, clpas )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_2d  ***
+      !! 
+      !! ** Purpose : verify the basin averaged properties of tracers and/or
+      !!              momentum equations at every time step frequency ntrd.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptrdu, ptrdv   ! U and V momentum trends 
+      INTEGER                     , INTENT(in   )           ::   ktrd           ! tracer trend index
+      CHARACTER(len=3)            , INTENT(in   )           ::   ctype          ! momentum ('DYN')
+      CHARACTER(len=3)            , INTENT(in   ), OPTIONAL ::   clpas          ! number of passage
+      !!
+      INTEGER          ::   ji, jj        ! loop indices
+      CHARACTER(len=3) ::   cpas          ! number of passage
+      REAL(wp)         ::   zmsku, zbtu   ! temporary scalars
+      REAL(wp)         ::   zmskv, zbtv   !    "         "
+      !!----------------------------------------------------------------------
+
+      cpas = 'fst'                           ! Control of optional arguments
+      IF( PRESENT(clpas) )  cpas = clpas
+
+      SELECT CASE( ctype )
+      !
+      CASE( 'DYN' )              ! Momentum
 
       ! 1. Mask trends
       ! --------------
-
-      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
+      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)
+            ptrdu(ji,jj) = ptrdu(ji,jj) * zmsku
+            ptrdv(ji,jj) = ptrdv(ji,jj) * zmskv
+         END DO
+      END DO
+      ptrdu(jpi, : ) = 0.e0      ;      ptrdv(jpi, : ) = 0.e0
+      ptrdu( : ,jpj) = 0.e0      ;      ptrdv( : ,jpj) = 0.e0
       
       ! 2. Basin averaged tracer/momentum trends
       ! ----------------------------------------
-
-      SELECT CASE( ctype )
-      !
-      CASE( 'DYN' )              ! Momentum
-         umo(ktrd) = 0.e0
-         vmo(ktrd) = 0.e0
-         !
-         SELECT CASE( ktrd )
-         !
-         CASE( jpdyn_trd_swf )         ! 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
+      umo(ktrd) = 0.e0
+      vmo(ktrd) = 0.e0
+      !
+      SELECT CASE( ktrd )
+      !
+      CASE( jpdyn_trd_swf )         ! surface forcing
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               umo(ktrd) = umo(ktrd) + ptrdu(ji,jj) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
+               vmo(ktrd) = vmo(ktrd) + ptrdv(ji,jj) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
+            END DO
+         END DO
             !
          CASE( jpdyn_trd_bfr )         ! 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)
+                  umo(ktrd) = umo(ktrd) + ptrdu(ji,jj)
+                  vmo(ktrd) = vmo(ktrd) + ptrdv(ji,jj)
                END DO
             END DO
             !
-         END SELECT
-         !
-      CASE( 'TRA' )              ! Tracers
-         IF( cpas == 'fst' )   THEN
-            tmo(ktrd) = 0.e0
-            smo(ktrd) = 0.e0
-         ENDIF
-         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/momentum square trends
       ! ----------------------------------------------
       ! c a u t i o n: field now
       
-      SELECT CASE( ctype )
-      !
-      CASE( 'DYN' )              ! Momentum
          hke(ktrd) = 0.e0
          DO jj = 1, jpj
@@ -157,29 +245,15 @@
                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
+               hke(ktrd) = hke(ktrd) + un(ji,jj,1) * ptrdu(ji,jj) * zbtu   &
+               &                     + vn(ji,jj,1) * ptrdv(ji,jj) * zbtv
             END DO
          END DO
          !
-      CASE( 'TRA' )              ! Tracers
-         IF( cpas == 'fst' )   THEN
-            t2(ktrd) = 0.e0
-            s2(ktrd) = 0.e0
-         ENDIF
-         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, clpas )
+   END SUBROUTINE trd_u2d
+
+
+   SUBROUTINE trd_u3d( ptrdu, ptrdv, ktrd, ctype, clpas )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE trd_3d  ***
@@ -188,25 +262,23 @@
       !!              momentum equations at every time step frequency ntrd.
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dx            ! Temperature or U trend 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dy            ! Salinity    or V trend
-      INTEGER,                          INTENT(in   ) ::   ktrd               ! momentum or tracer trend index
-      CHARACTER(len=3),                 INTENT(in   ) ::   ctype              ! momentum ('DYN') or tracers ('TRA') trends
-      CHARACTER(len=3),                 INTENT(in   ), OPTIONAL ::   clpas    ! number of passage
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::   ptrdu, ptrdv   ! U and V momentum trends
+      INTEGER,                          INTENT(in   )           ::   ktrd           ! momentum trend index
+      CHARACTER(len=3),                 INTENT(in   )           ::   ctype          ! momentum (='DYN')
+      CHARACTER(len=3),                 INTENT(in   ), OPTIONAL ::   clpas          ! number of passage
       !!
       INTEGER ::   ji, jj, jk
-      CHARACTER(len=3) ::   cpas                                              ! number of passage
-      REAL(wp) ::   zbt, zbtu, zbtv, zmsku, zmskv                             ! temporary scalars
-      !!----------------------------------------------------------------------
-
-      ! Control of optional arguments
-      cpas = 'fst'
+      CHARACTER(len=3) ::   cpas               ! number of passage
+      REAL(wp) ::   zbtu, zbtv, zmsku, zmskv   ! temporary scalars
+      !!----------------------------------------------------------------------
+
+      cpas = 'fst'                            ! Control of optional arguments
       IF( PRESENT(clpas) )  cpas = clpas
 
-      ! 1. Mask the trends
-      ! ------------------
-
       SELECT CASE( ctype )
       !
       CASE( 'DYN' )              ! Momentum        
+
+         ! 1. Mask the trends
+         ! ------------------
          DO jk = 1, jpk
             DO jj = 1, jpjm1
@@ -214,26 +286,16 @@
                   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
+                  ptrdu(ji,jj,jk) = ptrdu(ji,jj,jk) * zmsku
+                  ptrdv(ji,jj,jk) = ptrdv(ji,jj,jk) * zmskv
                END DO
             END DO
          END DO
-         ptrd3dx(jpi, : ,:) = 0.e0      ;      ptrd3dy(jpi, : ,:) = 0.e0
-         ptrd3dx( : ,jpj,:) = 0.e0      ;      ptrd3dy( : ,jpj,:) = 0.e0
+         ptrdu(jpi, : ,:) = 0.e0      ;      ptrdv(jpi, : ,:) = 0.e0
+         ptrdu( : ,jpj,:) = 0.e0      ;      ptrdv( : ,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(:,:)
-         END DO
-         !
-      END SELECT   
 
       ! 2. Basin averaged tracer/momentum trends
       ! ----------------------------------------
       
-      SELECT CASE( ctype )
-      !
-      CASE( 'DYN' )              ! Momentum
          umo(ktrd) = 0.e0
          vmo(ktrd) = 0.e0
@@ -243,26 +305,10 @@
                   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
+                  umo(ktrd) = umo(ktrd) + ptrdu(ji,jj,jk) * zbtu
+                  vmo(ktrd) = vmo(ktrd) + ptrdv(ji,jj,jk) * zbtv
                END DO
             END DO
          END DO
          !
-      CASE( 'TRA' )              ! Tracers
-         IF( cpas == 'fst' )   THEN
-            tmo(ktrd) = 0.e0
-            smo(ktrd) = 0.e0
-         ENDIF
-         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 trends
@@ -270,7 +316,4 @@
       ! c a u t i o n: field now
       
-      SELECT CASE( ctype )
-      !
-      CASE( 'DYN' )              ! Momentum
          hke(ktrd) = 0.e0
          DO jk = 1, jpk
@@ -279,22 +322,6 @@
                   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
-         IF( cpas == 'fst' )   THEN
-            t2(ktrd) = 0.e0
-            s2(ktrd) = 0.e0
-         ENDIF
-         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)
+                  hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * ptrdu(ji,jj,jk) * zbtu   &
+                  &                     + vn(ji,jj,jk) * ptrdv(ji,jj,jk) * zbtv
                END DO
             END DO
@@ -303,6 +330,5 @@
       END SELECT
       !
-   END SUBROUTINE trd_3d
-
+   END SUBROUTINE trd_u3d
 
 
@@ -313,9 +339,10 @@
       !! ** Purpose :   Read the namtrd namelist
       !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zmskt
+      INTEGER  ::   jk
 #if  defined key_trddyn
+      INTEGER  ::   ji, jj
       REAL(wp) ::   zmsku, zmskv
 #endif
+      REAL(wp), DIMENSION(jpi,jpj) ::   zsurf          ! 2D workspace array
       !!----------------------------------------------------------------------
 
@@ -328,11 +355,7 @@
       ! Total volume at t-points:
       tvolt = 0.e0
+      zsurf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
       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
+            tvolt = tvolt + SUM( tmask(:,:,jk) * zsurf(:,:) * fse3t(:,:,jk) )
       END DO
       IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
@@ -379,5 +402,6 @@
       !!
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   ze1e2w, zcof, zbe1ru, zbe2rv, zbtr, ztz, zth   !    "      scalars
+      REAL(wp) ::   zcof, zbe1ru, zbe2rv, zbtr, ztz, zth   !    "      scalars
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zsurf          ! 2D workspace array
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zkepe, zkx, zky, zkz   ! temporary arrays
       !!----------------------------------------------------------------------
@@ -400,13 +424,10 @@
 
          ! Density flux at w-point
+         zsurf(:,:) = 0.5 * e1t(:,:) * e2t(:,:) * tmask_i(:,:) /rau0
+!!gm better use bn2....
+         zkz(:,:,1) = 0.e0
          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
+            zkz(:,:,jk) = zsurf(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) )
+         END DO
          
          ! Density flux at u and v-points
@@ -588,8 +609,6 @@
          ! -------------------------------
          IF( lk_mpp ) THEN
-            CALL mpp_sum( tmo, jptot_tra )   
-            CALL mpp_sum( smo, jptot_tra )
-            CALL mpp_sum( t2 , jptot_tra )
-            CALL mpp_sum( s2 , jptot_tra )
+            CALL mpp_sum( tsmo, jpt_trd )   
+            CALL mpp_sum( ts2 , jpt_trd )
          ENDIF
 
@@ -601,21 +620,26 @@
             WRITE (numout,*)
             WRITE (numout,9400) kt
-            WRITE (numout,9401) (tmo(jpicpt_xad)+tmo(jpicpt_yad))/ tvolt, (smo(jpicpt_xad)+smo(jpicpt_yad))/ tvolt
-            WRITE (numout,9402)  tmo(jpicpt_zad) / tvolt, smo(jpicpt_zad) / tvolt
-            WRITE (numout,9403)  tmo(jpicpt_ldf) / tvolt, smo(jpicpt_ldf) / tvolt
-            WRITE (numout,9404)  tmo(jpicpt_zdf) / tvolt, smo(jpicpt_zdf) / tvolt
-            WRITE (numout,9405)  tmo(jpicpt_npc) / tvolt, smo(jpicpt_npc) / tvolt
-            WRITE (numout,9406)  tmo(jpicpt_dmp) / tvolt, smo(jpicpt_dmp) / tvolt
-            WRITE (numout,9407)  tmo(jpicpt_qsr) / tvolt
-            WRITE (numout,9408)  tmo(jpicpt_nsr) / tvolt, smo(jpicpt_nsr) / tvolt
+            WRITE (numout,9401) ( tsmo(jpt_trd_xad,jp_tem) + tsmo(jpt_trd_yad,jp_tem)  ) / tvolt,   &
+               &                ( tsmo(jpt_trd_xad,jp_sal) + tsmo(jpt_trd_yad,jp_sal)  ) / tvolt
+            WRITE (numout,9402)   tsmo(jpt_trd_zad,jp_tem) / tvolt,   tsmo(jpt_trd_zad,jp_sal) / tvolt
+            WRITE (numout,9403)   tsmo(jpt_trd_ldf,jp_tem) / tvolt,   tsmo(jpt_trd_ldf,jp_sal) / tvolt
+            WRITE (numout,9404)   tsmo(jpt_trd_zdf,jp_tem) / tvolt,   tsmo(jpt_trd_zdf,jp_sal) / tvolt
+            WRITE (numout,9405)   tsmo(jpt_trd_npc,jp_tem) / tvolt,   tsmo(jpt_trd_npc,jp_sal) / tvolt
+            WRITE (numout,9406)   tsmo(jpt_trd_dmp,jp_tem) / tvolt,   tsmo(jpt_trd_dmp,jp_sal) / tvolt
+            WRITE (numout,9407)   tsmo(jpt_trd_qsr,jp_tem) / tvolt
+            WRITE (numout,9408)   tsmo(jpt_trd_qns,jp_tem) / tvolt,   tsmo(jpt_trd_qns,jp_sal) / tvolt
             WRITE (numout,9409) 
-            WRITE (numout,9410) (  tmo(jpicpt_xad) + tmo(jpicpt_yad) + tmo(jpicpt_zad) + tmo(jpicpt_ldf) + tmo(jpicpt_zdf)   &
-            &                    + tmo(jpicpt_npc) + tmo(jpicpt_dmp) + tmo(jpicpt_qsr) + tmo(jpicpt_nsr) ) / tvolt,   &
-            &                   (  smo(jpicpt_xad) + smo(jpicpt_yad) + smo(jpicpt_zad) + smo(jpicpt_ldf) + smo(jpicpt_zdf)   &
-            &                    + smo(jpicpt_npc) + smo(jpicpt_dmp)                   + smo(jpicpt_nsr) ) / tvolt
-         ENDIF
-
-9400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
-              '              salinity',/' ============================')
+            WRITE (numout,9410) (  tsmo(jpt_trd_xad,jp_tem) + tsmo(jpt_trd_yad,jp_tem) + tsmo(jpt_trd_zad,jp_tem)     &
+               &                 + tsmo(jpt_trd_ldf,jp_tem) + tsmo(jpt_trd_zdf,jp_tem) + tsmo(jpt_trd_npc,jp_tem)     &
+               &                 + tsmo(jpt_trd_dmp,jp_tem) + tsmo(jpt_trd_qsr,jp_tem) + tsmo(jpt_trd_qns,jp_tem) )   &
+               &                / tvolt,   &
+               &                (  tsmo(jpt_trd_xad,jp_sal) + tsmo(jpt_trd_yad,jp_sal) + tsmo(jpt_trd_zad,jp_sal)     &
+               &                 + tsmo(jpt_trd_ldf,jp_sal) + tsmo(jpt_trd_zdf,jp_sal) + tsmo(jpt_trd_npc,jp_sal)     &
+               &                 + tsmo(jpt_trd_dmp,jp_sal)                           + tsmo(jpt_trd_qns,jp_sal) )   &
+               &                / tvolt
+
+
+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)
@@ -630,25 +654,28 @@
 
 
-         IF(lwp) THEN
             WRITE (numout,*)
             WRITE (numout,*)
             WRITE (numout,9420) kt
-            WRITE (numout,9421) ( t2(jpicpt_xad)+t2(jpicpt_yad) )/ tvolt, ( s2(jpicpt_xad)+s2(jpicpt_yad) )/ tvolt
-            WRITE (numout,9422)   t2(jpicpt_zad) / tvolt, s2(jpicpt_zad) / tvolt
-            WRITE (numout,9423)   t2(jpicpt_ldf) / tvolt, s2(jpicpt_ldf) / tvolt
-            WRITE (numout,9424)   t2(jpicpt_zdf) / tvolt, s2(jpicpt_zdf) / tvolt
-            WRITE (numout,9425)   t2(jpicpt_npc) / tvolt, s2(jpicpt_npc) / tvolt
-            WRITE (numout,9426)   t2(jpicpt_dmp) / tvolt, s2(jpicpt_dmp) / tvolt
-            WRITE (numout,9427)   t2(jpicpt_qsr) / tvolt
-            WRITE (numout,9428)   t2(jpicpt_nsr) / tvolt, s2(jpicpt_nsr) / tvolt
+            WRITE (numout,9421) ( ts2(jpt_trd_xad,jp_tem) + ts2(jpt_trd_yad,jp_tem) ) / tvolt,   &
+               &                ( ts2(jpt_trd_xad,jp_sal) + ts2(jpt_trd_yad,jp_sal) ) / tvolt
+            WRITE (numout,9422)   ts2(jpt_trd_zad,jp_tem) / tvolt, ts2(jpt_trd_zad,jp_sal) / tvolt
+            WRITE (numout,9423)   ts2(jpt_trd_ldf,jp_tem) / tvolt, ts2(jpt_trd_ldf,jp_sal) / tvolt
+            WRITE (numout,9424)   ts2(jpt_trd_zdf,jp_tem) / tvolt, ts2(jpt_trd_zdf,jp_sal) / tvolt
+            WRITE (numout,9425)   ts2(jpt_trd_npc,jp_tem) / tvolt, ts2(jpt_trd_npc,jp_sal) / tvolt
+            WRITE (numout,9426)   ts2(jpt_trd_dmp,jp_tem) / tvolt, ts2(jpt_trd_dmp,jp_sal) / tvolt
+            WRITE (numout,9427)   ts2(jpt_trd_qsr,jp_tem) / tvolt
+            WRITE (numout,9428)   ts2(jpt_trd_qns,jp_tem) / tvolt, ts2(jpt_trd_qns,jp_sal) / tvolt
             WRITE (numout,9429)
-            WRITE (numout,9430) (  t2(jpicpt_xad) + t2(jpicpt_yad) + t2(jpicpt_zad) + t2(jpicpt_ldf) + t2(jpicpt_zdf)   &
-            &                    + t2(jpicpt_npc) + t2(jpicpt_dmp) + t2(jpicpt_qsr) + t2(jpicpt_nsr) ) / tvolt,   &
-            &                   (  s2(jpicpt_xad) + s2(jpicpt_yad) + s2(jpicpt_zad) + s2(jpicpt_ldf) + s2(jpicpt_zdf)   &
-            &                    + s2(jpicpt_npc) + s2(jpicpt_dmp)                  + s2(jpicpt_nsr) ) / tvolt
-         ENDIF
-
-9420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
-            '               salinity', /, ' ===============================')
+            WRITE (numout,9430) (  ts2(jpt_trd_xad,jp_tem) + ts2(jpt_trd_yad,jp_tem) + ts2(jpt_trd_zad,jp_tem)     &
+               &                 + ts2(jpt_trd_ldf,jp_tem) + ts2(jpt_trd_zdf,jp_tem) + ts2(jpt_trd_npc,jp_tem)     &
+               &                 + ts2(jpt_trd_dmp,jp_tem) + ts2(jpt_trd_qsr,jp_tem) + ts2(jpt_trd_qns,jp_tem) )   &
+               &                / tvolt,   &
+               &                (  ts2(jpt_trd_xad,jp_sal) + ts2(jpt_trd_yad,jp_sal) + ts2(jpt_trd_zad,jp_sal)     &
+               &                 + ts2(jpt_trd_ldf,jp_sal) + ts2(jpt_trd_zdf,jp_sal) + ts2(jpt_trd_npc,jp_sal)     &
+               &                 + ts2(jpt_trd_dmp,jp_sal)                          + ts2(jpt_trd_qns,jp_sal) )   &
+               &                / tvolt
+
+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)
@@ -663,24 +690,22 @@
 
 
-         IF(lwp) THEN
             WRITE (numout,*)
             WRITE (numout,*)
             WRITE (numout,9440) kt
-            WRITE (numout,9441) ( tmo(jpicpt_xad)+tmo(jpicpt_yad)+tmo(jpicpt_zad) )/tvolt,   &
-            &                   ( smo(jpicpt_xad)+smo(jpicpt_yad)+smo(jpicpt_zad) )/tvolt
-            WRITE (numout,9442)   tmo(jpicpt_zl1)/tvolt,  smo(jpicpt_zl1)/tvolt
-            WRITE (numout,9443)   tmo(jpicpt_ldf)/tvolt,  smo(jpicpt_ldf)/tvolt
-            WRITE (numout,9444)   tmo(jpicpt_zdf)/tvolt,  smo(jpicpt_zdf)/tvolt
-            WRITE (numout,9445)   tmo(jpicpt_npc)/tvolt,  smo(jpicpt_npc)/tvolt
-            WRITE (numout,9446) ( t2(jpicpt_xad)+t2(jpicpt_yad)+t2(jpicpt_zad) )/tvolt,    &
-            &                   ( s2(jpicpt_xad)+s2(jpicpt_yad)+s2(jpicpt_zad) )/tvolt
-            WRITE (numout,9447)   t2(jpicpt_ldf)/tvolt,   s2(jpicpt_ldf)/tvolt
-            WRITE (numout,9448)   t2(jpicpt_zdf)/tvolt,   s2(jpicpt_zdf)/tvolt
-            WRITE (numout,9449)   t2(jpicpt_npc)/tvolt,   s2(jpicpt_npc)/tvolt
-         ENDIF
-
-9440     FORMAT(' tracer consistency at it= ',i6,   &
-            ' :         temperature','                salinity',/,   &
-            ' ==================================')
+            WRITE (numout,9441) ( tsmo(jpt_trd_xad,jp_tem)+tsmo(jpt_trd_yad,jp_tem)+tsmo(jpt_trd_zad,jp_tem) )/tvolt,  &
+            &                   ( tsmo(jpt_trd_xad,jp_sal)+tsmo(jpt_trd_yad,jp_sal)+tsmo(jpt_trd_zad,jp_sal) )/tvolt
+            WRITE (numout,9442)   tsmo(jpt_trd_zl1,jp_tem)/tvolt,  tsmo(jpt_trd_zl1,jp_sal)/tvolt
+            WRITE (numout,9443)   tsmo(jpt_trd_ldf,jp_tem)/tvolt,  tsmo(jpt_trd_ldf,jp_sal)/tvolt
+            WRITE (numout,9444)   tsmo(jpt_trd_zdf,jp_tem)/tvolt,  tsmo(jpt_trd_zdf,jp_sal)/tvolt
+            WRITE (numout,9445)   tsmo(jpt_trd_npc,jp_tem)/tvolt,  tsmo(jpt_trd_npc,jp_sal)/tvolt
+            WRITE (numout,9446) ( ts2(jpt_trd_xad,jp_tem)+ts2(jpt_trd_yad,jp_tem)+ts2(jpt_trd_zad,jp_tem) )/tvolt,    &
+            &                   ( ts2(jpt_trd_xad,jp_sal)+ts2(jpt_trd_yad,jp_sal)+ts2(jpt_trd_zad,jp_sal) )/tvolt
+            WRITE (numout,9447)   ts2(jpt_trd_ldf,jp_tem)/tvolt,   ts2(jpt_trd_ldf,jp_sal)/tvolt
+            WRITE (numout,9448)   ts2(jpt_trd_zdf,jp_tem)/tvolt,   ts2(jpt_trd_zdf,jp_sal)/tvolt
+            WRITE (numout,9449)   ts2(jpt_trd_npc,jp_tem)/tvolt,   ts2(jpt_trd_npc,jp_sal)/tvolt
+         ENDIF
+
+9440     FORMAT(' tracer consistency at it= ',i6, ' :         temperature','                salinity',   &
+                /, ' ==================================')
 9441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
 9442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
Index: branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdicp_oce.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdicp_oce.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdicp_oce.F90	(revision 790)
@@ -14,38 +14,4 @@
    PUBLIC
 
-   !! * Shared module variables
-#if  defined key_trdtra   &&   defined key_trddyn    ||   defined key_esopa
-   LOGICAL, PARAMETER ::   lk_trdtra = .TRUE.    !: tracers  trend flag
-   LOGICAL, PARAMETER ::   lk_trddyn = .TRUE.    !: momentum trend flag
-#elif  defined key_trdtra
-   LOGICAL, PARAMETER ::   lk_trdtra = .TRUE.    !: tracers  trend flag
-   LOGICAL, PARAMETER ::   lk_trddyn = .FALSE.   !: momentum trend flag
-#elif  defined key_trddyn
-   LOGICAL, PARAMETER ::   lk_trdtra = .FALSE.   !: tracers  trend flag
-   LOGICAL, PARAMETER ::   lk_trddyn = .TRUE.    !: momentum trend flag
-#else
-   LOGICAL, PARAMETER ::   lk_trdtra = .FALSE.   !: tracers  trend flag
-   LOGICAL, PARAMETER ::   lk_trddyn = .FALSE.   !: momentum trend flag
-#endif
-
-   !! Tracers trends diagnostics parameters
-   !!---------------------------------------------------------------------
-   INTEGER, PARAMETER ::   &  !: => tracer trends indexes <=
-        jpicpt_xad =  1,   &  !: x- horizontal advection
-        jpicpt_yad =  2,   &  !: y- horizontal advection
-        jpicpt_zad =  3,   &  !: z- vertical   advection
-        jpicpt_ldf =  4,   &  !: lateral       diffusion
-        jpicpt_zdf =  5,   &  !: vertical diffusion (Kz)
-        jpicpt_bbc =  6,   &  !: Bottom Boundary Condition (geoth. flux) 
-        jpicpt_bbl =  7,   &  !: Bottom Boundary Layer (diffusive/convective)
-        jpicpt_npc =  8,   &  !: static instability mixing
-        jpicpt_dmp =  9,   &  !: damping
-        jpicpt_qsr = 10,   &  !: penetrative solar radiation
-        jpicpt_nsr = 11,   &  !: non solar radiation
-        jpicpt_zl1 = 12       !: first level vertical flux
-
-   INTEGER, PARAMETER ::   &  !: => Total tracer trends indexes <=
-        jptot_tra = 12        !: change it when adding/removing one indice above
-   
    !! Momentum trends diagnostics parameters
    !!---------------------------------------------------------------------
@@ -66,25 +32,4 @@
         jptot_dyn = 11        !: change it when adding/removing one indice above
    
-#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
-   REAL(wp) ::   tvolu        !: volume of the whole ocean computed at u-points
-   REAL(wp) ::   tvolv        !: volume of the whole ocean computed at v-points
-
-   !! Tracers trends diagnostics variables
-   !!---------------------------------------------------------------------
-   REAL(wp), DIMENSION(jptot_tra) ::   tmo, smo         !: tracers trends average 
-   REAL(wp), DIMENSION(jptot_tra) ::   t2, s2           !: tracers square trends average 
-   
-   !! Momentum trends diagnostics variables
-   !!---------------------------------------------------------------------
-   REAL(wp), DIMENSION(jptot_dyn) ::   umo, vmo         !: momentum trends average 
-   REAL(wp), DIMENSION(jptot_dyn) ::   hke              !: momentum square trends average 
-   REAL(wp) ::   rpktrd   !: potential to kinetic energy conversion
-   REAL(wp) ::   peke     !: conversion potential energy - kinetic energy trend
-
-#endif
    !!----------------------------------------------------------------------
    !!  OPA 9.0 , LOCEAN-IPSL (2005) 
Index: branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdmod.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 790)
@@ -9,5 +9,7 @@
 #if  defined key_trdtra || defined key_trddyn || defined key_trdmld || defined key_trdvor || defined key_esopa
    !!----------------------------------------------------------------------
-   !!   trd_mod          : Call the trend to be computed
+   !!   trd_tra          : active tracer trend manager
+   !!   trd_tra_adv      : pre-treatment of the tracer advection trends
+   !!   trd_mod          : momentum trend manager
    !!   trd_mod_init     : Initialization step
    !!----------------------------------------------------------------------
@@ -29,5 +31,7 @@
    REAL(wp) ::   r2dt          ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0
 
-   PUBLIC trd_mod              ! called by all dynXX or traXX modules
+   PUBLIC trd_tra              ! called by all traXXX.F90 modules
+   PUBLIC trd_tra_adv          ! called by all traadv_XXX.F90 modules
+   PUBLIC trd_mod              ! called by all dynXXX.F90 modules
    PUBLIC trd_mod_init         ! called by opa.F90 module
 
@@ -42,4 +46,221 @@
 
 CONTAINS
+
+   SUBROUTINE trd_tra( kt, ktra, ktrd, ctype, ptrd2d, ptrd3d, cnbpas )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mod  ***
+      !! 
+      !! ** Purpose : Dispatch all trends computation, e.g. vorticity, mld or 
+      !!              integral constraints
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in   )                                   ::   kt      ! time step
+      INTEGER         , INTENT(in   )                                   ::   ktra    ! tracer index
+      INTEGER         , INTENT(in   )                                   ::   ktrd    ! tracer trend index
+      CHARACTER(len=3), INTENT(in   )                                   ::   ctype   ! tracers type 'TRA' or 'TRC'
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj)    , OPTIONAL ::   ptrd2d  ! Temperature or U trend 
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk), OPTIONAL ::   ptrd3d  ! Temperature or U trend 
+      CHARACTER(len=3), INTENT(in   )                        , OPTIONAL ::   cnbpas  ! number of passage
+      !!
+      CHARACTER(len=3) ::   ccpas                                  ! number of passage
+!     REAL(wp), DIMENSION(jpi,jpj) ::   z2dx, z2dy                 ! workspace arrays
+      !!----------------------------------------------------------------------
+
+      ccpas = 'fst'                            ! Control of optional arguments
+      IF( PRESENT(cnbpas) )  ccpas = cnbpas
+
+      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restarting with Euler time stepping)
+      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog)
+      ENDIF
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! I. Integral Constraints Properties for momentum and/or tracers trends
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      IF( ( mod(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend) )   THEN
+         !
+         IF( lk_trdtra .AND. ctype == 'TRA' )   THEN       ! active tracer trends
+            IF( PRESENT(ptrd2d) ) THEN   ;   CALL trd_icp( ptrd2d, ktra, ktrd, ctype, clpas=ccpas )
+            ELSE                         ;   CALL trd_icp( ptrd3d, ktra, ktrd, ctype, clpas=ccpas )
+            ENDIF
+         ENDIF
+!           SELECT CASE ( ktrd )
+!           CASE ( jptra_trd_ldf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_ldf, ctype )   ! lateral diff
+!           CASE ( jptra_trd_zdf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_zdf, ctype )   ! vertical diff (Kz)
+!           CASE ( jptra_trd_bbc )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_bbc, ctype )   ! bottom boundary cond
+!           CASE ( jptra_trd_bbl )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_bbl, ctype )   ! bottom boundary layer
+!           CASE ( jptra_trd_npc )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_npc, ctype )   ! static instability mixing
+!           CASE ( jptra_trd_dmp )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_dmp, ctype )   ! damping
+!           CASE ( jptra_trd_qsr )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_qsr, ctype )   ! penetrative solar radiat.
+!           CASE ( jptra_trd_nsr )   
+!              z2dx(:,:) = ptrdx(:,:,1)   ;   z2dy(:,:) = ptrdy(:,:,1)
+!              CALL trd_icp( z2dx, z2dy, jpicpt_nsr, ctype )                               ! non solar radiation
+!           CASE ( jptra_trd_xad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_xad, ctype )   ! x- horiz adv
+!           CASE ( jptra_trd_yad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_yad, ctype )   ! y- horiz adv
+!           CASE ( jptra_trd_zad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_zad, ctype, clpas=ccpas )    ! z- adv 
+!              CALL trd_icp( ptrdx, ptrdy, jpicpt_zad, ctype, clpas=ccpas )   
+!              ! compute the surface flux condition wn(:,:,1)*tn(:,:,1)
+!              z2dx(:,:) = wn(:,:,1)*tn(:,:,1)/fse3t(:,:,1)
+!              z2dy(:,:) = wn(:,:,1)*sn(:,:,1)/fse3t(:,:,1)
+!              CALL trd_icp( z2dx , z2dy , jpicpt_zl1, ctype )                             ! 1st z- vertical adv 
+!           END SELECT
+!        END IF
+         !
+      END IF
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! III. Mixed layer trends for active tracers
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      IF( lk_trdmld .AND. ctype == 'TRA' )   THEN
+         
+         !-----------------------------------------------------------------------------------------------
+         ! W.A.R.N.I.N.G :
+         ! jptra_trd_ldf : called by traldf.F90
+         !                 at this stage we store:
+         !                  - the lateral geopotential diffusion (here, lateral = horizontal)
+         !                  - and the iso-neutral diffusion if activated 
+         ! jptra_trd_zdf : called by trazdf.F90
+         !                 * in case of iso-neutral diffusion we store the vertical diffusion component in the 
+         !                   lateral trend including the K_z contrib, which will be removed later (see trd_mld)
+         !-----------------------------------------------------------------------------------------------
+
+!        SELECT CASE ( ktrd )
+!        CASE ( jpt_trd_zdf )
+!           IF( ln_traldf_iso )   THEN
+!              CALL trd_mld_zint( ptrdx, ptrdy, jpmld_ldf, '3D' )   ! vertical diffusion (K_z)
+!           ELSE
+!              CALL trd_mld_zint( ptrdx, ptrdy, jpmld_zdf, '3D' )   ! vertical diffusion (K_z)
+!           ENDIF
+!        CASE ( jpt_trd_nsr )
+!           ptrdx(:,:,2:jpk) = 0.e0   ;   ptrdy(:,:,2:jpk) = 0.e0
+!           CALL trd_mld_zint( ptrdx, ptrdy, jpmld_for, '2D' )                             ! air-sea : non solar flux
+!        CASE ( jpt_trd_bbc )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_bbc, '3D' )   ! bbc (geothermal flux)
+!        CASE DEFAULT
+!                                     CALL trd_mld_zint( ptrdx, ptrdy, ktra     , '3D' )   ! other 3D trends
+!        END SELECT
+
+      ENDIF
+
+   END SUBROUTINE trd_tra
+
+
+   SUBROUTINE trd_tra_adv( kt, ktra, ktrd, ctype, pf, pun, ptn, cnbpas )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_mod  ***
+      !! 
+      !! ** Purpose :   transformed the i-advective flux into the i-advective trends
+      !! ** Method  :   i-advective trends = un. di[T] = di[fi] - tn di[un]
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in   )                         ::   kt      ! time step
+      INTEGER         , INTENT(in   )                         ::   ktra    ! tracer index
+      INTEGER         , INTENT(in   )                         ::   ktrd    ! tracer trend index
+      CHARACTER(len=3), INTENT(in   )                         ::   ctype   ! tracers type 'TRA' or 'TRC'
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pf      ! advective flux in one direction
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun     ! now velocity  in one direction
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptn     ! now or before tracer 
+      CHARACTER(len=3), INTENT(in   )                        , OPTIONAL ::   cnbpas  ! number of passage
+      !!
+      INTEGER                          ::   ji, jj, jk      ! dummy loop indices
+      CHARACTER(len=3) ::   ccpas                           ! number of passage
+      REAL(wp)                         ::   zbtr, z_hdivn   ! 
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdt           ! 3D workspace 
+      !!----------------------------------------------------------------------
+
+      ccpas = 'fst'                            ! Control of optional arguments
+      IF( PRESENT(cnbpas) )  ccpas = cnbpas
+
+      ztrdt(:,:,:) = 0.e0
+      !
+      IF( ccpas == 'fst' ) THEN      ! first treatment : remove the divergence
+         SELECT CASE( ktrd ) 
+         CASE( jpt_trd_xad )      ! i-advective trend
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+#if defined key_zco
+                     zbtr    = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) )
+                     z_hdivn = (  e2u(ji  ,jj) * pun(ji  ,jj,jk)                              &
+                        &       - e2u(ji-1,jj) * pun(ji-1,jj,jk) )
+#else
+                     zbtr    = 1.e0/ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                     z_hdivn = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
+                        &       - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) )
+#endif
+                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji-1,jj,jk)  -  ptn(ji,jj,jk) * z_hdivn  )
+                  END DO
+               END DO
+            END DO
+            !
+         CASE( jpt_trd_yad )      ! j-advective trend
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+#if defined key_zco
+                     zbtr    = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) )
+                     z_hdivn = (  e1v(ji,jj  ) * pun(ji,jj  ,jk)                              &
+                        &       - e1v(ji,jj-1) * pun(ji,jj-1,jk) )
+#else
+                     zbtr    = 1.e0/ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                     z_hdivn = (  e1v(ji  ,jj) * fse3v(ji,jj  ,jk) * pun(ji,jj  ,jk)          &
+                        &       - e1v(ji-1,jj) * fse3v(ji,jj-1,jk) * pun(ji,jj-1,jk) )
+#endif
+                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji,jj-1,jk)  -  ptn(ji,jj,jk) * z_hdivn  )
+                  END DO
+               END DO
+            END DO
+            !
+         CASE( jpt_trd_zad )      ! z-advective trend
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zbtr    = 1.e0 / fse3t(ji,jj,jk)
+                     z_hdivn = pun(ji,jj,jk) - pun(ji,jj,jk+1)
+                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji,jj,jk+1)  -  ptn(ji,jj,jk) * z_hdivn  )
+                  END DO
+               END DO
+            END DO
+            !
+         END SELECT
+         !
+      ELSE                ! second call : just compute the trend   (TVD scheme)
+         SELECT CASE( ktrd ) 
+         CASE( jpt_trd_xad )      ! i-advective trend
+            DO jk = 1, jpkm1      
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zbtr    = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji-1,jj,jk)  )
+                  END DO
+               END DO
+            END DO
+            !
+         CASE( jpt_trd_yad )      ! j-advective trend
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zbtr    = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji,jj-1,jk)  )
+                  END DO
+               END DO
+            END DO
+            !
+         CASE( jpt_trd_zad )      ! z-advective trend
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     zbtr    = 1.e0 / fse3t(ji,jj,jk)
+                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji,jj,jk+1)  )
+                  END DO
+               END DO
+            END DO
+            !
+         END SELECT
+         !
+      ENDIF
+      !
+      CALL trd_tra( kt, ktra, ktrd, ctype, ptrd3d=ztrdt)       ! trend diagnostics
+      !
+   END SUBROUTINE trd_tra_adv
+
 
    SUBROUTINE trd_mod( ptrdx, ptrdy, ktrd, ctype, kt, cnbpas )
@@ -50,10 +271,10 @@
       !!              integral constraints
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt                                ! time step
-      INTEGER, INTENT( in ) ::   ktrd                              ! tracer trend index
-      CHARACTER(len=3), INTENT( in ) ::   ctype                    ! momentum or tracers trends type 'DYN'/'TRA'
-      CHARACTER(len=3), INTENT( in ), OPTIONAL ::   cnbpas         ! number of passage
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   ptrdx ! Temperature or U trend 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   ptrdy ! Salinity    or V trend
+      INTEGER         , INTENT(in   )                                   ::   kt      ! time step
+      INTEGER         , INTENT(in   )                                   ::   ktrd    ! tracer trend index
+      CHARACTER(len=3), INTENT(in   )                                   ::   ctype   ! momentum trends type ='DYN'
+      CHARACTER(len=3), INTENT(in   )                        , OPTIONAL ::   cnbpas  ! number of passage
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk)           ::   ptrdx   ! Temperature or U trend 
+      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk)           ::   ptrdy   ! Salinity    or V trend
       !!
       INTEGER ::   ji, ikbu, ikbum1
@@ -68,9 +289,8 @@
       z2dx(:,:)   = 0.e0   ;   z2dy(:,:)   = 0.e0                  ! initialization of workspace arrays
 
-      ! Control of optional arguments
-      ccpas = 'fst'
+      ccpas = 'fst'                            ! Control of optional arguments
       IF( PRESENT(cnbpas) )  ccpas = cnbpas
 
-      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restarting with Euler time stepping)
+      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (start with Euler time stepping)
       ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog)
       ENDIF
@@ -82,27 +302,4 @@
       IF( ( mod(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend) )   THEN
          !
-         IF( lk_trdtra .AND. ctype == 'TRA' )   THEN       ! active tracer trends
-            SELECT CASE ( ktrd )
-            CASE ( jptra_trd_ldf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_ldf, ctype )   ! lateral diff
-            CASE ( jptra_trd_zdf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_zdf, ctype )   ! vertical diff (Kz)
-            CASE ( jptra_trd_bbc )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_bbc, ctype )   ! bottom boundary cond
-            CASE ( jptra_trd_bbl )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_bbl, ctype )   ! bottom boundary layer
-            CASE ( jptra_trd_npc )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_npc, ctype )   ! static instability mixing
-            CASE ( jptra_trd_dmp )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_dmp, ctype )   ! damping
-            CASE ( jptra_trd_qsr )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_qsr, ctype )   ! penetrative solar radiat.
-            CASE ( jptra_trd_nsr )   
-               z2dx(:,:) = ptrdx(:,:,1)   ;   z2dy(:,:) = ptrdy(:,:,1)
-               CALL trd_icp( z2dx, z2dy, jpicpt_nsr, ctype )                               ! non solar radiation
-            CASE ( jptra_trd_xad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_xad, ctype )   ! x- horiz adv
-            CASE ( jptra_trd_yad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_yad, ctype )   ! y- horiz adv
-            CASE ( jptra_trd_zad )                                                         ! z- vertical adv 
-               CALL trd_icp( ptrdx, ptrdy, jpicpt_zad, ctype, clpas=ccpas )   
-               ! compute the surface flux condition wn(:,:,1)*tn(:,:,1)
-               z2dx(:,:) = wn(:,:,1)*tn(:,:,1)/fse3t(:,:,1)
-               z2dy(:,:) = wn(:,:,1)*sn(:,:,1)/fse3t(:,:,1)
-               CALL trd_icp( z2dx , z2dy , jpicpt_zl1, ctype )                             ! 1st z- vertical adv 
-            END SELECT
-         END IF
-
          IF( lk_trddyn .AND. ctype == 'DYN' )   THEN       ! momentum trends 
             !
@@ -200,46 +397,5 @@
          !
       ENDIF
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! III. Mixed layer trends for active tracers
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-      IF( lk_trdmld .AND. ctype == 'TRA' )   THEN
-         
-         !-----------------------------------------------------------------------------------------------
-         ! W.A.R.N.I.N.G :
-         ! jptra_trd_ldf : called by traldf.F90
-         !                 at this stage we store:
-         !                  - the lateral geopotential diffusion (here, lateral = horizontal)
-         !                  - and the iso-neutral diffusion if activated 
-         ! jptra_trd_zdf : called by trazdf.F90
-         !                 * in case of iso-neutral diffusion we store the vertical diffusion component in the 
-         !                   lateral trend including the K_z contrib, which will be removed later (see trd_mld)
-         !-----------------------------------------------------------------------------------------------
-
-         SELECT CASE ( ktrd )
-         CASE ( jptra_trd_xad )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_xad, '3D' )   ! merid. advection
-         CASE ( jptra_trd_yad )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_yad, '3D' )   ! zonal  advection
-         CASE ( jptra_trd_zad )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_zad, '3D' )   ! vertical advection
-         CASE ( jptra_trd_ldf )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_ldf, '3D' )   ! lateral diffusive
-         CASE ( jptra_trd_bbl )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_bbl, '3D' )   ! bottom boundary layer
-         CASE ( jptra_trd_zdf )
-            IF( ln_traldf_iso )   THEN
-               CALL trd_mld_zint( ptrdx, ptrdy, jpmld_ldf, '3D' )   ! vertical diffusion (K_z)
-            ELSE
-               CALL trd_mld_zint( ptrdx, ptrdy, jpmld_zdf, '3D' )   ! vertical diffusion (K_z)
-            ENDIF
-         CASE ( jptra_trd_dmp )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_dmp, '3D' )   ! internal 3D restoring (tradmp)
-         CASE ( jptra_trd_qsr )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_for, '3D' )   ! air-sea : penetrative sol radiat
-         CASE ( jptra_trd_nsr )
-            ptrdx(:,:,2:jpk) = 0.e0   ;   ptrdy(:,:,2:jpk) = 0.e0
-            CALL trd_mld_zint( ptrdx, ptrdy, jpmld_for, '2D' )                             ! air-sea : non penetr sol radiat
-         CASE ( jptra_trd_bbc )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_bbc, '3D' )   ! bottom bound cond (geoth flux)
-         CASE ( jptra_trd_atf )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_atf, '3D' )   ! asselin numerical
-         CASE ( jptra_trd_npc )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_npc, '3D' )   ! non penetr convect adjustment
-         END SELECT
-
-      ENDIF
-
+      !
    END SUBROUTINE trd_mod
 
@@ -254,10 +410,27 @@
 
 CONTAINS
+
+   SUBROUTINE trd_tra( kt, ktra, ktrd, ctype, ptrd2d, ptrd3d, cnbpas )
+      INTEGER ::   kt, ktra, ktrd
+      CHARACTER(len=3) ::   ctype   ! tracers type 'TRA' or 'TRC'
+      REAL, DIMENSION(:,:)  , OPTIONAL ::   ptrd2d  ! Temperature or U trend 
+      REAL, DIMENSION(:,:,:), OPTIONAL ::   ptrd3d  ! Temperature or U trend 
+      CHARACTER(len=3)      , OPTIONAL ::   cnbpas  ! number of passage
+      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', kt, ktra, ktrd, ctype,   &
+         &                                                               ptrd2d(1,1), ptrd3d(1,1,1), cnbpas
+   END SUBROUTINE trd_tra
+
+   SUBROUTINE trd_tra_adv( kt, ktra, ktrd, ctype, pf, pun, ptn, cnbpas )
+      INTEGER ::   kt, ktra, ktrd  
+      CHARACTER(len=3) ::   ctype
+      REAL, DIMENSION(:,:,:) ::   pf, pun, ptn
+      CHARACTER(len=3), OPTIONAL ::   cnbpas
+      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', kt, ktra, ktrd, ctype, pf(1,1,1),   &
+         &                                                               pun(1,1,1), ptn(1,1,1), cnbpas
+   END SUBROUTINE trd_tra_adv
+
    SUBROUTINE trd_mod(ptrd3dx, ptrd3dy, ktrd , ctype, kt, cnbpas)   ! 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
+      REAL, DIMENSION(:,:,:) ::   ptrd3dx, ptrd3dy    
+      INTEGER ::   ktrd, kt
       CHARACTER(len=3), INTENT( in ) ::  ctype                     ! momentum or tracers trends type
       CHARACTER(len=3), INTENT( in ), OPTIONAL ::   cnbpas         ! number of passage
Index: branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdmod_oce.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdmod_oce.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdmod_oce.F90	(revision 790)
@@ -32,18 +32,36 @@
    LOGICAL , PUBLIC ::   l_trddyn = .FALSE.             !: momentum trend flag
 # endif
+#if  defined key_trdtra   &&   defined key_trddyn    ||   defined key_esopa
+   LOGICAL, PARAMETER ::   lk_trdtra = .TRUE.    !: tracers  trend flag
+   LOGICAL, PARAMETER ::   lk_trddyn = .TRUE.    !: momentum trend flag
+#elif  defined key_trdtra
+   LOGICAL, PARAMETER ::   lk_trdtra = .TRUE.    !: tracers  trend flag
+   LOGICAL, PARAMETER ::   lk_trddyn = .FALSE.   !: momentum trend flag
+#elif  defined key_trddyn
+   LOGICAL, PARAMETER ::   lk_trdtra = .FALSE.   !: tracers  trend flag
+   LOGICAL, PARAMETER ::   lk_trddyn = .TRUE.    !: momentum trend flag
+#else
+   LOGICAL, PARAMETER ::   lk_trdtra = .FALSE.   !: tracers  trend flag
+   LOGICAL, PARAMETER ::   lk_trddyn = .FALSE.   !: momentum trend flag
+#endif
+
 
    !!* Active tracers trends indexes
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_xad =  1   !: x- horizontal advection
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_yad =  2   !: y- horizontal advection
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_zad =  3   !: z- vertical   advection
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_ldf =  4   !: lateral       diffusion
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_zdf =  5   !: vertical diffusion (Kz)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_bbc =  6   !: Bottom Boundary Condition (geoth. flux) 
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_bbl =  7   !: Bottom Boundary Layer (diffusive/convective)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_npc =  8   !: static instability mixing
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_dmp =  9   !: damping
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_qsr = 10   !: penetrative solar radiation
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_nsr = 11   !: non solar radiation
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_atf = 12   !: Asselin correction
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd = 13        !: change it when adding/removing one indice above
+
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_xad =  1   !: x- horizontal advection
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_yad =  2   !: y- horizontal advection
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_zad =  3   !: z- vertical   advection
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_ldf =  4   !: lateral       diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_zdf =  5   !: vertical diffusion (Kz)
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_bbc =  6   !: Bottom Boundary Condition (geoth. flux) 
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_bbl =  7   !: Bottom Boundary Layer (diffusive/convective)
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_npc =  8   !: static instability mixing
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_dmp =  9   !: damping
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_qsr = 10   !: penetrative solar radiation
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_qns = 11   !: non solar heat flux
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_atf = 12   !: Asselin correction
+   INTEGER, PUBLIC, PARAMETER ::   jpt_trd_zl1 = 13   !: first level vertical flux (lk_vvl=F)
+
    
    !!* Momentum trends indexes
Index: branches/dev_001_GM/NEMO/OPA_SRC/TRD/trends_manager.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/TRD/trends_manager.F90	(revision 790)
+++ branches/dev_001_GM/NEMO/OPA_SRC/TRD/trends_manager.F90	(revision 790)
@@ -0,0 +1,14 @@
+MODULE trends_manager
+   !!======================================================================
+   !!                       ***  MODULE  trends_manager  ***
+   !! Ocean diagnostics:  ocean tracers and dynamic trends
+   !!=====================================================================
+   !! History :  9.0  !  04-08  (C. Talandier) Original code
+   !!                 !  05-04  (C. Deltel)    Add Asselin trend in the ML budget
+   !!----------------------------------------------------------------------
+#if  defined key_trdtra || defined key_trddyn || defined key_trdmld || defined key_trdvor || defined key_esopa
+
+#endif
+
+   !!======================================================================
+END MODULE trends_manager
Index: branches/dev_001_GM/NEMO/OPA_SRC/ZDF/zdfkpp.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 790)
@@ -1291,5 +1291,5 @@
          ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:)
 !!bug gm jpttdzdf ==> jpttkpp
-         CALL trd_mod(ztrdt, ztrds, jptra_trd_zdf, 'TRA', kt)
+         CALL trd_mod(ztrdt, ztrds, jpt_trd_zdf, 'TRA', kt)
       ENDIF
 
Index: branches/dev_001_GM/NEMO/OPA_SRC/oce.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/oce.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/oce.F90	(revision 790)
@@ -4,17 +4,16 @@
    !! Ocean        :  dynamics and active tracers defined in memory 
    !!======================================================================
-   !! History :
-   !!   8.5  !  02-11  (G. Madec)  F90: Free form and module
-   !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
+   !! History :   1.0  !  02-11  (G. Madec)  F90: Free form and module
+   !!             2.0  !  05-11  (V. Garnier) Surface pressure gradient organization
    !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! NEMO/OPA 2.0 , LOCEAN-IPSL (2007) 
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_oce      ! ocean parameters
 
    IMPLICIT NONE
    PRIVATE
+
 
    !! Physics and algorithm flags
@@ -23,23 +22,20 @@
    INTEGER, PUBLIC ::   nn_dynhpg_rst   = 0        !: add dynhpg implicit variables in restart ot not
 
-   !! dynamics and tracer fields
-   !! --------------------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
-      ! before !  now      !  after  !      ! the after trends becomes the fields
-      ! fields !  fields   !  trends !      ! only in dyn(tra)_zdf and dyn(tra)_nxt
-      ub       ,  un       ,  ua     ,   &  !: i-horizontal velocity (m/s)
-      vb       ,  vn       ,  va     ,   &  !: j-horizontal velocity (m/s)
-                  wn       ,             &  !: vertical velocity (m/s)
-      rotb     ,  rotn     ,             &  !: relative vorticity (1/s)
-      hdivb    ,  hdivn    ,             &  !: horizontal divergence (1/s)
-      tb       ,  tn       ,  ta     ,   &  !: potential temperature (celcius)
-      sb       ,  sn       ,  sa            !: salinity (psu)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
-      rhd ,                              &  !: in situ density anomalie rhd=(rho-rau0)/rau0 (no units)
-      rhop,                              &  !: potential volumic mass (kg/m3)
-      rn2                                   !: brunt-vaisala frequency (1/s2)
+   !! dynamics and tracer fields              ! before ! now    ! after  !   the after trends becomes the fields
+   !! --------------------------              ! fields ! fields ! trends !   only in dyn(tra)_zdf and dyn(tra)_nxt
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   ub    ,  un    ,  ua    !: i-horizontal velocity (m/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   vb    ,  vn    ,  va    !: j-horizontal velocity (m/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::            wn             !: vertical velocity (m/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rotb  ,  rotn           !: relative vorticity (1/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   hdivb ,  hdivn          !: horizontal divergence (1/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   tb    ,  tn    ,  ta    !: potential temperature (celcius)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   sb    ,  sn    ,  sa    !: salinity (psu)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rhd                     !: in situ density anomalie rhd=(rho-rau0)/rau0
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rhop                    !: potential volumic mass (kg/m3)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   rn2                     !: brunt-vaisala frequency (1/s2)
 
       !! advection scheme choice
       !! -----------------------
+!!gm  only for trcadv...   be be suppressed
       CHARACTER(len=3), PUBLIC  ::   l_adv   !: 'ce2' centre scheme used
          !                                   !: 'tvd' TVD scheme used
@@ -49,27 +45,24 @@
    !! surface pressure gradient
    !! -------------------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      spgu, spgv             !: horizontal surface pressure gradient
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   spgu, spgv                  !: horizontal surface pressure gradient
 
    !! interpolated gradient (only used in zps case)
    !! ---------------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      gtu, gsu, gru,      &  !: t-, s- and rd horizontal gradient at u- and 
-      gtv, gsv, grv          !: v-points at bottom ocean level 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gtu, gsu, gru               !: t-, s- and rd horizontal gradient at u- and 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gtv, gsv, grv               !: v-points at bottom ocean level 
 
    !! free surface
    !! ------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      sshb , sshn ,        &  !: before, now sea surface height (meters)
-      sshu , sshv ,        &  !: sea surface height at u- and v- point
-      sshbb, ssha             !: before before sea surface height at t-point
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshb , sshn                 !: before, now sea surface height (meters)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshu , sshv                 !: sea surface height at u- and v- point
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sshbb, ssha                 !: before before sea surface height at t-point
 
 #if defined key_dynspg_rl   ||   defined key_esopa
    !! rigid-lid formulation
    !! ---------------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      bsfb, bsfn,         &  !: before, now barotropic streamfunction (m3/s)
-      bsfd                   !: now trend of barotropic streamfunction (m3/s2)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bsfb, bsfn                 !: before, now barotropic streamfunction [m3/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bsfd                       !: now barotropic streamfunction trend   [m3/s2]
 #endif
+
    !!----------------------------------------------------------------------
 END MODULE oce
Index: branches/dev_001_GM/NEMO/OPA_SRC/par_oce.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/par_oce.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/par_oce.F90	(revision 790)
@@ -9,5 +9,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 2.4 , LOCEAN-IPSL (2007) 
-   !! $Id:$ 
+   !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -98,4 +98,12 @@
 #endif
 
+   !!---------------------------------------------------------------------
+   !! Temperature/salinity indices
+   !!---------------------------------------------------------------------
+
+   INTEGER, PARAMETER ::   jp_tem = 1   !: index of temperature in the tsb, tsns &, tsa arrays
+   INTEGER, PARAMETER ::   jp_sal = 2   !: index of salinity    in the tsb, tsns &, tsa arrays
+
+
    !!======================================================================
 END MODULE par_oce
Index: branches/dev_001_GM/NEMO/OPA_SRC/prtctl.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/prtctl.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/prtctl.F90	(revision 790)
@@ -1,7 +1,9 @@
 MODULE prtctl
-   !!==============================================================================
+   !!======================================================================
    !!                       ***  MODULE prtctl   ***
    !! Ocean system   : print all SUM trends for each processor domain
-   !!==============================================================================
+   !!======================================================================
+   !! History :  2.0  !  05-07  (C. Talandier) original code
+   !!----------------------------------------------------------------------
    USE dom_oce          ! ocean space and time domain variables
    USE in_out_manager   ! I/O manager
@@ -11,8 +13,6 @@
    PRIVATE
 
-   !! * Module declaration
    INTEGER, DIMENSION(:), ALLOCATABLE :: numid
-   INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::   &  !:
-      nlditl , nldjtl ,   &  !: first, last indoor index for each i-domain
+   INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::   nlditl , nldjtl ,   &  !: first, last indoor index for each i-domain
       nleitl , nlejtl ,   &  !: first, last indoor index for each j-domain
       nimpptl, njmpptl,   &  !: i-, j-indexes for each processor
@@ -20,18 +20,17 @@
       ibonitl, ibonjtl
 
-   REAL(wp), DIMENSION(:), ALLOCATABLE ::   &  !:
-      t_ctll , s_ctll ,   &  !: previous trend values
+   REAL(wp), DIMENSION(:), ALLOCATABLE ::   t_ctll , s_ctll ,   &  !: previous trend values
       u_ctll , v_ctll
 
    INTEGER ::   ktime        !: time step
 
-   !! * Routine accessibility
    PUBLIC prt_ctl         ! called by all subroutines
    PUBLIC prt_ctl_info    ! called by all subroutines
    PUBLIC prt_ctl_init    ! called by opa.F90
+
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    !!----------------------------------------------------------------------
 
@@ -74,9 +73,5 @@
       !!                    kdim    : k- direction for 3D arrays 
       !!                    clinfo3 : additional information 
-      !!
-      !! History :
-      !!   9.0  !  05-07  (C. Talandier) original code
-      !!----------------------------------------------------------------------
-      !! * Arguments
+      !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:)  , INTENT(in), OPTIONAL :: tab2d_1
       REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: tab3d_1
@@ -90,9 +85,8 @@
       INTEGER                   , INTENT(in), OPTIONAL :: kdim
       CHARACTER (len=*)         , INTENT(in), OPTIONAL :: clinfo3
-
-      !! * Local declarations
-      INTEGER :: overlap, jn, sind, eind, kdir,j_id
+      !!
       CHARACTER (len=15) :: cl2
-      REAL(wp) :: zsum1, zsum2, zvctl1, zvctl2
+      INTEGER  ::   overlap, jn, sind, eind, kdir,j_id
+      REAL(wp) ::   zsum1, zsum2, zvctl1, zvctl2
       REAL(wp), DIMENSION(jpi,jpj) :: ztab2d_1, ztab2d_2
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmask1, zmask2, ztab3d_1, ztab3d_2
@@ -125,10 +119,8 @@
       IF( PRESENT(mask2)   )  zmask2  (:,:,:)= mask2  (:,:,:)
 
-      IF( lk_mpp )   THEN
-         ! processor number
+      IF( lk_mpp )   THEN      ! processor number
          sind = narea
          eind = narea
-      ELSE
-         ! processors total number
+      ELSE                     ! processors total number
          sind = 1
          eind = ijsplt
Index: branches/dev_001_GM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/dev_001_GM/NEMO/OPA_SRC/step.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/OPA_SRC/step.F90	(revision 790)
@@ -4,20 +4,21 @@
    !! Time-stepping    : manager of the ocean, tracer and ice time stepping
    !!======================================================================
-   !! History :        !  91-03  ()  Original code
-   !!                  !  91-11  (G. Madec)
-   !!                  !  92-06  (M. Imbard)  add a first output record
-   !!                  !  96-04  (G. Madec)  introduction of dynspg
-   !!                  !  96-04  (M.A. Foujols)  introduction of passive tracer
-   !!             8.0  !  97-06  (G. Madec)  new architecture of call
-   !!             8.2  !  97-06  (G. Madec, M. Imbard, G. Roullet)  free surface
-   !!             8.2  !  99-02  (G. Madec, N. Grima)  hpg implicit
-   !!             8.2  !  00-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
-   !!             9.0  !  02-06  (G. Madec)  free form, suppress macro-tasking
-   !!             " "  !  04-08  (C. Talandier) New trends organization
-   !!             " "  !  05-01  (C. Ethe) Add the KPP closure scheme
-   !!             " "  !  05-11  (V. Garnier) Surface pressure gradient organization
-   !!             " "  !  05-11  (G. Madec)  Reorganisation of tra and dyn calls
-   !!             " "  !  06-01  (L. Debreu, C. Mazauric)  Agrif implementation
-   !!             " "  !  06-07  (S. Masson)  restart using iom
+   !! History :        !  1991-03  ()  Original code
+   !!                  !  1991-11  (G. Madec)
+   !!                  !  1992-06  (M. Imbard)  add a first output record
+   !!                  !  1996-04  (G. Madec)  introduction of dynspg
+   !!                  !  1996-04  (M.A. Foujols)  introduction of passive tracer
+   !!             8.0  !  1997-06  (G. Madec)  new architecture of call
+   !!             8.2  !  1997-06  (G. Madec, M. Imbard, G. Roullet)  free surface
+   !!             8.2  !  1999-02  (G. Madec, N. Grima)  hpg implicit
+   !!             8.2  !  2000-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
+   !!   NEMO      1.0  !  2002-06  (G. Madec)  free form, suppress macro-tasking
+   !!              -   !  2004-08  (C. Talandier) New trends organization
+   !!              -   !  2005-01  (C. Ethe) Add the KPP closure scheme
+   !!              -   !  2005-11  (V. Garnier) Surface pressure gradient organization
+   !!              -   !  2005-11  (G. Madec)  Reorganisation of tra and dyn calls
+   !!             2.0  !  2006-01  (L. Debreu, C. Mazauric)  Agrif implementation
+   !!              -   !  2006-07  (S. Masson)  restart using iom
+   !!             2.4  !  2008-01  (G. Madec) Merge TRA-TRC
    !!----------------------------------------------------------------------
    !!   stp            : OPA system time-stepping
@@ -323,5 +324,8 @@
       IF( lk_trabbl_adv  )   CALL tra_bbl_adv( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
 
-      IF( lk_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends
+      IF( lk_tradmp      )   THEN                           ! internal damping trends
+                             CALL tra_dmp    ( kstp, 'TRA', jp_tem, tb, ta )      ! temperature
+                             CALL tra_dmp    ( kstp, 'TRA', jp_sal, sb, sa )      ! salinity
+      ENDIF
 
                              CALL tra_adv    ( kstp )       ! horizontal & vertical advection
Index: branches/dev_001_GM/NEMO/TOP_SRC/LOBSTER/trcsed.F90
===================================================================
--- branches/dev_001_GM/NEMO/TOP_SRC/LOBSTER/trcsed.F90	(revision 786)
+++ branches/dev_001_GM/NEMO/TOP_SRC/LOBSTER/trcsed.F90	(revision 790)
@@ -28,5 +28,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
-   !! $Id:$ 
+   !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -90,4 +90,5 @@
 # endif
 # if defined key_trc_diaadd
+!!gm  bug introduced: no more mask below jpkb ! 
                trc2d(ji,jj,8) = trc2d(ji,jj,8) + ztra * fse3t(ji,jj,jk) * 86400.
 # endif
@@ -97,5 +98,5 @@
 
 #if defined key_trc_diabio
-      CALL lbc_lnk (trbio(:,:,1,8), 'T', 1. )    ! Lateral boundary conditions on trcbio
+      CALL lbc_lnk( trbio(:,:,1,8), 'T', 1. )    ! Lateral boundary conditions on trcbio
 #endif
 #if defined key_trc_diaadd
