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 3325)
+++ branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml	(revision 3326)
@@ -117,31 +117,57 @@
     <field id="Mmx_v"          description="Mm current barotrope along j-axis harmonic real part "      unit="m/s"        />
     <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" />
-    < 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" />
+    <!-- T & S trends (ln_tra_trd=T) -->
+    <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. (NOT key_vvl)"                unit="degC/s"                   />
+    <field id="strd_sad"     description="salinity   -trend: surface adv. (NOT key_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  (ln_traldf_iso=T)"   unit="degC/s" axis_ref="deptht" />
+    <field id="strd_zdfp"    description="salinity   -trend: pure vert. diffusion"                      unit="psu/s"  axis_ref="deptht" />
+    <      * internal restoring only >
+    <field id="ttrd_dmp"     description="temperature-trend: internal restoring       (ln_tradmp=T) "   unit="degC/s" axis_ref="deptht" />
+    <field id="strd_dmp"     description="salinity   -trend: internal restoring  (ln_tradmp=T) "        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 convection"                unit="degC/s" axis_ref="deptht" />
+    <field id="strd_npc"     description="salinity   -trend: non-penetrative convection"                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" />
     <>
-    <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" />
+    <!-- Kinetic energy trends (ln_PE_trd=T) -->
+    <field id="petrd_xad"     description="PE-trend: i-advection"                                       unit="XX/s" axis_ref="deptht" />
+    <field id="petrd_yad"     description="PE-trend: j-advection"                                       unit="XX/s" axis_ref="deptht" />
+    <field id="petrd_zad"     description="PE-trend: k-advection"                                       unit="XX/s" axis_ref="deptht" />
+    <field id="petrd_sad"     description="PE-trend: surface advection                 (NOT key_vvl)"   unit="XX/s"                   />
+    <field id="petrd_ldf"     description="PE-trend: lateral  diffusion"                                unit="XX/s" axis_ref="deptht" />
+    <field id="petrd_zdf"     description="PE-trend: vertical diffusion"                                unit="XX/s" axis_ref="deptht" />
+    <field id="petrd_zdfp"    description="PE-trend: pure vert. diffusion          (ln_traldf_iso=T)"   unit="XX/s" axis_ref="deptht" />
+    <field id="petrd_dmp"     description="PE-trend: internal restoring                (ln_tradmp=T)"   unit="XX/s" axis_ref="deptht" />
+    <field id="petrd_bbl"     description="PE-trend: bottom boundary layer"                             unit="XX/s" axis_ref="deptht" />
+    <field id="petrd_npc"     description="PE-trend: non-penetrative conv."                             unit="XX/s" axis_ref="deptht" />
+    <field id="petrd_qns"     description="PE-trend: non-solar flux + runoff"                           unit="XX/s" axis_ref="deptht" />
+    <field id="petrd_qsr"     description="PE-trend: solar penetr. heating"                             unit="XX/s" axis_ref="deptht" />
+    <field id="petrd_bbc"     description="PE-trend: geothermal heating"                                unit="XX/s" axis_ref="deptht" />
+    <>  
+    <!-- variables available with key_trddyn -->
+    <field id="ketrd_hpg"     description="PE-trend: hydrostatic pressure gradient"      unit="YY"                />
+    <field id="ketrd_spg"     description="PE-trend: surface     pressure gradient"      unit="YY"                />
+    <field id="ketrd_rvo"     description="PE-trend: relative  vorticity or metric term" unit="YY"                />
+    <field id="ketrd_pvo"     description="PE-trend: planetary vorticity"                unit="YY"                />
+    <field id="ketrd_keg"     description="PE-trend: KE gradient         or hor. adv."   unit="YY"                />
+    <field id="ketrd_zad"     description="PE-trend: vertical  advection"                unit="YY"                />
+    <field id="ketrd_ldf"     description="PE-trend: lateral   diffusion"                unit="YY"                />
+    <field id="ketrd_zdf"     description="PE-trend: vertical  diffusion"                unit="YY"                />
+    <field id="ketrd_tau"     description="PE-trend: wind stress "                       unit="YY" axis_ref="none"/>
+    <field id="ketrd_bfr"     description="PE-trend: bottom friction"                    unit="YY"                />   
       </group>
 
Index: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 3325)
+++ branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 3326)
@@ -18,4 +18,5 @@
    !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
    !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_alpbet used in ldfslp
+   !!            3.5  ! 2012-03  (G. Madec)  add eos_drau_dtds used in trdpen
    !!----------------------------------------------------------------------
 
@@ -28,4 +29,5 @@
    !!   eos_bn2        : Compute the Brunt-Vaisala frequency
    !!   eos_alpbet     : calculates the in situ thermal/haline expansion ratio
+   !!   eos_drau_dtds  : calculates the partial derivative of in situ density with respect to T and S
    !!   tfreez         : Compute the surface freezing temperature
    !!   eos_init       : set eos parameters (namelist)
@@ -51,9 +53,10 @@
    END INTERFACE
 
-   PUBLIC   eos        ! called by step, istate, tranpc and zpsgrd modules
-   PUBLIC   eos_init   ! called by istate module
-   PUBLIC   bn2        ! called by step module
-   PUBLIC   eos_alpbet ! called by ldfslp module
-   PUBLIC   tfreez     ! called by sbcice_... modules
+   PUBLIC   eos            ! called by step, istate, tranpc and zpsgrd modules
+   PUBLIC   eos_init       ! called by istate module
+   PUBLIC   bn2            ! called by step module
+   PUBLIC   eos_alpbet     ! called by ldfslp module
+   PUBLIC   tfreez         ! called by sbcice_... modules
+   PUBLIC   eos_drau_dtds  ! called by trdpen module
 
    !                                          !!* Namelist (nameos) *
@@ -678,4 +681,62 @@
 
 
+   SUBROUTINE eos_drau_dtds( pts, pdr_dt, pdr_ds )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE eos_alpbet  ***
+      !!
+      !! ** Purpose :   Calculates the in situ thermal/haline expansion ratio at T-points
+      !!
+      !! ** Method  :   calculates alpha / beta ratio at T-points
+      !!       * nn_eos = 0  : UNESCO sea water properties
+      !!                       The alpha/beta ratio is returned as 3-D array palpbet using the polynomial
+      !!                       polynomial expression of McDougall (1987).
+      !!                       Scalar beta0 is returned = 1.
+      !!       * nn_eos = 1  : linear equation of state (temperature only)
+      !!                       The ratio is undefined, so we return alpha as palpbet
+      !!                       Scalar beta0 is returned = 0.
+      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
+      !!                       The alpha/beta ratio is returned as ralpbet
+      !!                       Scalar beta0 is returned = 1.
+      !!
+      !! ** Action  : - palpbet : thermal/haline expansion ratio at T-points
+      !!            :   beta0   : 1. or 0.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts      ! pot. temperature & salinity
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pdr_dt   ! partial derivative of in situ density
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pdr_ds   ! with respect to T and S, resp.
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zt, zs, zh   ! local scalars
+      !!----------------------------------------------------------------------
+      !
+      IF ( nn_timing == 1 )   CALL timing_start('eos_drau_dtds')
+      !
+      SELECT CASE ( nn_eos )
+      !
+      CASE ( 0 )               ! Jackett and McDougall (1994) formulation
+         !
+         CALL ctl_stop('eos_drau_dtds: full eos not coded yet')
+         !
+         !
+      CASE ( 1 )              !==  Linear formulation = F( temperature )  ==!
+         rn_beta(:,:,:) = rn_alpha
+         rn_beta(:,:,:) = 0._wp
+         !
+      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==!
+         pdr_dt(:,:,:) = rn_alpha
+         pdr_ds(:,:,:) = rn_beta
+         !
+      CASE DEFAULT
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
+         nstop = nstop + 1
+         !
+      END SELECT
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('eos_drau_dtds')
+      !
+   END SUBROUTINE eos_drau_dtds
+
+
    FUNCTION tfreez( psal ) RESULT( ptf )
       !!----------------------------------------------------------------------
Index: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 3325)
+++ branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 3326)
@@ -13,20 +13,20 @@
 
    !!----------------------------------------------------------------------
-   !!   tra_qsr      : trend due to the solar radiation penetration
-   !!   tra_qsr_init : solar radiation penetration initialization
+   !!   tra_qsr       : trend due to the solar radiation penetration
+   !!   tra_qsr_init  : solar radiation penetration initialization
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers
-   USE dom_oce         ! ocean space and time domain
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE trc_oce         ! share SMS/Ocean variables
-   USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
-   USE phycst          ! physical constants
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE iom             ! I/O manager
-   USE fldread         ! read input fields
-   USE restart         ! ocean restart
-   USE lib_mpp         ! MPP library
+   USE oce            ! ocean dynamics and active tracers
+   USE dom_oce        ! ocean space and time domain
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE trc_oce        ! share SMS/Ocean variables
+   USE trd_oce        ! trends: ocean variables
+   USE trdtra         ! trends manager: tracers 
+   USE phycst         ! physical constants
+   USE in_out_manager ! I/O manager
+   USE prtctl         ! Print control
+   USE iom            ! I/O manager
+   USE fldread        ! read input fields
+   USE restart        ! ocean restart
+   USE lib_mpp        ! MPP library
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
@@ -48,10 +48,9 @@
    REAL(wp), PUBLIC ::   rn_si1     = 23.0_wp   !: deepest depth of extinction (water type I)       (2 bands)
    
-   ! Module variables
-   REAL(wp) ::   xsi0r                           !: inverse of rn_si0
-   REAL(wp) ::   xsi1r                           !: inverse of rn_si1
+   INTEGER , PUBLIC ::   nksr   !: levels below which the light cannot penetrate ( depth larger than 391 m)
+
+   REAL(wp)                  ::   xsi0r, xsi1r        ! inverse of rn_si0 and rn_si1, resp.
+   REAL(wp), DIMENSION(3,61) ::   rkrgb               ! tabulated attenuation coefficients for RGB absorption
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
-   INTEGER, PUBLIC ::   nksr              ! levels below which the light cannot penetrate ( depth larger than 391 m)
-   REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption
 
    !! * Substitutions
@@ -92,5 +91,4 @@
       !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
       !!----------------------------------------------------------------------
-      !
       INTEGER, INTENT(in) ::   kt     ! ocean time-step
       !
@@ -116,5 +114,5 @@
       ENDIF
 
-      IF( l_trdtra ) THEN      ! Save ta and sa trends
+      IF( l_trdtra ) THEN      ! Save temperature trend
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
@@ -141,5 +139,4 @@
       !                                        Compute now qsr tracer content field
       !                                        ************************************
-      
       !                                           ! ============================================== !
       IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
@@ -168,5 +165,5 @@
             IF( nn_chldta == 1 .OR. lk_vvl ) THEN            !*  Variable Chlorophyll or ocean volume
                !
-               IF( nn_chldta == 1 ) THEN                             !*  Variable Chlorophyll
+               IF( nn_chldta == 1 ) THEN                             !- Variable Chlorophyll
                   !
                   CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step
@@ -184,18 +181,18 @@
                      END DO
                   END DO
-               ELSE                                            ! Variable ocean volume but constant chrlorophyll
-                  zchl = 0.05                                     ! constant chlorophyll
+               ELSE                                                  !- Variable ocean volume but constant chrlorophyll
+                  zchl = 0.05                                           ! constant chlorophyll
                   irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 )
-                  zekb(:,:) = rkrgb(1,irgb)                       ! Separation in R-G-B depending of the chlorophyll 
+                  zekb(:,:) = rkrgb(1,irgb)                             ! Separation in R-G-B depending of the chlorophyll 
                   zekg(:,:) = rkrgb(2,irgb)
                   zekr(:,:) = rkrgb(3,irgb)
                ENDIF
                !
-               zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B
-               ze0(:,:,1) = rn_abs  * qsr(:,:)
-               ze1(:,:,1) = zcoef * qsr(:,:)
-               ze2(:,:,1) = zcoef * qsr(:,:)
-               ze3(:,:,1) = zcoef * qsr(:,:)
-               zea(:,:,1) =         qsr(:,:)
+               zcoef  = ( 1. - rn_abs ) / 3.e0                       !- equi-partition in R-G-B
+               ze0(:,:,1) = rn_abs * qsr(:,:)
+               ze1(:,:,1) =  zcoef * qsr(:,:)
+               ze2(:,:,1) =  zcoef * qsr(:,:)
+               ze3(:,:,1) =  zcoef * qsr(:,:)
+               zea(:,:,1) =          qsr(:,:)
                !
                DO jk = 2, nksr+1
Index: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90
===================================================================
--- branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90	(revision 3325)
+++ branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90	(revision 3326)
@@ -88,4 +88,7 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zke                 ! 3D workspace 
       !!----------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpi, jpj     , z2dx, z2dy, zke2d )
+      CALL wrk_alloc( jpi, jpj, jpk, zke )
       !
       CALL lbc_lnk( putrd, 'U', -1. )   ;   CALL lbc_lnk( pvtrd, 'V', -1. )      ! lateral boundary conditions
Index: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdpen.F90
===================================================================
--- branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdpen.F90	(revision 3326)
+++ branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdpen.F90	(revision 3326)
@@ -0,0 +1,148 @@
+MODULE trdpen
+   !!======================================================================
+   !!                       ***  MODULE  trdpen  ***
+   !! Ocean diagnostics:  Potential ENnergy trends
+   !!=====================================================================
+   !! History :  3.5  !  2012-02  (G. Madec) original code 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   trd_pen       : compute and output Potential Energy trends from T & S trends
+   !!   trd_pen_init  : initialisation of PE trends
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean domain 
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE zdf_oce        ! ocean vertical physics
+   USE trd_oce        ! trends: ocean variables
+   USE eosbn2         ! equation of state and related derivatives
+   USE ldftra_oce     ! ocean active tracers lateral physics
+   USE zdfddm         ! vertical physics: double diffusion
+   USE phycst         ! physical constants
+   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
+   PRIVATE
+
+   PUBLIC   trd_pen        ! called by all trdtra module
+   PUBLIC   trd_pen_init   ! called by all nemogcm module
+
+   INTEGER ::   nkstp   ! current time step 
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   drau_dt, drau_ds   ! partial derivative of rau with respect to T and S
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "zdfddm_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: trdtra.F90 3318 2012-02-25 15:50:01Z gm $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION trd_pen_alloc()
+      !!---------------------------------------------------------------------
+      !!                  ***  FUNCTION trd_tra_alloc  ***
+      !!---------------------------------------------------------------------
+      ALLOCATE( drau_dt(jpi,jpj,jpk) , drau_ds(jpi,jpj,jpk) , STAT= trd_pen_alloc )
+      !
+      IF( lk_mpp             )   CALL mpp_sum ( trd_pen_alloc )
+      IF( trd_pen_alloc /= 0 )   CALL ctl_warn('trd_pen_alloc: failed to allocate arrays')
+   END FUNCTION trd_pen_alloc
+
+
+   SUBROUTINE trd_pen( ptrdx, ptrdy, ktrd, kt, pdt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_tra_mng  ***
+      !! 
+      !! ** Purpose :   Dispatch all trends computation, e.g. 3D output, integral
+      !!                constraints, barotropic vorticity, kinetic enrgy, 
+      !!                potential energy, and/or mixed layer budget.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   ptrdx, ptrdy   ! Temperature & Salinity trends
+      INTEGER                   , INTENT(in) ::   ktrd           ! tracer trend index
+      INTEGER                   , INTENT(in) ::   kt             ! time step index
+      REAL(wp)                  , INTENT(in) ::   pdt            ! time step [s]
+      !
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2d   ! 2D workspace 
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpe   ! 3D workspace 
+      !!----------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zpe )
+      !
+      IF ( nn_eos == 0 .AND.  kt /= nkstp ) THEN   ! full eos: set partial derivatives at the 1st call of kt time step
+         nkstp = kt
+         CALL eos_drau_dtds( tsn, drau_dt, drau_ds )
+      ENDIF
+      !
+      DO jk = 1, jpkm1
+         zpe(:,:,jk) = grav * fsde3w(:,:,jk) * (  drau_dt(:,:,jk) * ptrdx(:,:,jk)   &
+            &                                   + drau_ds(:,:,jk) * ptrdx(:,:,jk)  )
+      END DO
+
+      SELECT CASE ( ktrd )
+      CASE ( jptra_xad  )   ;   CALL iom_put( "petrd_xad", zke )   ! zonal    advection
+      CASE ( jptra_yad  )   ;   CALL iom_put( "petrd_yad", zke )   ! merid.   advection
+      CASE ( jptra_zad  )   ;   CALL iom_put( "petrd_zad", zke )   ! vertical advection
+                                IF( .NOT.lk_vvl ) THEN                   ! cst volume : adv flux through z=0 surface
+                                   CALL wrk_alloc( jpi, jpj, z2d )
+                                   z2d(:,:) = wn(:,:,1) * (   drau_dt(:,:,1) * tsn(:,:,1,jp_tem) 
+                                      &                     + drau_ds(:,:,1) * tsn(:,:,1,jp_sal)  ) / fse3t(:,:,1)
+                                   CALL iom_put( "petrd_sad" , z2d )
+                                   CALL wrk_dealloc( jpi, jpj, z2d )
+                                ENDIF
+      CASE ( jptra_ldf  )   ;   CALL iom_put( "petrd_ldf" , zpe )   ! lateral  diffusion
+      CASE ( jptra_zdf  )   ;   CALL iom_put( "petrd_ldf" , zpe )   ! lateral  diffusion (K_z)
+      CASE ( jptra_zdfp )   ;   CALL iom_put( "petrd_zdfp", zpe )   ! vertical diffusion (K_z)
+      CASE ( jptra_dmp  )   ;   CALL iom_put( "petrd_dmp" , zpe )   ! internal 3D restoring (tradmp)
+      CASE ( jptra_bbl  )   ;   CALL iom_put( "petrd_bbl" , zpe )   ! bottom boundary layer
+      CASE ( jptra_npc  )   ;   CALL iom_put( "petrd_npc" , zpe )   ! non penetr convect adjustment
+      CASE ( jptra_nsr  )   ;   CALL iom_put( "petrd_for" , zpe )   ! surface forcing + runoff (ln_rnf=T)
+      CASE ( jptra_qsr  )   ;   CALL iom_put( "petrd_qsr" , zke )   ! air-sea : penetrative sol radiat
+      CASE ( jptra_bbc  )   ;   CALL iom_put( "petrd_bbc", zke )   ! bottom bound cond (geoth flux)
+      CASE ( jptra_atf  )   ;   CALL iom_put( "petrd_atf", zke )   ! asselin time filter (last trend)
+                                IF( .NOT.lk_vvl ) THEN                   ! cst volume : ssh term (otherwise include in e3t variation)
+                                   CALL wrk_alloc( jpi, jpj, z2d )
+                                   z2d(:,:) = ( ssha(:,:) - sshb(:,:) )   &
+                                      &     * (   drau_dt(:,:,1) * tsn(:,:,1,jp_tem) 
+                                      &         + drau_ds(:,:,1) * tsn(:,:,1,jp_sal)  ) / ( fse3t(:,:,1) * pdt )
+                                   CALL iom_put( "petrd_sad" , z2d )
+                                   CALL wrk_dealloc( jpi, jpj, z2d )
+                                ENDIF
+      END SELECT
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zpe )
+      !
+   END SUBROUTINE trd_pen
+
+
+   SUBROUTINE trd_pen_init
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trd_pen_init  ***
+      !! 
+      !! ** Purpose :   initialisation of 3D Kinetic Energy trend diagnostic
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'trd_pen_init : 3D Potential ENergy trends'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+      ENDIF
+      !                           ! allocate box volume arrays
+      IF ( trd_pen_alloc() /= 0 )   CALL ctl_stop('trd_pen_alloc: failed to allocate arrays')
+      !
+      IF ( lk_vvl               )   CALL ctl_stop('trd_pen_init : PE trends not coded for variable volume')
+      !
+      r1_2_rau0 = 0.5_wp / rau0
+      !
+   END SUBROUTINE trd_pen_init
+
+   !!======================================================================
+END MODULE trdtra
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 3325)
+++ branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 3326)
@@ -122,33 +122,33 @@
             END SELECT
          ELSE                            ! other trends: mask and send T & S trends to trd_tra_mng
-            ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
-            CALL trd_tra_mng( trdt, ztrds, ktrd, kt )  
-         ENDIF
-         !
-         IF( ktrd == jptra_zdfp ) THEN     ! diagnose the "PURE" Kz trend (here: just before the swap)
             !
-            IF( ln_traldf_iso ) THEN      ! iso-neutral diffusion only otherwise jptra_zdf is "PURE"
+            IF( ktrd == jptra_zdfp ) THEN     ! diagnose the "PURE" Kz trend (here: just before the swap)
                !
-               CALL wrk_alloc( jpi, jpj, jpk, zwt, zws, ztrdt )
+               IF( ln_traldf_iso ) THEN       ! iso-neutral diffusion only otherwise jptra_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_tra_mng( ztrdt, ztrds, jptra_zdfp, kt )  
+                  !
+                  CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt )
+                  !
+               ENDIF
                !
-               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_tra_mng( ztrdt, ztrds, jptra_zdfp, kt )  
-               !
-               CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt )
-               !
+            ELSE                              ! 
+               ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
+               CALL trd_tra_mng( trdt, ztrds, ktrd, kt )  
             ENDIF
-            !
-         ENDIF
       ENDIF
 
@@ -305,7 +305,7 @@
       SELECT CASE( ktrd )
       CASE( jptra_xad  )   ;   CALL iom_put( "ttrd_xad" , ptrdx )        ! x- horizontal advection
-                                   CALL iom_put( "strd_xad" , ptrdy )
+                               CALL iom_put( "strd_xad" , ptrdy )
       CASE( jptra_yad  )   ;   CALL iom_put( "ttrd_yad" , ptrdx )        ! y- horizontal advection
-                                   CALL iom_put( "strd_yad" , ptrdy )
+                               CALL iom_put( "strd_yad" , ptrdy )
       CASE( jptra_zad  )   ;   CALL iom_put( "ttrd_zad" , ptrdx )        ! z- vertical   advection
                                CALL iom_put( "strd_zad" , ptrdy )
