Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml	(revision 3317)
@@ -118,27 +118,30 @@
     <field id="Mmy_v"          description="Mm current barotrope along j-axis harmonic imaginary part " unit="m/s"        />
     <!-- variables available with key_trdtra -->
-    <field id="ttrd_xad"     description="temperature-trend: i-advection"           unit="degC/s" axis_ref="deptht" />
-    <field id="strd_xad"     description="salinity   -trend: i-advection"           unit="psu/s"  axis_ref="deptht" />
-    <field id="ttrd_yad"     description="temperature-trend: j-advection"           unit="degC/s" axis_ref="deptht" />
-    <field id="strd_yad"     description="salinity   -trend: j-advection"           unit="psu/s"  axis_ref="deptht" />
-    <field id="ttrd_zad"     description="temperature-trend: k-advection"           unit="degC/s" axis_ref="deptht" />
-    <field id="strd_zad"     description="salinity   -trend: k-advection"           unit="psu/s"  axis_ref="deptht" />
-    <field id="ttrd_sad"     description="temperature-trend: surface adv. (no-vvl)" unit="degC/s"                   />
-    <field id="strd_sad"     description="salinity   -trend: surface adv. (no-vvl)" unit="psu/s"                    />
-    <field id="ttrd_ldf"     description="temperature-trend: lateral  diffusion"    unit="degC/s" axis_ref="deptht" />
-    <field id="strd_ldf"     description="salinity   -trend: lateral  diffusion"    unit="psu/s"  axis_ref="deptht" />
-    <field id="ttrd_zdf"     description="temperature-trend: vertical diffusion"    unit="degC/s" axis_ref="deptht" />
-    <field id="strd_zdf"     description="salinity   -trend: vertical diffusion"    unit="psu/s"  axis_ref="deptht" />
-    <field id="ttrd_zdfp"    description="temperature-trend: pure vert. diffusion"  unit="degC/s" axis_ref="deptht" />
-    <field id="strd_zdfp"    description="salinity   -trend: pure vert. diffusion"  unit="psu/s"  axis_ref="deptht" />
-    <field id="ttrd_dmp"     description="temperature-trend: interior restoring"    unit="degC/s" axis_ref="deptht" />
-    <field id="strd_dmp"     description="salinity   -trend: interior restoring"    unit="psu/s"  axis_ref="deptht" />
-    <field id="ttrd_bbl"     description="temperature-trend: bottom boundary layer" unit="degC/s" axis_ref="deptht" />
-    <field id="strd_bbl"     description="salinity   -trend: bottom boundary layer" unit="psu/s"  axis_ref="deptht" />
-    <field id="ttrd_npc"     description="temperature-trend: non-penetrative conv." unit="degC/s" axis_ref="deptht" />
-    <field id="strd_npc"     description="salinity   -trend: non-penetrative conv." unit="psu/s"  axis_ref="deptht" />
-    <field id="ttrd_qns"     description="temperature-trend: non-solar heating"     unit="degC/s"                   />
-    <field id="ttrd_qsr"     description="temperature-trend: solar penetr. heating" unit="degC/s"  axis_ref="deptht" />
-    <field id="ttrd_bbc"     description="temperature-trend: geothermal heating"    unit="degC/s"  axis_ref="deptht" />
+    <field id="ttrd_xad"     description="temperature-trend: i-advection"              unit="degC/s" axis_ref="deptht" />
+    <field id="strd_xad"     description="salinity   -trend: i-advection"              unit="psu/s"  axis_ref="deptht" />
+    <field id="ttrd_yad"     description="temperature-trend: j-advection"              unit="degC/s" axis_ref="deptht" />
+    <field id="strd_yad"     description="salinity   -trend: j-advection"              unit="psu/s"  axis_ref="deptht" />
+    <field id="ttrd_zad"     description="temperature-trend: k-advection"              unit="degC/s" axis_ref="deptht" />
+    <field id="strd_zad"     description="salinity   -trend: k-advection"              unit="psu/s"  axis_ref="deptht" />
+    <field id="ttrd_sad"     description="temperature-trend: surface adv. (no-vvl)"    unit="degC/s"                   />
+    <field id="strd_sad"     description="salinity   -trend: surface adv. (no-vvl)"    unit="psu/s"                    />
+    <field id="ttrd_ldf"     description="temperature-trend: lateral  diffusion"       unit="degC/s" axis_ref="deptht" />
+    <field id="strd_ldf"     description="salinity   -trend: lateral  diffusion"       unit="psu/s"  axis_ref="deptht" />
+    <field id="ttrd_zdf"     description="temperature-trend: vertical diffusion"       unit="degC/s" axis_ref="deptht" />
+    <field id="strd_zdf"     description="salinity   -trend: vertical diffusion"       unit="psu/s"  axis_ref="deptht" />
+    < ln_traldf_iso=T only (iso-neutral diffusion) >
+    <field id="ttrd_zdfp"    description="temperature-trend: pure vert. diffusion"     unit="degC/s" axis_ref="deptht" />
+    <field id="strd_zdfp"    description="salinity   -trend: pure vert. diffusion"     unit="psu/s"  axis_ref="deptht" />
+    <>
+    <field id="ttrd_dmp"     description="temperature-trend: interior restoring"       unit="degC/s" axis_ref="deptht" />
+    <field id="strd_dmp"     description="salinity   -trend: interior restoring"       unit="psu/s"  axis_ref="deptht" />
+    <field id="ttrd_bbl"     description="temperature-trend: bottom boundary layer"    unit="degC/s" axis_ref="deptht" />
+    <field id="strd_bbl"     description="salinity   -trend: bottom boundary layer"    unit="psu/s"  axis_ref="deptht" />
+    <field id="ttrd_npc"     description="temperature-trend: non-penetrative conv."    unit="degC/s" axis_ref="deptht" />
+    <field id="strd_npc"     description="salinity   -trend: non-penetrative conv."    unit="psu/s"  axis_ref="deptht" />
+    <field id="ttrd_qns"     description="temperature-trend: non-solar flux + runoff"  unit="degC/s" axis_ref="deptht" />
+    <field id="strd_cdt"     description="salinity   -trend: C/D term       + runoff"  unit="degC/s" axis_ref="deptht" />
+    <field id="ttrd_qsr"     description="temperature-trend: solar penetr. heating"    unit="degC/s" axis_ref="deptht" />
+    <field id="ttrd_bbc"     description="temperature-trend: geothermal heating"       unit="degC/s" axis_ref="deptht" />
       </group>
 
Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 3317)
@@ -57,7 +57,10 @@
       IF( nn_timing == 1 )  CALL timing_start('dyn_bfr')
       !
+!!gm issue: better to put the logical in step to control the call of zdf_bfr
+!!          ==> change the logical from ln_bfrimp to ln_bfr_exp !!
       IF( .NOT.ln_bfrimp) THEN     ! only for explicit bottom friction form
                                     ! implicit bfr is implemented in dynzdf_imp
 
+!!gm bug : time step is only rdt (not 2 rdt if euler start !)
         zm1_2dt = - 1._wp / ( 2._wp * rdt )
 
Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 3317)
@@ -31,4 +31,6 @@
    USE dynadv          ! dynamics: vector invariant versus flux form
    USE domvvl          ! variable volume
+   USE trdmod         ! ocean dynamics trends 
+   USE trdmod_oce     ! ocean variables trends
    USE obc_oce         ! ocean open boundary conditions
    USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
@@ -292,8 +294,7 @@
 
       IF( ln_3D_trd_d ) THEN             ! 3D output: asselin filter trends on momentum
-        zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
-        zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
-        CALL iom_put( "utrd_atf", zua )        ! asselin filter trends on momentum 
-        CALL iom_put( "vtrd_atf", zva )
+         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
+         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
+         CALL trd_mod( zua, zva, jpdyn_trd_atf, 'DYN', kt )
       ENDIF
 
Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 3317)
@@ -94,18 +94,18 @@
       IF( ln_bfrimp ) THEN
 # if defined key_vectopt_loop
-      DO jj = 1, 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+         DO jj = 1, 1
+            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
 # else
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
 # endif
-            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
-            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-            zavmu(ji,jj) = avmu(ji,jj,ikbu+1)
-            zavmv(ji,jj) = avmv(ji,jj,ikbv+1)
-            avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
-            avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
-         END DO
-      END DO
+               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
+               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
+               zavmu(ji,jj) = avmu(ji,jj,ikbu+1)
+               zavmv(ji,jj) = avmv(ji,jj,ikbv+1)
+               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
+               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
+            END DO
+         END DO
       ENDIF
 
@@ -284,16 +284,16 @@
       IF( ln_bfrimp ) THEN
 # if defined key_vectopt_loop
-      DO jj = 1, 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+         DO jj = 1, 1
+            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
 # else
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
 # endif
-            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
-            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-            avmu(ji,jj,ikbu+1) = zavmu(ji,jj)
-            avmv(ji,jj,ikbv+1) = zavmv(ji,jj)
-         END DO
-      END DO
+               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
+               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
+               avmu(ji,jj,ikbu+1) = zavmu(ji,jj)
+               avmv(ji,jj,ikbv+1) = zavmv(ji,jj)
+            END DO
+         END DO
       ENDIF
       !
Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 3317)
@@ -47,5 +47,5 @@
    PUBLIC   dtacof_zoom  ! routine called by in both tradmp.F90 and trcdmp.F90
 
-   !                                !!* Namelist namtra_dmp : T & S newtonian damping *
+   !                                         !!* Namelist namtra_dmp : T & S newtonian damping *
    LOGICAL, PUBLIC ::   ln_tradmp = .TRUE.    !: internal damping flag
    INTEGER         ::   nn_hdmp   =   -1      ! = 0/-1/'latitude' for damping over T and S
@@ -111,4 +111,5 @@
       !
       CALL wrk_alloc( jpi, jpj, jpk, jpts,  zts_dta )
+      !
       !                           !==   input T-S data at kt   ==!
       CALL dta_tsd( kt, zts_dta )            ! read and interpolates T-S data at kt
Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 3317)
@@ -29,4 +29,5 @@
    USE zdf_oce         ! ???
    USE domvvl          ! variable volume
+   USE ldftra_oce      ! ocean active tracers lateral physics
    USE dynspg_oce      ! surface     pressure gradient variables
    USE dynhpg          ! hydrostatic pressure gradient 
@@ -132,4 +133,8 @@
          ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
          ztrds(:,:,:) = tsn(:,:,:,jp_sal)
+         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend 
+            CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_zdfp, ztrdt )
+            CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_zdfp, ztrds )
+         ENDIF
       ENDIF
 
@@ -168,6 +173,5 @@
          &                       tab3d_2=tsn(:,:,:,jp_sal), clinfo2=       ' Sn: ', mask2=tmask )
       !
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_nxt')
+      IF( nn_timing == 1 )   CALL timing_stop('tra_nxt')
       !
    END SUBROUTINE tra_nxt
Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 3317)
@@ -104,5 +104,5 @@
       !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated
       !!                with the tracer surface boundary condition 
-      !!              - save the trend it in ttrd ('key_trdtra')
+      !!              - send trends to trdtra module (l_trdtra=T)
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
@@ -123,5 +123,5 @@
       zsrau = 1. / rau0             ! initialization
 
-      IF( l_trdtra )   THEN                    !* Save ta and sa trends
+      IF( l_trdtra ) THEN                    !* Save ta and sa trends
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
@@ -145,14 +145,14 @@
               & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
             IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file'
-            zfact = 0.5e0
+            zfact = 0.5_wp
             CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend
             CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend
          ELSE                                         ! No restart or restart not found: Euler forward time stepping
-            zfact = 1.e0
+            zfact = 1._wp
             sbc_tsc_b(:,:,:) = 0.e0
          ENDIF
       ELSE                                         ! Swap of forcing fields
          !                                         ! ----------------------
-         zfact = 0.5e0
+         zfact = 0.5_wp
          sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
       ENDIF
@@ -206,12 +206,11 @@
       !----------------------------------------
       !
-      zfact = 0.5e0
 
       ! Effect on (t,s) due to river runoff (dilution effect automatically applied via vertical tracer advection) 
       IF( ln_rnf ) THEN  
+         zfact = 0.5_wp
          DO jj = 2, jpj 
             DO ji = fs_2, fs_jpim1
-               zdep = 1. / h_rnf(ji,jj)
-               zdep = zfact * zdep  
+               zdep = zfact / h_rnf(ji,jj)
                IF ( rnf(ji,jj) /= 0._wp ) THEN
                   DO jk = 1, nk_rnf(ji,jj)
@@ -225,6 +224,7 @@
          END DO  
       ENDIF  
-!!gm  It should be useless
+!!gm  It should be useless  ==>> to be suppressed !
       CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )    ;    CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
+!!gm : end of useless
 
       IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 3317)
@@ -7,20 +7,21 @@
    !!            3.5  !  2012-02 (G. Madec)  add 3D tracer zdf trend output using iom
    !!----------------------------------------------------------------------
-#if  defined key_trdtra   ||   defined key_trddyn   ||   defined key_esopa
+
    !!----------------------------------------------------------------------
-   !!   'key_trdtra'  or                  active tracers trends diagnostics
-   !!   'key_trddyn'                            momentum trends diagnostics
-   !!----------------------------------------------------------------------
-   !!   trd_icp          : compute the basin averaged properties for tra/dyn 
-   !!   trd_dwr          : print dynmaic trends in ocean.output file
-   !!   trd_twr          : print tracers trends in ocean.output file
-   !!   trd_icp_init     : initialization step
+   !!   trd_budget     : domain averaged budget of trends (including kinetic energy and T^2 trends)
+   !!   trd_icp        : compute the basin averaged properties for tra/dyn 
+   !!   trd_dwr        : print dynmaic trends in ocean.output file
+   !!   trd_twr        : print tracers trends in ocean.output file
+   !!   trd_icp_init   : initialization step
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE phycst          ! physical constants
    USE trdmod_oce      ! ocean variables trends
    USE ldftra_oce      ! ocean active tracers: lateral physics
    USE ldfdyn_oce      ! ocean dynamics: lateral physics
    USE zdf_oce         ! ocean vertical physics
+   USE zdfbfr          ! bottom friction
    USE zdfddm          ! ocean vertical physics: double diffusion
    USE eosbn2          ! equation of state
@@ -34,12 +35,21 @@
    PRIVATE
 
-   INTERFACE trd_icp
-      MODULE PROCEDURE trd_2d, trd_3d
-   END INTERFACE
-
-   PUBLIC   trd_icp       ! called by trdmod.F90
+   PUBLIC   trd_budget    ! called by trdmod.F90
    PUBLIC   trd_dwr       ! called by step.F90
    PUBLIC   trd_twr       ! called by step.F90
    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
+   REAL(wp) ::   rpktrd   ! potential to kinetic energy conversion
+   REAL(wp) ::   peke     ! conversion potential energy - kinetic energy trend
+
+   !                     !!! domain averaged trends
+   REAL(wp), DIMENSION(jptot_tra) ::   tmo, smo   ! temperature and salinity trends 
+   REAL(wp), DIMENSION(jptot_tra) ::   t2 , s2    ! T^2 and S^2 trends 
+   REAL(wp), DIMENSION(jptot_dyn) ::   umo, vmo   ! momentum trends 
+   REAL(wp), DIMENSION(jptot_dyn) ::   hke        ! kinetic energy trends (u^2+v^2) 
 
    !! * Substitutions
@@ -54,146 +64,124 @@
 CONTAINS
 
-   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype )
+   SUBROUTINE trd_budget( ptrdx, ptrdy, ktrd, ctype, kt )
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_2d  ***
+      !!                  ***  ROUTINE trd_budget  ***
       !! 
-      !! ** Purpose :   compute and print the domain averaged properties of tracers 
-      !!              and/or momentum equations at each nn_trd time step.
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ptrd2dx   ! Temperature or U trend 
-      REAL(wp), DIMENSION(:,:), 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
+      !! ** Purpose : integral constraint diagnostics for momentum and/or tracer trends
+      !!              
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend 
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
+      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
+      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type 'DYN'/'TRA'
+      INTEGER                   , INTENT(in   ) ::   kt      ! time step
       !!
-      INTEGER ::   ji, jj   ! loop indices
-      !!----------------------------------------------------------------------
-
-      SELECT CASE( ctype )    !==  Mask trends  ==!
-      !
-      CASE( 'DYN' )                    ! Momentum
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               ptrd2dx(ji,jj) = ptrd2dx(ji,jj) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,1)
-               ptrd2dy(ji,jj) = ptrd2dy(ji,jj) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,1)
+      INTEGER ::   ji, jj, jk      ! dummy loop indices
+      INTEGER ::   ikbu, ikbv      ! local integers
+      REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace 
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
+
+      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
+         !
+      	SELECT CASE( ctype )
+      	!
+      	CASE( 'TRA' )          !==  Tracers (T & S)  ==!
+      	   DO jk = 1, jpkm1       ! global sum of mask volume trend and trend*T (including interior mask)
+      	      DO jj = 1, jpj
+      	         DO ji = 1, jpi        
+      	            zvm = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+      	            zvt = ptrdx(ji,jj,jk) * zvm
+      	            zvs = ptrdy(ji,jj,jk) * zvm
+         	         tmo(ktrd) = tmo(ktrd) + zvt   
+         	         smo(ktrd) = smo(ktrd) + zvs
+         	         t2 (ktrd) = t2(ktrd)  + zvt * tsn(ji,jj,jk,jp_tem)
+         	         s2 (ktrd) = s2(ktrd)  + zvs * tsn(ji,jj,jk,jp_sal)
+                  END DO
+               END DO
             END DO
-         END DO
-         ptrd2dx(jpi, : ) = 0._wp      ;      ptrd2dy(jpi, : ) = 0._wp
-         ptrd2dx( : ,jpj) = 0._wp      ;      ptrd2dy( : ,jpj) = 0._wp
-         !
-      CASE( 'TRA' )                    ! Tracers
-         ptrd2dx(:,:) = ptrd2dx(:,:) * tmask_i(:,:)
-         ptrd2dy(:,:) = ptrd2dy(:,:) * tmask_i(:,:)
-         !
-      END SELECT
-      
-      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum trends  ==!
-      !
-      CASE( 'DYN' )                    ! Momentum
-         umo(ktrd) = 0._wp
-         vmo(ktrd) = 0._wp
-         !
-         SELECT CASE( ktrd )
-         CASE( jpdyn_trd_swf )         ! surface forcing
-            umo(ktrd) = SUM( ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1) )
-            vmo(ktrd) = SUM( ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1) )
+            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
+            IF( .NOT.lk_vvl .AND. ktrd == jptra_trd_zad ) THEN  
+               tmo(jptra_trd_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * e1e2t(:,:) )
+               smo(jptra_trd_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
+               t2 (jptra_trd_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) * e1e2t(:,:)  )
+               s2 (jptra_trd_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
+            ENDIF
+            !
+            IF( ktrd == jptra_trd_atf ) THEN     ! last trend (asselin time filter)
+               ! 
+               CALL trd_twr( kt )   ! print the results in ocean.output
+               !                
+	            tmo(:) = 0._wp       ! prepare the next time step (domain averaged array reset to zero)
+   	         smo(:) = 0._wp
+               t2 (:) = 0._wp
+               s2 (:) = 0._wp
+               !
+            ENDIF
+            !
+         CASE( 'DYN' )          !==  Momentum and KE  ==!        
+            DO jk = 1, jpkm1
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1
+                     zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
+                        &                  * e1u    (ji  ,jj  ) * e2u    (ji,jj) * fse3u(ji,jj,jk)
+                     zvs = ptrdy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
+                        &                  * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
+                     umo(ktrd) = umo(ktrd) + zvt
+                     vmo(ktrd) = vmo(ktrd) + zvs
+                     hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * zvt + vn(ji,jj,jk) * zvs
+                  END DO
+               END DO
+            END DO
+            !                 
+            IF( ktrd == jpdyn_trd_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
+               z1_2rau0 = 0.5_wp / rau0
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1
+                     zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
+                        &                       * z1_2rau0 * e1u    (ji  ,jj  ) * e2u    (ji,jj)
+                     zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
+                        &                       * z1_2rau0 * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
+                     umo(jpdyn_trd_tau) = umo(jpdyn_trd_tau) + zvt
+                     vmo(jpdyn_trd_tau) = vmo(jpdyn_trd_tau) + zvs
+                     hke(jpdyn_trd_tau) = hke(jpdyn_trd_tau) + un(ji,jj,1) * zvt + vn(ji,jj,1) * zvs
+                  END DO
+               END DO
+            ENDIF
+            !                       
+            IF( ktrd == jpdyn_trd_atf ) THEN     ! last trend (asselin time filter)
+               !
+               IF( ln_bfrimp ) THEN                   ! implicit bfr case: compute separately the bottom friction 
+                  z1_2rau0 = 0.5_wp / rau0
+                  DO jj = 1, jpjm1
+                     DO ji = 1, jpim1
+                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
+                        ikbv = mbkv(ji,jj)
+                        zvt = bfrua(ji,jj) * un(ji,jj,ikbu) * e1u(ji,jj) * e2v(ji,jj)
+                        zvs = bfrva(ji,jj) * vn(ji,jj,ikbv) * e1v(ji,jj) * e2v(ji,jj)
+                        umo(jpdyn_trd_bfr) = umo(jpdyn_trd_bfr) + zvt
+                        vmo(jpdyn_trd_bfr) = vmo(jpdyn_trd_bfr) + zvs
+                        hke(jpdyn_trd_bfr) = hke(jpdyn_trd_bfr) + un(ji,jj,ikbu) * zvt + vn(ji,jj,ikbv) * zvs
+                     END DO
+                  END DO
+               ENDIF
+               ! 
+               CALL trd_dwr( kt )                     ! print the results in ocean.output
+               !                
+               umo(:) = 0._wp                         ! reset for the next time step
+               vmo(:) = 0._wp
+               hke(:) = 0._wp
+               !
+            ENDIF
+            !
          END SELECT
          !
-      CASE( 'TRA' )              ! Tracers
-         tmo(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
-         smo(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
-      END SELECT
-      
-      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum square trends  ==!   (now field)
-      !
-      CASE( 'DYN' )              ! Momentum
-         hke(ktrd) = SUM(   un(:,:,1) * ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1)   &
-            &             + vn(:,:,1) * ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1)   )
-         !
-      CASE( 'TRA' )              ! Tracers
-         t2(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) * tsn(:,:,1,jp_tem) )
-         s2(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) * tsn(:,:,1,jp_sal) )
-         !      
-      END SELECT
-      !
-   END SUBROUTINE trd_2d
-
-
-   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd, ctype )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_3d  ***
-      !! 
-      !! ** Purpose : verify the basin averaged properties of tracers and/or
-      !!              momentum equations at every time step frequency nn_trd.
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrd3dx   ! Temperature or U trend 
-      REAL(wp), DIMENSION(:,:,:), 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
-      !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      !!----------------------------------------------------------------------
-
-      SELECT CASE( ctype )    !==  Mask the trends  ==!
-      !
-      CASE( 'DYN' )              ! Momentum        
-         DO jk = 1, jpkm1
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
-                  ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         ptrd3dx(jpi, : ,:) = 0._wp      ;      ptrd3dy(jpi, : ,:) = 0._wp
-         ptrd3dx( : ,jpj,:) = 0._wp      ;      ptrd3dy( : ,jpj,:) = 0._wp
-         !
-      CASE( 'TRA' )              ! Tracers
-         DO jk = 1, jpkm1
-            ptrd3dx(:,:,jk) = ptrd3dx(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
-            ptrd3dy(:,:,jk) = ptrd3dy(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
-         END DO
-         !
-      END SELECT   
-
-      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum trends  ==!
-      !
-      CASE( 'DYN' )              ! Momentum
-         umo(ktrd) = 0._wp
-         vmo(ktrd) = 0._wp
-         DO jk = 1, jpkm1
-            umo(ktrd) = umo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk) )
-            vmo(ktrd) = vmo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk) )
-         END DO
-         !
-      CASE( 'TRA' )              ! Tracers
-         tmo(ktrd) = 0._wp
-         smo(ktrd) = 0._wp
-         DO jk = 1, jpkm1
-            tmo(ktrd) = tmo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
-            smo(ktrd) = smo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
-         END DO
-         !
-      END SELECT
-
-      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum square trends  ==!   (now field)
-      !
-      CASE( 'DYN' )              ! Momentum
-         hke(ktrd) = 0._wp
-         DO jk = 1, jpkm1
-            hke(ktrd) = hke(ktrd) + SUM(   un(:,:,jk) * ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk)   &
-               &                         + vn(:,:,jk) * ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk)   )
-         END DO
-         !
-      CASE( 'TRA' )              ! Tracers
-         t2(ktrd) = 0._wp
-         s2(ktrd) = 0._wp
-         DO jk = 1, jpkm1
-            t2(ktrd) = t2(ktrd) + SUM( ptrd3dx(:,:,jk) * tsn(:,:,jk,jp_tem) * e1e2t(:,:) * fse3t(:,:,jk) )
-            s2(ktrd) = s2(ktrd) + SUM( ptrd3dy(:,:,jk) * tsn(:,:,jk,jp_sal) * e1e2t(:,:) * fse3t(:,:,jk) )
-         END DO
-         !
-      END SELECT
-      !
-   END SUBROUTINE trd_3d
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
+      !
+   END SUBROUTINE trd_budget
 
 
@@ -262,6 +250,4 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)  ::  zkx, zky, zkz, zkepe  
       !!----------------------------------------------------------------------
-
-      IF( ln_3D_trd_t .AND. ln_3D_trd_d )   RETURN            ! do nothing if 3D output with IOM
 
       CALL wrk_alloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
@@ -333,24 +319,23 @@
             WRITE (numout,*)
             WRITE (numout,9500) kt
-            WRITE (numout,9501) umo(jpicpd_hpg) / tvolu, vmo(jpicpd_hpg) / tvolv
-            WRITE (numout,9502) umo(jpicpd_keg) / tvolu, vmo(jpicpd_keg) / tvolv
-            WRITE (numout,9503) umo(jpicpd_rvo) / tvolu, vmo(jpicpd_rvo) / tvolv
-            WRITE (numout,9504) umo(jpicpd_pvo) / tvolu, vmo(jpicpd_pvo) / tvolv
-            WRITE (numout,9505) umo(jpicpd_ldf) / tvolu, vmo(jpicpd_ldf) / tvolv
-            WRITE (numout,9506) umo(jpicpd_had) / tvolu, vmo(jpicpd_had) / tvolv
-            WRITE (numout,9507) umo(jpicpd_zad) / tvolu, vmo(jpicpd_zad) / tvolv
-            WRITE (numout,9508) umo(jpicpd_zdf) / tvolu, vmo(jpicpd_zdf) / tvolv
-            WRITE (numout,9509) umo(jpicpd_spg) / tvolu, vmo(jpicpd_spg) / tvolv
-            WRITE (numout,9510) umo(jpicpd_swf) / tvolu, vmo(jpicpd_swf) / tvolv
-            WRITE (numout,9511) umo(jpicpd_dat) / tvolu, vmo(jpicpd_dat) / tvolv
-            WRITE (numout,9512) umo(jpicpd_bfr) / tvolu, vmo(jpicpd_bfr) / tvolv
+            WRITE (numout,9501) umo(jpdyn_trd_hpg) / tvolu, vmo(jpdyn_trd_hpg) / tvolv
+            WRITE (numout,9509) umo(jpdyn_trd_spg) / tvolu, vmo(jpdyn_trd_spg) / tvolv
+            WRITE (numout,9502) umo(jpdyn_trd_keg) / tvolu, vmo(jpdyn_trd_keg) / tvolv
+            WRITE (numout,9503) umo(jpdyn_trd_rvo) / tvolu, vmo(jpdyn_trd_rvo) / tvolv
+            WRITE (numout,9504) umo(jpdyn_trd_pvo) / tvolu, vmo(jpdyn_trd_pvo) / tvolv
+            WRITE (numout,9507) umo(jpdyn_trd_zad) / tvolu, vmo(jpdyn_trd_zad) / tvolv
+            WRITE (numout,9505) umo(jpdyn_trd_ldf) / tvolu, vmo(jpdyn_trd_ldf) / tvolv
+            WRITE (numout,9508) umo(jpdyn_trd_zdf) / tvolu, vmo(jpdyn_trd_zdf) / tvolv
+            WRITE (numout,9510) umo(jpdyn_trd_tau) / tvolu, vmo(jpdyn_trd_tau) / tvolv
+            WRITE (numout,9511) umo(jpdyn_trd_bfr) / tvolu, vmo(jpdyn_trd_bfr) / tvolv
+            WRITE (numout,9512) umo(jpdyn_trd_atf) / tvolu, vmo(jpdyn_trd_atf) / tvolv
             WRITE (numout,9513)
             WRITE (numout,9514)                                                 &
-            &     (  umo(jpicpd_hpg) + umo(jpicpd_keg) + umo(jpicpd_rvo) + umo(jpicpd_pvo) + umo(jpicpd_ldf)   &
-            &      + umo(jpicpd_had) + umo(jpicpd_zad) + umo(jpicpd_zdf) + umo(jpicpd_spg) + umo(jpicpd_dat)   &
-            &      + umo(jpicpd_swf) + umo(jpicpd_bfr) ) / tvolu,   &
-            &     (  vmo(jpicpd_hpg) + vmo(jpicpd_keg) + vmo(jpicpd_rvo) + vmo(jpicpd_pvo) + vmo(jpicpd_ldf)   &
-            &      + vmo(jpicpd_had) + vmo(jpicpd_zad) + vmo(jpicpd_zdf) + vmo(jpicpd_spg) + vmo(jpicpd_dat)   &
-            &      + vmo(jpicpd_swf) + vmo(jpicpd_bfr) ) / tvolv
+            &     (  umo(jpdyn_trd_hpg) + umo(jpdyn_trd_spg) + umo(jpdyn_trd_keg) + umo(jpdyn_trd_rvo)   &
+            &      + umo(jpdyn_trd_pvo) + umo(jpdyn_trd_zad) + umo(jpdyn_trd_ldf) + umo(jpdyn_trd_zdf)   &
+            &      + umo(jpdyn_trd_tau) + umo(jpdyn_trd_bfr) + umo(jpdyn_trd_atf) ) / tvolu,   &
+            &     (  vmo(jpdyn_trd_hpg) + vmo(jpdyn_trd_spg) + vmo(jpdyn_trd_keg) + vmo(jpdyn_trd_rvo)   &
+            &      + vmo(jpdyn_trd_pvo) + vmo(jpdyn_trd_zad) + vmo(jpdyn_trd_ldf) + vmo(jpdyn_trd_zdf)   &
+            &      + vmo(jpdyn_trd_tau) + vmo(jpdyn_trd_bfr) + vmo(jpdyn_trd_atf) ) / tvolv
          ENDIF
 
@@ -361,11 +346,10 @@
  9504    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
  9505    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
- 9506    FORMAT(' horizontal advection       u= ', e20.13, '    v= ', e20.13)
  9507    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
  9508    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
  9509    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
- 9510    FORMAT(' surface wind forcing       u= ', e20.13, '    v= ', e20.13)
- 9511    FORMAT(' dampimg term               u= ', e20.13, '    v= ', e20.13)
- 9512    FORMAT(' bottom flux                u= ', e20.13, '    v= ', e20.13)
+ 9510    FORMAT(' surface wind stress        u= ', e20.13, '    v= ', e20.13)
+ 9511    FORMAT(' bottom friction            u= ', e20.13, '    v= ', e20.13)
+ 9512    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
  9513    FORMAT(' -----------------------------------------------------------------------------')
  9514    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
@@ -375,36 +359,34 @@
             WRITE (numout,*)
             WRITE (numout,9520) kt
-            WRITE (numout,9521) hke(jpicpd_hpg) / tvolt
-            WRITE (numout,9522) hke(jpicpd_keg) / tvolt
-            WRITE (numout,9523) hke(jpicpd_rvo) / tvolt
-            WRITE (numout,9524) hke(jpicpd_pvo) / tvolt
-            WRITE (numout,9525) hke(jpicpd_ldf) / tvolt
-            WRITE (numout,9526) hke(jpicpd_had) / tvolt
-            WRITE (numout,9527) hke(jpicpd_zad) / tvolt
-            WRITE (numout,9528) hke(jpicpd_zdf) / tvolt
-            WRITE (numout,9529) hke(jpicpd_spg) / tvolt
-            WRITE (numout,9530) hke(jpicpd_swf) / tvolt
-            WRITE (numout,9531) hke(jpicpd_dat) / tvolt
-            WRITE (numout,9532) hke(jpicpd_bfr) / tvolt
+            WRITE (numout,9521) hke(jpdyn_trd_hpg) / tvolt
+            WRITE (numout,9529) hke(jpdyn_trd_spg) / tvolt
+            WRITE (numout,9522) hke(jpdyn_trd_keg) / tvolt
+            WRITE (numout,9523) hke(jpdyn_trd_rvo) / tvolt
+            WRITE (numout,9524) hke(jpdyn_trd_pvo) / tvolt
+            WRITE (numout,9527) hke(jpdyn_trd_zad) / tvolt
+            WRITE (numout,9525) hke(jpdyn_trd_ldf) / tvolt
+            WRITE (numout,9528) hke(jpdyn_trd_zdf) / tvolt
+            WRITE (numout,9530) hke(jpdyn_trd_tau) / tvolt
+            WRITE (numout,9531) hke(jpdyn_trd_bfr) / tvolt
+            WRITE (numout,9532) hke(jpdyn_trd_atf) / tvolt
             WRITE (numout,9533)
             WRITE (numout,9534)   &
-            &     (  hke(jpicpd_hpg) + hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_pvo) + hke(jpicpd_ldf)   &
-            &      + hke(jpicpd_had) + hke(jpicpd_zad) + hke(jpicpd_zdf) + hke(jpicpd_spg) + hke(jpicpd_dat)   &
-            &      + hke(jpicpd_swf) + hke(jpicpd_bfr) ) / tvolt
+            &     (  hke(jpdyn_trd_hpg) + hke(jpdyn_trd_spg) + hke(jpdyn_trd_keg) + hke(jpdyn_trd_rvo)   &
+            &      + hke(jpdyn_trd_pvo) + hke(jpdyn_trd_zad) + hke(jpdyn_trd_ldf) + hke(jpdyn_trd_zdf)   &
+            &      + hke(jpdyn_trd_tau) + hke(jpdyn_trd_bfr) + hke(jpdyn_trd_atf)  ) / tvolt
          ENDIF
 
  9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
  9521    FORMAT(' pressure gradient         u2= ', e20.13)
+ 9529    FORMAT(' surface pressure gradient u2= ', e20.13)
  9522    FORMAT(' ke gradient               u2= ', e20.13)
  9523    FORMAT(' relative vorticity term   u2= ', e20.13)
  9524    FORMAT(' coriolis term             u2= ', e20.13)
+ 9527    FORMAT(' vertical advection        u2= ', e20.13)
  9525    FORMAT(' horizontal diffusion      u2= ', e20.13)
- 9526    FORMAT(' horizontal advection      u2= ', e20.13)
- 9527    FORMAT(' vertical advection        u2= ', e20.13)
  9528    FORMAT(' vertical diffusion        u2= ', e20.13)
- 9529    FORMAT(' surface pressure gradient u2= ', e20.13)
- 9530    FORMAT(' surface wind forcing      u2= ', e20.13)
- 9531    FORMAT(' dampimg term              u2= ', e20.13)
- 9532    FORMAT(' bottom flux               u2= ', e20.13)
+ 9530    FORMAT(' surface wind stress       u2= ', e20.13)
+ 9531    FORMAT(' bottom friction           u2= ', e20.13)
+ 9532    FORMAT(' Asselin time filter       u2= ', e20.13)
  9533    FORMAT(' --------------------------------------------------')
  9534    FORMAT(' total trend               u2= ', e20.13)
@@ -414,12 +396,12 @@
             WRITE (numout,*)
             WRITE (numout,9540) kt
-            WRITE (numout,9541) ( hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
-            WRITE (numout,9542) ( hke(jpicpd_keg) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
-            WRITE (numout,9543) ( hke(jpicpd_pvo) ) / tvolt
-            WRITE (numout,9544) ( hke(jpicpd_rvo) ) / tvolt
-            WRITE (numout,9545) ( hke(jpicpd_spg) ) / tvolt
-            WRITE (numout,9546) ( hke(jpicpd_ldf) ) / tvolt
-            WRITE (numout,9547) ( hke(jpicpd_zdf) ) / tvolt
-            WRITE (numout,9548) ( hke(jpicpd_hpg) ) / tvolt, rpktrd / tvolt
+            WRITE (numout,9541) ( hke(jpdyn_trd_keg) + hke(jpdyn_trd_rvo) + hke(jpdyn_trd_zad) ) / tvolt
+            WRITE (numout,9542) ( hke(jpdyn_trd_keg) + hke(jpdyn_trd_zad) ) / tvolt
+            WRITE (numout,9543) ( hke(jpdyn_trd_pvo) ) / tvolt
+            WRITE (numout,9544) ( hke(jpdyn_trd_rvo) ) / tvolt
+            WRITE (numout,9545) ( hke(jpdyn_trd_spg) ) / tvolt
+            WRITE (numout,9546) ( hke(jpdyn_trd_ldf) ) / tvolt
+            WRITE (numout,9547) ( hke(jpdyn_trd_zdf) ) / tvolt
+            WRITE (numout,9548) ( hke(jpdyn_trd_hpg) ) / tvolt, rpktrd / tvolt
             WRITE (numout,*)
             WRITE (numout,*)
@@ -427,11 +409,11 @@
 
  9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
- 9541    FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
- 9542    FORMAT(' 0 = ke gradient + horizontal + vertical advection         : ', e20.13)
- 9543    FORMAT(' 0 = coriolis term  (true if key_vorenergy or key_combined): ', e20.13)
- 9544    FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.)        : ', e20.13)
- 9545    FORMAT(' 0 = surface pressure gradient                             : ', e20.13)
- 9546    FORMAT(' 0 > horizontal diffusion                                  : ', e20.13)
- 9547    FORMAT(' 0 > vertical diffusion                                    : ', e20.13)
+ 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
+ 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
+ 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
+ 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
+ 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
+ 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
+ 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
  9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
          !
@@ -455,36 +437,5 @@
       !
       INTEGER  ::   jk   ! loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:)  ::   zwt, zws, ztrdt, ztrds   ! 3D workspace
-      !!----------------------------------------------------------------------
-
-
-      IF( ln_3D_trd_t ) THEN      ! 3D output: treat the vertical diffusion trends (if iso)
-         !
-         CALL wrk_alloc( jpi, jpj, jpk, zwt, zws, ztrdt, ztrds )
-         !
-         IF( ln_traldf_iso ) THEN      ! iso-neutral diffusion : re-compute the PURE vertical diffusive trend
-            !                                 !  zdf trends using now field (called after the swap)
-            zwt(:,:, 1 ) = 0._wp   ;   zws(:,:, 1 ) = 0._wp            ! vertical diffusive fluxes
-            zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp
-            DO jk = 2, jpk
-               zwt(:,:,jk) =   avt(:,:,jk) * ( tsn(:,:,jk-1,jp_tem) - tsn(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk)
-               zws(:,:,jk) = fsavs(:,:,jk) * ( tsn(:,:,jk-1,jp_sal) - tsn(:,:,jk,jp_sal) ) / fse3w(:,:,jk) * tmask(:,:,jk)
-            END DO
-            !
-            ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
-            DO jk = 1, jpkm1
-               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / fse3t(:,:,jk)
-               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / fse3t(:,:,jk)
-            END DO
-            CALL iom_put( "ttrd_zdfp", ztrdt )        ! PURE vertical diffusion (no isoneutral contribution)
-            CALL iom_put( "strd_zdfp", ztrds )
-         ENDIF
-         !
-         CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt, ztrds )
-         !
-         RETURN                     ! do nothing else if 3D output with IOM
-         !
-      ENDIF
-
+      !!----------------------------------------------------------------------
 
       ! I. Tracers trends
@@ -509,18 +460,18 @@
             WRITE (numout,*)
             WRITE (numout,9400) kt
-            WRITE (numout,9401)  tmo(jpicpt_xad) / tvolt, smo(jpicpt_xad) / tvolt
-            WRITE (numout,9411)  tmo(jpicpt_yad) / tvolt, 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)  tmo(jptra_trd_xad) / tvolt, smo(jptra_trd_xad) / tvolt
+            WRITE (numout,9411)  tmo(jptra_trd_yad) / tvolt, smo(jptra_trd_yad) / tvolt
+            WRITE (numout,9402)  tmo(jptra_trd_zad) / tvolt, smo(jptra_trd_zad) / tvolt
+            WRITE (numout,9403)  tmo(jptra_trd_ldf) / tvolt, smo(jptra_trd_ldf) / tvolt
+            WRITE (numout,9404)  tmo(jptra_trd_zdf) / tvolt, smo(jptra_trd_zdf) / tvolt
+            WRITE (numout,9405)  tmo(jptra_trd_npc) / tvolt, smo(jptra_trd_npc) / tvolt
+            WRITE (numout,9406)  tmo(jptra_trd_dmp) / tvolt, smo(jptra_trd_dmp) / tvolt
+            WRITE (numout,9407)  tmo(jptra_trd_qsr) / tvolt
+            WRITE (numout,9408)  tmo(jptra_trd_nsr) / tvolt, smo(jptra_trd_nsr) / 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
+            WRITE (numout,9410) (  tmo(jptra_trd_xad) + tmo(jptra_trd_yad) + tmo(jptra_trd_zad) + tmo(jptra_trd_ldf) + tmo(jptra_trd_zdf)   &
+            &                    + tmo(jptra_trd_npc) + tmo(jptra_trd_dmp) + tmo(jptra_trd_qsr) + tmo(jptra_trd_nsr) ) / tvolt,   &
+            &                   (  smo(jptra_trd_xad) + smo(jptra_trd_yad) + smo(jptra_trd_zad) + smo(jptra_trd_ldf) + smo(jptra_trd_zdf)   &
+            &                    + smo(jptra_trd_npc) + smo(jptra_trd_dmp)                   + smo(jptra_trd_nsr) ) / tvolt
          ENDIF
 
@@ -544,18 +495,18 @@
             WRITE (numout,*)
             WRITE (numout,9420) kt
-            WRITE (numout,9421)   t2(jpicpt_xad) / tvolt, s2(jpicpt_xad) / tvolt
-            WRITE (numout,9431)   t2(jpicpt_yad) / tvolt, 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)   t2(jptra_trd_xad) / tvolt, s2(jptra_trd_xad) / tvolt
+            WRITE (numout,9431)   t2(jptra_trd_yad) / tvolt, s2(jptra_trd_yad) / tvolt
+            WRITE (numout,9422)   t2(jptra_trd_zad) / tvolt, s2(jptra_trd_zad) / tvolt
+            WRITE (numout,9423)   t2(jptra_trd_ldf) / tvolt, s2(jptra_trd_ldf) / tvolt
+            WRITE (numout,9424)   t2(jptra_trd_zdf) / tvolt, s2(jptra_trd_zdf) / tvolt
+            WRITE (numout,9425)   t2(jptra_trd_npc) / tvolt, s2(jptra_trd_npc) / tvolt
+            WRITE (numout,9426)   t2(jptra_trd_dmp) / tvolt, s2(jptra_trd_dmp) / tvolt
+            WRITE (numout,9427)   t2(jptra_trd_qsr) / tvolt
+            WRITE (numout,9428)   t2(jptra_trd_nsr) / tvolt, s2(jptra_trd_nsr) / 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
+            WRITE (numout,9430) (  t2(jptra_trd_xad) + t2(jptra_trd_yad) + t2(jptra_trd_zad) + t2(jptra_trd_ldf) + t2(jptra_trd_zdf)   &
+            &                    + t2(jptra_trd_npc) + t2(jptra_trd_dmp) + t2(jptra_trd_qsr) + t2(jptra_trd_nsr) ) / tvolt,   &
+            &                   (  s2(jptra_trd_xad) + s2(jptra_trd_yad) + s2(jptra_trd_zad) + s2(jptra_trd_ldf) + s2(jptra_trd_zdf)   &
+            &                    + s2(jptra_trd_npc) + s2(jptra_trd_dmp)                  + s2(jptra_trd_nsr) ) / tvolt
          ENDIF
 
@@ -579,15 +530,15 @@
             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
+            WRITE (numout,9441) ( tmo(jptra_trd_xad)+tmo(jptra_trd_yad)+tmo(jptra_trd_zad) )/tvolt,   &
+            &                   ( smo(jptra_trd_xad)+smo(jptra_trd_yad)+smo(jptra_trd_zad) )/tvolt
+            WRITE (numout,9442)   tmo(jptra_trd_sad)/tvolt,  smo(jptra_trd_sad)/tvolt
+            WRITE (numout,9443)   tmo(jptra_trd_ldf)/tvolt,  smo(jptra_trd_ldf)/tvolt
+            WRITE (numout,9444)   tmo(jptra_trd_zdf)/tvolt,  smo(jptra_trd_zdf)/tvolt
+            WRITE (numout,9445)   tmo(jptra_trd_npc)/tvolt,  smo(jptra_trd_npc)/tvolt
+            WRITE (numout,9446) ( t2(jptra_trd_xad)+t2(jptra_trd_yad)+t2(jptra_trd_zad) )/tvolt,    &
+            &                   ( s2(jptra_trd_xad)+s2(jptra_trd_yad)+s2(jptra_trd_zad) )/tvolt
+            WRITE (numout,9447)   t2(jptra_trd_ldf)/tvolt,   s2(jptra_trd_ldf)/tvolt
+            WRITE (numout,9448)   t2(jptra_trd_zdf)/tvolt,   s2(jptra_trd_zdf)/tvolt
+            WRITE (numout,9449)   t2(jptra_trd_npc)/tvolt,   s2(jptra_trd_npc)/tvolt
          ENDIF
 
@@ -609,37 +560,4 @@
    END SUBROUTINE trd_twr
 
-#   else
-   !!----------------------------------------------------------------------
-   !!   Default case :                                         Empty module
-   !!----------------------------------------------------------------------
-   INTERFACE trd_icp
-      MODULE PROCEDURE trd_2d, trd_3d
-   END INTERFACE
-
-CONTAINS
-   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype )       ! Empty routine
-      REAL, DIMENSION(:,:) ::   ptrd2dx, ptrd2dy
-      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
-      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
-      WRITE(*,*) 'trd_2d: You should not have seen this print! error ?', &
-          &       ptrd2dx(1,1), ptrd2dy(1,1), ktrd, ctype
-   END SUBROUTINE trd_2d
-   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd , ctype )       ! Empty routine
-      REAL, DIMENSION(:,:,:) ::   ptrd3dx, ptrd3dy
-      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
-      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
-      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', &
-          &       ptrd3dx(1,1,1), ptrd3dy(1,1,1), ktrd, ctype
-   END SUBROUTINE trd_3d
-   SUBROUTINE trd_icp_init               ! Empty routine
-   END SUBROUTINE trd_icp_init
-   SUBROUTINE trd_dwr( kt )          ! Empty routine
-      WRITE(*,*) 'trd_dwr: You should not have seen this print! error ?', kt
-   END SUBROUTINE trd_dwr
-   SUBROUTINE trd_twr( kt )          ! Empty routine
-      WRITE(*,*) 'trd_twr: You should not have seen this print! error ?', kt
-   END SUBROUTINE trd_twr
-#endif
-
    !!======================================================================
 END MODULE trdicp
Index: anches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp_oce.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp_oce.F90	(revision 3316)
+++ 	(revision )
@@ -1,88 +1,0 @@
-MODULE trdicp_oce
-   !!======================================================================
-   !!                   ***  MODULE trdicp_oce  ***
-   !! Ocean trends :   set tracer and momentum trend variables
-   !!======================================================================
-   !! History :  1.0  !  2004-08 (C. Talandier) New trends organization
-   !!            3.5  !  2012-02 (G. Madec)  add 3D tracer zdf trend output using iom
-   !!----------------------------------------------------------------------
-   USE par_oce                 ! ocean parameters
-
-   IMPLICIT NONE
-   PUBLIC
-
-#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
-
-   !                                        !!! => tracer trends indexes <=
-   INTEGER, PARAMETER ::   jpicpt_xad =  1   !: x- horizontal advection
-   INTEGER, PARAMETER ::   jpicpt_yad =  2   !: y- horizontal advection
-   INTEGER, PARAMETER ::   jpicpt_zad =  3   !: z- vertical   advection
-   INTEGER, PARAMETER ::   jpicpt_ldf =  4   !: lateral       diffusion
-   INTEGER, PARAMETER ::   jpicpt_zdf =  5   !: vertical diffusion (Kz)
-   INTEGER, PARAMETER ::   jpicpt_bbc =  6   !: Bottom Boundary Condition (geoth. flux) 
-   INTEGER, PARAMETER ::   jpicpt_bbl =  7   !: Bottom Boundary Layer (diffusive/convective)
-   INTEGER, PARAMETER ::   jpicpt_npc =  8   !: static instability mixing
-   INTEGER, PARAMETER ::   jpicpt_dmp =  9   !: damping
-   INTEGER, PARAMETER ::   jpicpt_qsr = 10   !: penetrative solar radiation
-   INTEGER, PARAMETER ::   jpicpt_nsr = 11   !: non solar radiation
-   INTEGER, PARAMETER ::   jpicpt_zl1 = 12   !: first level vertical flux
-
-   !                                        !!! => Total tracer trends indexes <=
-   INTEGER, PARAMETER ::   jptot_tra  = 12   !: change it when adding/removing one indice above
-   
-   !                                        !!! => dynamic trends indexes <=
-   INTEGER, PARAMETER ::   jpicpd_hpg =  1   !: hydrostatic pressure gradient 
-   INTEGER, PARAMETER ::   jpicpd_keg =  2   !: kinetic energy gradient
-   INTEGER, PARAMETER ::   jpicpd_rvo =  3   !: relative vorticity
-   INTEGER, PARAMETER ::   jpicpd_pvo =  4   !: planetary vorticity
-   INTEGER, PARAMETER ::   jpicpd_ldf =  5   !: lateral diffusion
-   INTEGER, PARAMETER ::   jpicpd_had =  6   !: horizontal advection
-   INTEGER, PARAMETER ::   jpicpd_zad =  7   !: vertical advection
-   INTEGER, PARAMETER ::   jpicpd_zdf =  8   !: vertical diffusion
-   INTEGER, PARAMETER ::   jpicpd_spg =  9   !: surface pressure gradient
-   INTEGER, PARAMETER ::   jpicpd_dat = 10   !: damping term
-   INTEGER, PARAMETER ::   jpicpd_swf = 11   !: surface wind forcing
-   INTEGER, PARAMETER ::   jpicpd_bfr = 12   !: bottom friction 
-
-   !                                        !!! => Total dynamic trends indexes <=
-   INTEGER, PARAMETER ::   jptot_dyn  = 12   !: 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
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!======================================================================
-END MODULE trdicp_oce
Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod.F90	(revision 3317)
@@ -11,23 +11,23 @@
 #if  defined key_trdtra || defined key_trddyn || defined key_trdmld || defined key_trdvor || defined key_esopa
    !!----------------------------------------------------------------------
-   !!   trd_mod          : manage the type of trend diagnostics
-   !!   trd_3Diom        : output 3D momentum and/or tracer trends using IOM
-   !!   trd_budget       : domain averaged budget of trends (including kinetic energy and tracer variance trends)
-   !!   trd_mod_init     : Initialization step
-   !!----------------------------------------------------------------------
-   USE oce                     ! ocean dynamics and tracers variables
-   USE dom_oce                 ! ocean space and time domain variables
-   USE zdf_oce                 ! ocean vertical physics variables
-   USE trdmod_oce              ! ocean variables trends
-   USE ldftra_oce              ! ocean active tracers lateral physics
-   USE sbc_oce                 ! surface boundary condition: ocean
-   USE phycst                  ! physical constants
-   USE trdvor                  ! ocean vorticity trends 
-   USE trdicp                  ! ocean bassin integral constraints properties
-   USE trdmld                  ! ocean active mixed layer tracers trends 
-   USE in_out_manager          ! I/O manager
-   USE iom					   ! I/O manager library
-   USE lib_mpp                 ! MPP library
-   USE wrk_nemo                ! Memory allocation
+   !!   trd_mod       : manage the type of trend diagnostics
+   !!   trd_3Diom     : output 3D momentum and/or tracer trends using IOM
+   !!   trd_mod_init  : Initialization step
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain variables
+   USE zdf_oce        ! ocean vertical physics variables
+   USE trdmod_oce     ! ocean variables trends
+   USE zdfbfr         ! bottom friction
+   USE ldftra_oce     ! ocean active tracers lateral physics
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE phycst         ! physical constants
+   USE trdvor         ! ocean vorticity trends 
+   USE trdicp         ! ocean bassin integral constraints properties
+   USE trdmld         ! ocean active mixed layer tracers trends 
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! Memory allocation
 
    IMPLICIT NONE
@@ -89,18 +89,18 @@
       IF( lk_trdvor .AND. ctype == 'DYN' ) THEN    ! II. Vorticity trends
          !                                         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-         SELECT CASE ( ktrd ) 
-         CASE ( jpdyn_trd_hpg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_prg )   ! Hydrostatique Pressure Gradient 
-         CASE ( jpdyn_trd_keg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_keg )   ! KE Gradient 
-         CASE ( jpdyn_trd_rvo )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_rvo )   ! Relative Vorticity 
-         CASE ( jpdyn_trd_pvo )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_pvo )   ! Planetary Vorticity Term 
-         CASE ( jpdyn_trd_ldf )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_ldf )   ! Horizontal Diffusion 
-         CASE ( jpdyn_trd_zad )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_zad )   ! Vertical Advection 
-         CASE ( jpdyn_trd_spg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_spg )   ! Surface Pressure Grad. 
-         CASE ( jpdyn_trd_zdf )                                                      ! Vertical Diffusion 
+         SELECT CASE( ktrd ) 
+         CASE( jpdyn_trd_hpg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_prg )   ! Hydrostatique Pressure Gradient 
+         CASE( jpdyn_trd_keg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_keg )   ! KE Gradient 
+         CASE( jpdyn_trd_rvo )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_rvo )   ! Relative Vorticity 
+         CASE( jpdyn_trd_pvo )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_pvo )   ! Planetary Vorticity Term 
+         CASE( jpdyn_trd_ldf )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_ldf )   ! Horizontal Diffusion 
+         CASE( jpdyn_trd_zad )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_zad )   ! Vertical Advection 
+         CASE( jpdyn_trd_spg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_spg )   ! Surface Pressure Grad. 
+         CASE( jpdyn_trd_zdf )                                                      ! Vertical Diffusion 
             ztswu(:,:) = 0.e0   ;   ztswv(:,:) = 0.e0
             DO jj = 2, jpjm1                                                             ! wind stress trends
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ztswu(ji,jj) = ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 )
-                  ztswv(ji,jj) = ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 )
+                  ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 )
+                  ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 )
                END DO
             END DO
@@ -157,78 +157,4 @@
 
 
-   SUBROUTINE trd_budget( ptrdx, ptrdy, ktrd, ctype, kt )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE trd_budget  ***
-      !! 
-      !! ** Purpose : integral constraint diagnostics for momentum and/or tracer trends
-      !!              
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend 
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
-      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
-      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type 'DYN'/'TRA'
-      INTEGER                   , INTENT(in   ) ::   kt      ! time step
-      !!
-      INTEGER ::   ji, jj   ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace 
-      !!----------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
-
-      IF( MOD(kt,nn_trd) == 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)                          ! non solar radiation
-                                         z2dy(:,:) = ptrdy(:,:,1)
-                                         CALL trd_icp( z2dx , z2dy , jpicpt_nsr, ctype )
-            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 )   ! z- vertical adv 
-                                         !                                                 ! surface flux
-                                         IF( lk_vvl ) THEN                                      ! variable volume = zero
-                                            z2dx(:,:) = 0._wp
-                                            z2dy(:,:) = 0._wp
-                                         ELSE                                                   ! constant volume = wn*tsn/e3t
-                                            z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / fse3t(:,:,1)
-                                            z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / fse3t(:,:,1)
-                                         ENDIF
-                                         CALL trd_icp( z2dx , z2dy , jpicpt_zl1, ctype ) 
-            END SELECT
-         ENDIF
-
-         IF( lk_trddyn .AND. ctype == 'DYN' ) THEN       ! momentum trends 
-            !
-            SELECT CASE ( ktrd )
-            CASE( jpdyn_trd_hpg )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_hpg, ctype )   ! hydrost. pressure gradient
-            CASE( jpdyn_trd_spg )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_spg, ctype )   ! surface pressure grad.
-            CASE( jpdyn_trd_pvo )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_pvo, ctype )   ! planetary vorticity
-            CASE( jpdyn_trd_rvo )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_rvo, ctype )   ! relative  vorticity or metric term
-            CASE( jpdyn_trd_keg )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_keg, ctype )   ! KE gradient         or hor. advection
-            CASE( jpdyn_trd_zad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_zad, ctype )   ! vertical  advection 
-            CASE( jpdyn_trd_ldf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_ldf, ctype )   ! lateral  diffusion 
-            CASE( jpdyn_trd_zdf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_zdf, ctype )   ! vertical diffusion (icluding bfr & tau)
-                                        ztswu(:,:) = ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(:,:,1) * rau0 )
-                                        ztswv(:,:) = ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(:,:,1) * rau0 )
-                                        CALL trd_icp( ztswu, ztswv, jpicpd_swf, ctype )   ! wind stress trends
-            CASE( jpdyn_trd_bfr )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_bfr, ctype )   ! bottom friction trends
-            END SELECT
-            !
-         ENDIF
-         !
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
-      !
-   END SUBROUTINE trd_budget
-
-
    SUBROUTINE trd_3Diom( ptrdx, ptrdy, ktrd, ctype, kt )
       !!---------------------------------------------------------------------
@@ -244,6 +170,7 @@
       !!
       INTEGER ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy, ztswu, ztswv    ! 2D workspace 
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3dx, z3dy                  ! 3D workspace 
+      INTEGER ::   ikbu, ikbv   ! local integers
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy, ztswu, ztswv   ! 2D workspace 
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3dx, z3dy                 ! 3D workspace 
       !!----------------------------------------------------------------------
 
@@ -253,29 +180,35 @@
          !
          SELECT CASE( ktrd )
-         CASE( jptra_trd_xad )   ;   CALL iom_put( "ttrd_xad", ptrdx )        ! x- horizontal advection
-                                     CALL iom_put( "strd_xad", ptrdy )
-         CASE( jptra_trd_yad )   ;   CALL iom_put( "ttrd_yad", ptrdx )        ! y- horizontal advection
-                                     CALL iom_put( "strd_yad", ptrdy )
-         CASE( jptra_trd_zad )   ;   CALL iom_put( "ttrd_zad", ptrdx )        ! z- vertical   advection
-                                     CALL iom_put( "strd_zad", ptrdy )
-                                     IF( .NOT.lk_vvl ) THEN                   ! cst volume : adv flux through z=0 surface
-                                        z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / fse3t(:,:,1)
-                                        z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / fse3t(:,:,1)
-                                        CALL iom_put( "ttrd_sad", z2dx )
-                                        CALL iom_put( "strd_sad", z2dy )
-                                     ENDIF
-         CASE( jptra_trd_ldf )   ;   CALL iom_put( "ttrd_ldf", ptrdx )        ! lateral diffusion
-                                     CALL iom_put( "strd_ldf", ptrdy )
-         CASE( jptra_trd_zdf )   ;   CALL iom_put( "ttrd_zdf", ptrdx )        ! vertical diffusion (including Kz contribution)
-                                     CALL iom_put( "strd_zdf", ptrdy )
-         CASE( jptra_trd_dmp )   ;   CALL iom_put( "ttrd_dmp", ptrdx )        ! internal restoring (damping)
-                                     CALL iom_put( "strd_dmp", ptrdy )
-         CASE( jptra_trd_bbl )   ;   CALL iom_put( "ttrd_bbl", ptrdx )        ! bottom boundary layer
-                                     CALL iom_put( "strd_bbl", ptrdy )
-         CASE( jptra_trd_npc )   ;   CALL iom_put( "ttrd_npc", ptrdx )        ! static instability mixing
-                                     CALL iom_put( "strd_npc", ptrdy )
-         CASE( jptra_trd_qsr )   ;   CALL iom_put( "ttrd_qsr", ptrdx )        ! penetrative solar radiat. (only on temperature)
-         CASE( jptra_trd_nsr )   ;   CALL iom_put( "ttrd_qns", ptrdx(:,:,1) ) ! non-solar     radiation   (only on temperature)
-         CASE( jptra_trd_bbc )   ;   CALL iom_put( "ttrd_bbc", ptrdx )        ! geothermal heating   (only on temperature)
+         CASE( jptra_trd_xad  )   ;   CALL iom_put( "ttrd_xad" , ptrdx )        ! x- horizontal advection
+                                      CALL iom_put( "strd_xad" , ptrdy )
+         CASE( jptra_trd_yad  )   ;   CALL iom_put( "ttrd_yad" , ptrdx )        ! y- horizontal advection
+                                      CALL iom_put( "strd_yad" , ptrdy )
+         CASE( jptra_trd_zad  )   ;   CALL iom_put( "ttrd_zad" , ptrdx )        ! z- vertical   advection
+                                      CALL iom_put( "strd_zad" , ptrdy )
+                                      IF( .NOT.lk_vvl ) THEN                   ! cst volume : adv flux through z=0 surface
+                                         z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / fse3t(:,:,1)
+                                         z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / fse3t(:,:,1)
+                                         CALL iom_put( "ttrd_sad", z2dx )
+                                         CALL iom_put( "strd_sad", z2dy )
+                                      ENDIF
+         CASE( jptra_trd_ldf  )   ;   CALL iom_put( "ttrd_ldf" , ptrdx )        ! lateral diffusion
+                                      CALL iom_put( "strd_ldf" , ptrdy )
+         CASE( jptra_trd_zdf  )   ;   CALL iom_put( "ttrd_zdf" , ptrdx )        ! vertical diffusion (including Kz contribution)
+                                      CALL iom_put( "strd_zdf" , ptrdy )
+         CASE( jptra_trd_zdfp )   ;   CALL iom_put( "ttrd_zdfp", ptrdx )        ! PURE vertical diffusion (no isoneutral contribution)
+                                      CALL iom_put( "strd_zdfp", ptrdy )
+         CASE( jptra_trd_dmp  )   ;   CALL iom_put( "ttrd_dmp" , ptrdx )        ! internal restoring (damping)
+                                      CALL iom_put( "strd_dmp" , ptrdy )
+         CASE( jptra_trd_bbl  )   ;   CALL iom_put( "ttrd_bbl" , ptrdx )        ! bottom boundary layer
+                                      CALL iom_put( "strd_bbl" , ptrdy )
+         CASE( jptra_trd_npc  )   ;   CALL iom_put( "ttrd_npc" , ptrdx )        ! static instability mixing
+                                      CALL iom_put( "strd_npc" , ptrdy )
+         CASE( jptra_trd_nsr  )   ;   CALL iom_put( "ttrd_qns" , ptrdx )        ! surface forcing + runoff (ln_rnf=T)
+                                      CALL iom_put( "strd_cdt" , ptrdy )
+         CASE( jptra_trd_qsr  )   ;   CALL iom_put( "ttrd_qsr" , ptrdx )        ! penetrative solar radiat. (only on temperature)
+         CASE( jptra_trd_bbc  )   ;   CALL iom_put( "ttrd_bbc" , ptrdx )        ! geothermal heating   (only on temperature)
+         CASE( jptra_trd_atf  )   ;   CALL iom_put( "ttrd_atf" , ptrdx )        ! asselin time Filter
+                                      CALL iom_put( "strd_atf" , ptrdy )
+
          END SELECT
       ENDIF
@@ -322,6 +255,27 @@
                                      CALL iom_put( "utrd_tau", z2dx )
                                      CALL iom_put( "vtrd_tau", z2dy )
-         CASE( jpdyn_trd_bfr )   ;   CALL iom_put( "utrd_bfr", ptrdx )    ! bottom friction term
-                                     CALL iom_put( "vtrd_bfr", ptrdy )
+         CASE( jpdyn_trd_bfr )
+            IF( .NOT.ln_bfrimp )     CALL iom_put( "utrd_bfr", ptrdx )    ! bottom friction (explicit case)
+            IF( .NOT.ln_bfrimp )     CALL iom_put( "vtrd_bfr", ptrdy )
+!!gm only valid if ln_bfrimp=T otherwise the bottom stress as to be recomputed....
+
+         CASE( jpdyn_trd_atf )   ;   CALL iom_put( "utrd_atf", ptrdx )    ! asselin filter trends 
+                                     CALL iom_put( "vtrd_atf", ptrdy )
+            IF( ln_bfrimp ) THEN                                          ! bottom friction (implicit case)
+					z3dx(:,:,:) = 0._wp   ;   z3dy(:,:,:) = 0._wp                 ! after velocity known (now filed at this stage)
+               DO jk = 1, jpkm1
+                  DO jj = 2, jpjm1
+                     DO ji = 2, jpim1
+                        ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
+                        ikbv = mbkv(ji,jj)
+                        z3dx(ji,jj,jk) = bfrua(ji,jj) * un(ji,jj,ikbu) / fse3u(ji,jj,ikbu)
+                        z3dy(ji,jj,jk) = bfrva(ji,jj) * vn(ji,jj,ikbv) / fse3v(ji,jj,ikbv)
+                     END DO
+                  END DO
+               END DO
+                                     CALL iom_put( "utrd_bfr", z3dx )    ! bottom friction (implicit)
+                                     CALL iom_put( "vtrd_bfr", z3dy )
+            ENDIF
+            !
          END SELECT
          !
Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod_oce.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod_oce.F90	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod_oce.F90	(revision 3317)
@@ -6,5 +6,5 @@
    !! History :  1.0  !  2004-08  (C. Talandier) Original code
    !!----------------------------------------------------------------------
-   USE trdicp_oce              ! ocean momentum/tracers bassin properties trends variables
+   USE par_oce                 ! ocean parameters
    USE trdmld_oce              ! ocean active mixed layer tracers trends variables
    USE trdvor_oce              ! ocean vorticity trends variables
@@ -12,4 +12,18 @@
    IMPLICIT NONE
    PUBLIC
+
+#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
 
    !                                                     !!* Namelist namtrd:  diagnostics on dynamics/tracer trends *
@@ -46,35 +60,42 @@
    LOGICAL , PUBLIC ::   l_trdtrc = .FALSE.               !: tracers  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
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_bbc =  6     !: Bottom Boundary Condition (geoth. heating) 
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_bbl =  7     !: Bottom Boundary Layer (diffusive and/or advective)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_npc =  8     !: non-penetrative convection treatment
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_dmp =  9     !: internal restoring (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 time filter
+   !                                                      !!!* Active tracers trends indexes
+   INTEGER, PUBLIC, PARAMETER ::   jptot_tra      = 14     !: Total trend nb: change it when adding/removing one indice below
+   !                               ===================     !  
+   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_sad  =  4     !: z- vertical   advection
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_ldf  =  5     !: lateral       diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_zdf  =  6     !: vertical      diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_zdfp =  7     !: "PURE" vert.  diffusion (ln_traldf_iso=T)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_bbc  =  8     !: Bottom Boundary Condition (geoth. heating) 
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_bbl  =  9     !: Bottom Boundary Layer (diffusive and/or advective)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_npc  = 10     !: non-penetrative convection treatment
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_dmp  = 11     !: internal restoring (damping)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_qsr  = 12     !: penetrative solar radiation
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_nsr  = 13     !: non solar radiation / C/D on salinity  (+runoff if ln_rnf=T)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_atf  = 14     !: Asselin time filter
    !
-   !                                                     !!!* Passive tracers trends indexes (use if "key_top" defined)
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_sms  = 13    !: sources m. sinks
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_radn = 14    !: corr. trn<0 in trcrad
-   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_radb = 15    !: corr. trb<0 in trcrad (like atf)
+   !                                                      !!!* Passive tracers trends indices (use if "key_top" defined)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_sms  = 13     !: sources m. sinks
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_radn = 14     !: corr. trn<0 in trcrad
+   INTEGER, PUBLIC, PARAMETER ::   jptra_trd_radb = 15     !: corr. trb<0 in trcrad (like atf)
    !
-   !                                                     !!!* Momentum trends indexes
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_hpg =  1     !: hydrostatic pressure gradient 
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_spg =  2     !: surface     pressure gradient
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_keg =  3     !: kinetic energy gradient  or horizontal advection
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_rvo =  4     !: relative  vorticity      or metric term
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_pvo =  5     !: planetary vorticity
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_zad =  6     !: vertical advection
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_ldf =  7     !: horizontal diffusion   
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_zdf =  8     !: vertical   diffusion
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_swf =  9     !: surface stress
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_bfr = 10     !: bottom  stress 
-
+   !                                                      !!!* Momentum trends indices
+   INTEGER, PUBLIC, PARAMETER ::   jptot_dyn      = 11     !: Total trend nb: change it when adding/removing one indice below
+   !                               ===================     !  
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_hpg  =  1     !: hydrostatic pressure gradient 
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_spg  =  2     !: surface     pressure gradient
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_keg  =  3     !: kinetic energy gradient  or horizontal advection
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_rvo  =  4     !: relative  vorticity      or metric term
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_pvo  =  5     !: planetary vorticity
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_zad  =  6     !: vertical advection
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_ldf  =  7     !: horizontal diffusion   
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_zdf  =  8     !: vertical   diffusion
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_tau  =  9     !: surface stress
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_bfr  = 10     !: bottom  stress 
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_trd_atf  = 11     !: Asselin time filter
+   !
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 3317)
@@ -12,8 +12,12 @@
    !!   trd_tra_adv   : transform a div(U.T) trend into a U.grad(T) trend
    !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean domain 
+   USE zdf_oce        ! ocean vertical physics
+   USE zdfddm         ! vertical physics: double diffusion
    USE trdmod_oce     ! ocean active mixed layer tracers trends 
    USE trdmod         ! ocean active mixed layer tracers trends 
    USE trdmod_trc     ! ocean passive mixed layer tracers trends 
+   USE ldftra_oce     ! ocean active tracers lateral physics
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
@@ -29,4 +33,5 @@
    !! * Substitutions
 #  include "domzgr_substitute.h90"
+#  include "zdfddm_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -38,7 +43,7 @@
 
    INTEGER FUNCTION trd_tra_alloc()
-      !!----------------------------------------------------------------------------
+      !!---------------------------------------------------------------------
       !!                  ***  FUNCTION trd_tra_alloc  ***
-      !!----------------------------------------------------------------------------
+      !!---------------------------------------------------------------------
       ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , STAT= trd_tra_alloc )
       !
@@ -69,13 +74,14 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   ptra    ! now tracer variable
       !
-      REAL(wp), POINTER, DIMENSION(:,:,:)  ::  ztrds
+      INTEGER  ::   jk   ! loop indices
+      REAL(wp), POINTER, DIMENSION(:,:,:)  ::   zwt, zws, ztrdt, ztrds   ! 3D workspace
       !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi, jpj, jpk, ztrds )
-      !
+      !      
       IF( .NOT. ALLOCATED( trdtx ) ) THEN      ! allocate trdtra arrays
          IF( trd_tra_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trd_tra : unable to allocate arrays' )
       ENDIF
-      
+
       IF( ctype == 'TRA' .AND. ktra == jp_tem ) THEN   !==  Temperature trend  ==!
          !
@@ -112,4 +118,29 @@
          ENDIF
          !
+         IF( ktrd == jptra_trd_zdfp ) THEN     ! diagnose the "PURE" Kz trend (here: just before the swap)
+            !
+            IF( ln_traldf_iso ) THEN      ! iso-neutral diffusion only otherwise jptra_trd_zdf is "PURE"
+               !
+               CALL wrk_alloc( jpi, jpj, jpk, zwt, zws, ztrdt )
+               !
+               zwt(:,:, 1 ) = 0._wp   ;   zws(:,:, 1 ) = 0._wp            ! vertical diffusive fluxes
+               zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp
+               DO jk = 2, jpk
+                  zwt(:,:,jk) =   avt(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk)
+                  zws(:,:,jk) = fsavs(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / fse3w(:,:,jk) * tmask(:,:,jk)
+               END DO
+               !
+               ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
+               DO jk = 1, jpkm1
+                  ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / fse3t(:,:,jk)
+                  ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / fse3t(:,:,jk) 
+               END DO
+               CALL trd_mod( ztrdt, ztrds, jptra_trd_zdfp, ctype, kt )  
+               !
+               CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt )
+               !
+            ENDIF
+            !
+         ENDIF
       ENDIF
 
Index: /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 3316)
+++ /branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 3317)
@@ -25,11 +25,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rotb ,  rotn           !: relative vorticity           [s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hdivb,  hdivn          !: horizontal divergence        [s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn            !: 4D T-S fields        [Celcius,psu] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn   , tsa    !: 4D T-S fields        [Celcius,psu] 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rn2b ,  rn2            !: brunt-vaisala frequency**2   [s-2]
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  tsa             !: 4D T-S trends fields & work array 
-   !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0  [no units]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhop   !: potential volumic mass                           [kg/m3]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0  [no units]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rhop   !: potential volumic mass                           [kg/m3]
 
    !! free surface                                      !  before  ! now    ! after  !
