Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/NST_SRC/agrif_user.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/NST_SRC/agrif_user.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/NST_SRC/agrif_user.F90	(revision 2082)
@@ -307,6 +307,6 @@
 #if defined key_top
         ! Check passive tracer cell
-        IF( ndttrc .ne. 1 ) THEN
-           WRITE(*,*) 'ndttrc should be equal to 1'
+        IF( nn_dttrc .ne. 1 ) THEN
+           WRITE(*,*) 'nn_dttrc should be equal to 1'
         ENDIF
 #endif
@@ -444,6 +444,6 @@
          ENDIF
         ! Check passive tracer cell
-        IF( ndttrc .ne. 1 ) THEN
-           WRITE(*,*) 'ndttrc should be equal to 1'
+        IF( nn_dttrc .ne. 1 ) THEN
+           WRITE(*,*) 'nn_dttrc should be equal to 1'
         ENDIF
 
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/dtadyn.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/dtadyn.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/dtadyn.F90	(revision 2082)
@@ -184,16 +184,14 @@
          
          IF( lk_ldfslp ) THEN
-            ! Computes slopes
-            ! Caution : here tn, sn and avt are used as workspace
-            tn (:,:,:) = tdta  (:,:,:,2)
-            sn (:,:,:) = sdta  (:,:,:,2)
-            avt(:,:,:) = avtdta(:,:,:,2)
-         
-            CALL eos( tn, sn, rhd, rhop )   ! Time-filtered in situ density 
-            CALL bn2( tn, sn, rn2 )         ! before Brunt-Vaisala frequency
+            ! Computes slopes. Caution : here tsn and avt are used as workspace
+            tsn (:,:,:,jp_tem) = tdta  (:,:,:,2)
+            tsn (:,:,:,jp_sal) = sdta  (:,:,:,2)
+            avt(:,:,:)         = avtdta(:,:,:,2)
+         
+            CALL eos( tsn, rhd, rhop )   ! Time-filtered in situ density 
+            CALL bn2( tsn, rn2 )         ! before Brunt-Vaisala frequency
             IF( ln_zps )   &
-               &   CALL zps_hde( kt, tn , sn , rhd,  &  ! Partial steps: before Horizontal DErivative
-               &                 gtu, gsu, gru,  &  ! of t, s, rd at the bottom ocean level
-               &                 gtv, gsv, grv )
+               &   CALL zps_hde( kt, jpts, tsn, gtsu, gtsv,  &  ! Partial steps: before Horizontal DErivative
+               &                           rhd, gru , grv   )    ! of t, s, rd at the bottom ocean level
                    CALL zdf_mxl( kt )              ! mixed layer depth
                    CALL ldf_slp( kt, rhd, rn2 )
@@ -213,16 +211,14 @@
          
          IF( lk_ldfslp ) THEN
-            ! Computes slopes
-            ! Caution : here tn, sn and avt are used as workspace
-            tn (:,:,:) = tdta  (:,:,:,2)
-            sn (:,:,:) = sdta  (:,:,:,2)
-            avt(:,:,:) = avtdta(:,:,:,2)
-
-            CALL eos( tn, sn, rhd, rhop )   ! Time-filtered in situ density 
-            CALL bn2( tn, sn, rn2 )         ! before Brunt-Vaisala frequency
+            ! Computes slopes. Caution : here tsn and avt are used as workspace
+            tsn (:,:,:,jp_tem) = tdta  (:,:,:,2)
+            tsn (:,:,:,jp_sal) = sdta  (:,:,:,2)
+            avt(:,:,:)         = avtdta(:,:,:,2)
+         
+            CALL eos( tsn, rhd, rhop )   ! Time-filtered in situ density 
+            CALL bn2( tsn, rn2 )         ! before Brunt-Vaisala frequency
             IF( ln_zps )   &
-               &   CALL zps_hde( kt, tn , sn , rhd,  &  ! Partial steps: before Horizontal DErivative
-               &                 gtu, gsu, gru,  &  ! of t, s, rd at the bottom ocean level
-               &                 gtv, gsv, grv )
+               &   CALL zps_hde( kt, jpts, tsn, gtsu, gtsv,  &  ! Partial steps: before Horizontal DErivative
+               &                           rhd, gru , grv   )    ! of t, s, rd at the bottom ocean level
                    CALL zdf_mxl( kt )              ! mixed layer depth
                    CALL ldf_slp( kt, rhd, rn2 )
@@ -261,16 +257,14 @@
 
          IF( lk_ldfslp ) THEN
-            ! Computes slopes
-            ! Caution : here tn, sn and avt are used as workspace
-            tn (:,:,:) = tdta  (:,:,:,2)
-            sn (:,:,:) = sdta  (:,:,:,2)
-            avt(:,:,:) = avtdta(:,:,:,2)
-
-            CALL eos( tn, sn, rhd, rhop )   ! Time-filtered in situ density 
-            CALL bn2( tn, sn, rn2 )         ! before Brunt-Vaisala frequency
+            ! Computes slopes. Caution : here tsn and avt are used as workspace
+            tsn (:,:,:,jp_tem) = tdta  (:,:,:,2)
+            tsn (:,:,:,jp_sal) = sdta  (:,:,:,2)
+            avt(:,:,:)         = avtdta(:,:,:,2)
+         
+            CALL eos( tsn, rhd, rhop )   ! Time-filtered in situ density 
+            CALL bn2( tsn, rn2 )         ! before Brunt-Vaisala frequency
             IF( ln_zps )   &
-               &   CALL zps_hde( kt, tn , sn , rhd,  &  ! Partial steps: before Horizontal DErivative
-               &                 gtu, gsu, gru,  &  ! of t, s, rd at the bottom ocean level
-               &                 gtv, gsv, grv )
+               &   CALL zps_hde( kt, jpts, tsn, gtsu, gtsv,  &  ! Partial steps: before Horizontal DErivative
+               &                           rhd, gru , grv   )    ! of t, s, rd at the bottom ocean level
                    CALL zdf_mxl( kt )              ! mixed layer depth
                    CALL ldf_slp( kt, rhd, rn2 )
@@ -311,5 +305,5 @@
       
       ! In any case, we need rhop
-      CALL eos( tn, sn, rhd, rhop ) 
+      CALL eos( tsn, rhd, rhop ) 
       
 #if ! defined key_degrad && defined key_traldf_c2d
@@ -320,6 +314,5 @@
       ! Compute bbl coefficients if needed
       IF( lk_trabbl ) THEN
-         tb(:,:,:) = tn(:,:,:)
-         sb(:,:,:) = sn(:,:,:)
+         tsb(:,:,:,:) = tsn(:,:,:,:)
          CALL bbl( kt, 'TRC')
       END IF
@@ -735,7 +728,7 @@
       !!----------------------------------------------------------------------
       
-      tn (:,:,:) = tdta  (:,:,:,2)
-      sn (:,:,:) = sdta  (:,:,:,2)
-      avt(:,:,:) = avtdta(:,:,:,2)
+      tsn(:,:,:,jp_tem) = tdta  (:,:,:,2)
+      tsn(:,:,:,jp_sal) = sdta  (:,:,:,2)
+      avt(:,:,:)        = avtdta(:,:,:,2)
       
       un (:,:,:) = udta  (:,:,:,2) 
@@ -796,7 +789,7 @@
       zweighm1 = 1. - pweigh
       
-      tn (:,:,:) = zweighm1 * tdta  (:,:,:,1) + pweigh * tdta  (:,:,:,2)
-      sn (:,:,:) = zweighm1 * sdta  (:,:,:,1) + pweigh * sdta  (:,:,:,2)
-      avt(:,:,:) = zweighm1 * avtdta(:,:,:,1) + pweigh * avtdta(:,:,:,2)
+      tsn(:,:,:,jp_tem) = zweighm1 * tdta  (:,:,:,1) + pweigh * tdta  (:,:,:,2)
+      tsn(:,:,:,jp_sal) = zweighm1 * sdta  (:,:,:,1) + pweigh * sdta  (:,:,:,2)
+      avt(:,:,:)        = zweighm1 * avtdta(:,:,:,1) + pweigh * avtdta(:,:,:,2)
       
       un (:,:,:) = zweighm1 * udta  (:,:,:,1) + pweigh * udta  (:,:,:,2) 
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/istate.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/istate.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/istate.F90	(revision 2082)
@@ -59,6 +59,5 @@
       ;   hdivn(:,:,:) = 0.e0   ;
 
-      ;   tn   (:,:,:) = 0.e0   ;   ta   (:,:,:) = 0.e0
-      ;   sn   (:,:,:) = 0.e0   ;   sa   (:,:,:) = 0.e0
+      ;   tsn  (:,:,:,:) = 0.e0 
 
       rhd  (:,:,:) = 0.e0
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/opa.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/opa.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OFF_SRC/opa.F90	(revision 2082)
@@ -223,4 +223,5 @@
 
 
+      !                                     ! General initialization
                             CALL     phy_cst    ! Physical constants
                             CALL     eos_init   ! Equation of state
@@ -234,14 +235,16 @@
       IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
          &                  CALL zdf_ddm_init   ! double diffusive mixing
-      !                                     ! Lateral physics
 #if ! defined key_degrad
                             CALL ldf_tra_init   ! Lateral ocean tracer physics
 #endif
       IF( lk_ldfslp )       CALL ldf_slp_init   ! slope of lateral mixing
+
       !                                     ! Active tracers
                             CALL tra_qsr_init   ! penetrative solar radiation qsr
       IF( lk_trabbl     )   CALL tra_bbl_init   ! advective (and/or diffusive) bottom boundary layer scheme
 
-                            CALL     trc_ini    ! Passive tracers
+      !                                     ! Passive tracers
+                            CALL     trc_init   ! Passive tracers initialization
+      !                                     ! Dynamics
                             CALL dta_dyn_init   ! Initialization for the dynamics
                             CALL     iom_init       ! iom_put initialization
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 2082)
@@ -67,4 +67,5 @@
       REAL(wp), DIMENSION(jpi,jpj    ) ::   zarea_ssh, zbotpres
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zrhd, zrhop
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::   ztsn
       !!--------------------------------------------------------------------
 
@@ -82,5 +83,7 @@
 
       !                                         ! thermosteric ssh
-      CALL eos( tn, sn0, zrhd )                       ! now in situ density using initial salinity
+      ztsn(:,:,:,jp_tem) = tn (:,:,:)
+      ztsn(:,:,:,jp_sal) = sn0(:,:,:)
+      CALL eos( ztsn, zrhd )                       ! now in situ density using initial salinity
       !
       zbotpres(:,:) = 0.e0                            ! no atmospheric surface pressure, levitating sea-ice
@@ -96,5 +99,5 @@
       
       !                                         ! steric sea surface height
-      CALL eos( tn, sn, zrhd, zrhop )                 ! now in situ and potential density
+      CALL eos( tsn, zrhd, zrhop )                 ! now in situ and potential density
       zrhop(:,:,jpk) = 0.e0
       CALL iom_put( 'rhop', zrhop )
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/IOM/restart.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/IOM/restart.F90	(revision 2082)
@@ -26,5 +26,7 @@
    USE zdfmxl          ! mixed layer depth
    USE trdmld_oce      ! ocean active mixed layer tracers trends variables
-
+#if defined key_zdfkpp
+   USE traswap
+#endif
    IMPLICIT NONE
    PRIVATE
@@ -205,5 +207,6 @@
          CALL iom_get( numror, jpdom_autoglo, 'rhd' , rhd  )       ! now    in situ density anomaly
       ELSE
-         CALL eos( tn, sn, rhd )   ! compute rhd
+         CALL tra_swap
+         CALL eos( tsn, rhd )   ! compute rhd
       ENDIF
 #endif
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 2082)
@@ -68,5 +68,5 @@
 CONTAINS
 
-   SUBROUTINE eos_insitu( ptem, psal, prd )
+   SUBROUTINE eos_insitu( pts, prd )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE eos_insitu  ***
@@ -98,5 +98,5 @@
       !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
       !!      Note that no boundary condition problem occurs in this routine
-      !!      as (ptem,psal) are defined over the whole domain.
+      !!      as pts are defined over the whole domain.
       !!
       !! ** Action  :   compute prd , the in situ density (no units)
@@ -104,7 +104,7 @@
       !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptem   ! potential temperature  [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   psal   ! salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   prd    ! in situ density 
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      !                                                               ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   prd   ! in situ density 
       !!
       INTEGER  ::   ji, jj, jk           ! dummy loop indices
@@ -123,11 +123,11 @@
          zrau0r = 1.e0 / rau0
 !CDIR NOVERRCHK
-         zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
+         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
          !  
          DO jk = 1, jpkm1
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zt = ptem  (ji,jj,jk)
-                  zs = psal  (ji,jj,jk)
+                  zt = pts   (ji,jj,jk,jp_tem)
+                  zs = pts   (ji,jj,jk,jp_sal)
                   zh = fsdept(ji,jj,jk)        ! depth
                   zsr= zws   (ji,jj,jk)        ! square root salinity
@@ -169,10 +169,10 @@
       CASE( 1 )                !==  Linear formulation function of temperature only  ==!
          DO jk = 1, jpkm1
-            prd(:,:,jk) = ( 0.0285 - rn_alpha * ptem(:,:,jk) ) * tmask(:,:,jk)
+            prd(:,:,jk) = ( 0.0285 - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
          END DO
          !
       CASE( 2 )                !==  Linear formulation function of temperature and salinity  ==!
          DO jk = 1, jpkm1
-            prd(:,:,jk) = ( rn_beta  * psal(:,:,jk) - rn_alpha * ptem(:,:,jk) ) * tmask(:,:,jk)
+            prd(:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
          END DO
          !
@@ -184,5 +184,5 @@
 
 
-   SUBROUTINE eos_insitu_pot( ptem, psal, prd, prhop )
+   SUBROUTINE eos_insitu_pot( pts, prd, prhop )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_pot  ***
@@ -230,8 +230,8 @@
       !!                Brown and Campana, Mon. Weather Rev., 1978
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptem   ! potential temperature  [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   psal   ! salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   prd    ! in situ density 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      !                                                               ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
 
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -246,11 +246,11 @@
          zrau0r = 1.e0 / rau0
 !CDIR NOVERRCHK
-         zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
+         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
          !  
          DO jk = 1, jpkm1
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zt = ptem  (ji,jj,jk)
-                  zs = psal  (ji,jj,jk)
+                  zt = pts   (ji,jj,jk,jp_tem)
+                  zs = pts   (ji,jj,jk,jp_sal)
                   zh = fsdept(ji,jj,jk)        ! depth
                   zsr= zws   (ji,jj,jk)        ! square root salinity
@@ -295,12 +295,12 @@
       CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
          DO jk = 1, jpkm1
-            prd  (:,:,jk) = ( 0.0285 - rn_alpha * ptem(:,:,jk) )        * tmask(:,:,jk)
-            prhop(:,:,jk) = ( 1.e0   +          prd (:,:,jk) ) * rau0 * tmask(:,:,jk)
+            prd  (:,:,jk) = ( 0.0285 - rn_alpha * pts(:,:,jk,jp_sal) )        * tmask(:,:,jk)
+            prhop(:,:,jk) = ( 1.e0   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk)
          END DO
          !
       CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
          DO jk = 1, jpkm1
-            prd  (:,:,jk) = ( rn_beta  * psal(:,:,jk) - rn_alpha * ptem(:,:,jk) )        * tmask(:,:,jk)
-            prhop(:,:,jk) = ( 1.e0   + prd (:,:,jk) )                         * rau0 * tmask(:,:,jk)
+            prd  (:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk)
+            prhop(:,:,jk) = ( 1.e0     + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk)
          END DO
          !
@@ -312,5 +312,5 @@
 
 
-   SUBROUTINE eos_insitu_2d( ptem, psal, pdep, prd )
+   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_2d  ***
@@ -342,5 +342,5 @@
       !!              prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
       !!      Note that no boundary condition problem occurs in this routine
-      !!      as (ptem,psal) are defined over the whole domain.
+      !!      as pts are defined over the whole domain.
       !!
       !! ** Action  : - prd , the in situ density (no units)
@@ -348,8 +348,8 @@
       !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature  [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pdep   ! depth                  [m]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   prd    ! in situ density 
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      !                                                           ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                  [m]
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
       !!
       INTEGER  ::   ji, jj                    ! dummy loop indices
@@ -369,5 +369,5 @@
 !CDIR NOVERRCHK
             DO ji = 1, fs_jpim1   ! vector opt.
-               zws(ji,jj) = SQRT( ABS( psal(ji,jj) ) )
+               zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )
             END DO
          END DO
@@ -375,6 +375,6 @@
             DO ji = 1, fs_jpim1   ! vector opt.
                zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask
-               zt    = ptem (ji,jj)            ! interpolated T
-               zs    = psal (ji,jj)            ! interpolated S
+               zt    = pts  (ji,jj,jp_tem)            ! interpolated T
+               zs    = pts  (ji,jj,jp_sal)            ! interpolated S
                zsr   = zws  (ji,jj)            ! square root of interpolated S
                zh    = pdep (ji,jj)            ! depth at the partial step level
@@ -416,5 +416,5 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               prd(ji,jj) = ( 0.0285 - rn_alpha * ptem(ji,jj) ) * tmask(ji,jj,1)
+               prd(ji,jj) = ( 0.0285 - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)
             END DO
          END DO
@@ -423,5 +423,5 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               prd(ji,jj) = ( rn_beta * psal(ji,jj) - rn_alpha * ptem(ji,jj) ) * tmask(ji,jj,1) 
+               prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
             END DO
          END DO
@@ -434,5 +434,5 @@
 
 
-   SUBROUTINE eos_bn2( ptem, psal, pn2 )
+   SUBROUTINE eos_bn2( pts, pn2 )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_bn2  ***
@@ -464,7 +464,7 @@
       !! References :   McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987.
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptem   ! potential temperature   [Celcius]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   psal   ! salinity                [psu]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pn2    ! Brunt-Vaisala frequency [s-1]
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      !                                                               ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pn2    ! Brunt-Vaisala frequency [s-1]
       !!
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -485,6 +485,6 @@
                DO ji = 1, jpi
                   zgde3w = grav / fse3w(ji,jj,jk)
-                  zt = 0.5 * ( ptem(ji,jj,jk) + ptem(ji,jj,jk-1) )          ! potential temperature at w-point
-                  zs = 0.5 * ( psal(ji,jj,jk) + psal(ji,jj,jk-1) ) - 35.0   ! salinity anomaly (s-35) at w-point
+                  zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )          ! potential temperature at w-point
+                  zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0   ! salinity anomaly (s-35) at w-point
                   zh = fsdepw(ji,jj,jk)                                     ! depth in meters  at w-point
                   !
@@ -515,11 +515,11 @@
                      !
                   pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
-                     &          * ( zalbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) )   &
-                     &                     - ( psal(ji,jj,jk-1) - psal(ji,jj,jk) ) )
+                     &          * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   &
+                     &                     - ( pts(ji,jj,jk-1njp_sal) - pts(ji,jj,jk,jp_sal) ) )
 #if defined key_zdfddm
                   !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!!
-                  zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s])
+                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s])
                   IF ( ABS( zds) <= 1.e-20 ) zds = 1.e-20
-                  rrau(ji,jj,jk) = zalbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds
+                  rrau(ji,jj,jk) = zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
 #endif
                END DO
@@ -529,11 +529,11 @@
       CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
          DO jk = 2, jpkm1
-            pn2(:,:,jk) = grav * rn_alpha * ( ptem(:,:,jk-1) - ptem(:,:,jk) ) / fse3w(:,:,jk) * tmask(:,:,jk)
+            pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk)
          END DO
          !
       CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==!
          DO jk = 2, jpkm1
-            pn2(:,:,jk) = grav * (  rn_alpha * ( ptem(:,:,jk-1) - ptem(:,:,jk) )      &
-               &                  - rn_beta  * ( psal(:,:,jk-1) - psal(:,:,jk) )  )   &
+            pn2(:,:,jk) = grav * (  rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) )      &
+               &                  - rn_beta  * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) )  )   &
                &               / fse3w(:,:,jk) * tmask(:,:,jk)
          END DO 
@@ -542,7 +542,7 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) )
+                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  
                   IF ( ABS( zds ) <= 1.e-20 ) zds = 1.e-20
-                  rrau(ji,jj,jk) = ralpbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds
+                  rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
                END DO
             END DO
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv.F90	(revision 2082)
@@ -42,4 +42,7 @@
    INTEGER ::   nadv   ! choice of the type of advection scheme
 
+   REAL(wp), DIMENSION(jpk) ::   r2dt  ! vertical profile time-step, = 2 rdttra
+      !                                ! except at nit000 (=rdttra) if neuler=0
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -66,4 +69,10 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk)   ::  zun, zvn, zwn   ! effective transport
       !!----------------------------------------------------------------------
+      !                                          ! set time step
+      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
+         r2dt(:) =  rdttra(:)                          ! = rdtra (restarting with Euler time stepping)
+      ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1
+         r2dt(:) = 2. * rdttra(:)                      ! = 2 rdttra (leapfrog)
+      ENDIF
       !                                               !==  effective transport  ==!
       DO jk = 1, jpkm1
@@ -82,28 +91,28 @@
 
       SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==!
-      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  2nd order centered
-      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD 
-      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, 'TRA', zun, zvn, zwn, tsb,      tsa, jpts )   !  MUSCL 
-      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  MUSCL2 
-      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS 
-      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST 
+      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, 'TRA',       zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  2nd order centered
+      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD 
+      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb,      tsa, jpts )   !  MUSCL 
+      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  MUSCL2 
+      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS 
+      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST 
       !
       CASE (-1 )                                      !==  esopa: test all possibility with control print  ==!
-         CALL tra_adv_cen2  ( kt , 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_cen2  ( kt, 'TRA',       zun, zvn, zwn, tsb, tsn, tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv0 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_tvd   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_tvd   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv1 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_muscl ( kt, 'TRA', zun, zvn, zwn, tsb,      tsa, jpts )          
+         CALL tra_adv_muscl ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb,      tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv3 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_muscl2( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_muscl2( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv4 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_ubs   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_ubs   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv5 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_adv_qck   ( kt, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_qck   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv6 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
@@ -179,5 +188,4 @@
          IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used'
          IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used'
-         IF( nadv ==  7 )   WRITE(numout,*) '         SMOLAR    scheme is used'
          IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme'
       ENDIF
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_cen2.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_cen2.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_cen2.F90	(revision 2082)
@@ -20,5 +20,5 @@
    !!                  area (set for orca 2 and 4 only)
    !!----------------------------------------------------------------------
-   USE oce, ONLY: tn, sn  ! now ocean temperature and salinity
+   USE oce, ONLY: tsn  ! now ocean temperature and salinity
    USE dom_oce         ! ocean space and time domain
    USE eosbn2          ! equation of state
@@ -32,4 +32,5 @@
    USE zdf_oce         ! ocean vertical physics
    USE restart         ! ocean restart
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -134,7 +135,8 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
+         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
          IF(lwp) WRITE(numout,*)
@@ -159,11 +161,11 @@
 !!gm  not strickly exact : the freezing point should be computed at each ocean levels...
 !!gm  not a big deal since cen2 is no more used in global ice-ocean simulations
-      ztfreez(:,:) = tfreez( sn(:,:,1) )
+      ztfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) )
       DO jk = 1, jpk
          DO jj = 1, jpj
             DO ji = 1, jpi
                !                                        ! below ice covered area (if tn < "freezing"+0.1 )
-               IF( tn(ji,jj,jk) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0
-               ELSE                                              ;   zice = 0.e0
+               IF( tsn(ji,jj,jk,jp_tem) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0
+               ELSE                                                      ;   zice = 0.e0
                ENDIF
                zind(ji,jj,jk) = MAX (   &
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_eiv.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_eiv.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_eiv.F90	(revision 2082)
@@ -22,4 +22,5 @@
    USE in_out_manager  ! I/O manager
    USE iom
+   USE trc_oce         ! share passive tracers/Ocean variables
 # if defined key_diaeiv
    USE phycst          ! physical constants
@@ -80,7 +81,7 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection :'
+         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection on ', cdtype,' :'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
 # if defined key_diaeiv 
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_muscl.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 2082)
@@ -24,4 +24,5 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
 
@@ -44,6 +45,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_muscl( kt, cdtype, pun, pvn, pwn, &
-      &                                       ptb, pta, kjpt   )
+   SUBROUTINE tra_adv_muscl( kt, cdtype, p2dt, pun, pvn, pwn, &
+      &                                        ptb, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE tra_adv_muscl  ***
@@ -67,4 +68,5 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb           ! before and now tracer fields
@@ -75,11 +77,11 @@
       REAL(wp) ::   zv, z0v, zzwy
       REAL(wp) ::   zw, z0w
-      REAL(wp) ::   ztra, zbtr, z2, zdt, zalpha
+      REAL(wp) ::   ztra, zbtr, zdt, zalpha
       REAL(wp), DIMENSION (jpi,jpj,jpk) :: zslpx, zslpy   ! 3D workspace
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 .AND. lwp ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'tra_adv : MUSCL advection scheme'
+         WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
          WRITE(numout,*) '~~~~~~~'
          !
@@ -88,8 +90,4 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
-      ELSE                                        ;    z2 = 2.
-      ENDIF
-      !
       !                                                     ! ===========
       DO jn = 1, kjpt                                       ! tracer loop
@@ -139,5 +137,5 @@
          !                                             !-- MUSCL horizontal advective fluxes
          DO jk = 1, jpkm1                                     ! interior values
-            zdt  = z2 * rdttra(jk)
+            zdt  = p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -230,5 +228,5 @@
          !
          DO jk = 1, jpkm1                                     ! interior values
-            zdt  = z2 * rdttra(jk)
+            zdt  = p2dt(jk)
             DO jj = 2, jpjm1      
                DO ji = fs_2, fs_jpim1   ! vector opt.
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 2082)
@@ -22,4 +22,5 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
 
@@ -43,6 +44,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_muscl2( kt, cdtype, pun, pvn, pwn, &
-      &                                   ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_muscl2( kt, cdtype, p2dt, pun, pvn, pwn, &
+      &                                         ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_adv_muscl2  ***
@@ -67,4 +68,5 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
@@ -75,11 +77,11 @@
       REAL(wp) ::   zv, z0v, zzwy
       REAL(wp) ::   zw, z0w
-      REAL(wp) ::   ztra, zbtr, z2, zdt, zalpha
+      REAL(wp) ::   ztra, zbtr, zdt, zalpha
       REAL(wp), DIMENSION (jpi,jpj,jpk) ::  zslpx, zslpy   ! 3D workspace
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 .AND. lwp ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme'
+         WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme on ', cdtype
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
          !
@@ -88,8 +90,4 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2 = 1.
-      ELSE                                        ;   z2 = 2.
-      ENDIF
-      !
       !
       DO jn = 1, kjpt                                            ! tracer loop
@@ -139,5 +137,5 @@
         !                                             !-- MUSCL horizontal advective fluxes
          DO jk = 1, jpkm1                                     ! interior values
-            zdt  = z2 * rdttra(jk)
+            zdt  = p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -255,5 +253,5 @@
          !
          DO jk = 1, jpkm1                                     ! interior values
-            zdt  = z2 * rdttra(jk)
+            zdt  = p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_qck.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 2082)
@@ -25,4 +25,5 @@
    USE in_out_manager  ! I/O manager
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -45,6 +46,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_qck ( kt, cdtype, pun, pvn, pwn, &
-      &                                 ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_qck ( kt, cdtype, p2dt, pun, pvn, pwn, &
+      &                                       ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_qck  ***
@@ -85,14 +86,13 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend 
-      !!
-      REAL(wp) ::   z2                            !    temporary scalar
-      !!----------------------------------------------------------------------
-
-      IF( kt == nit000 ) THEN
+      !!----------------------------------------------------------------------
+
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme'
+         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
          IF(lwp) WRITE(numout,*)
@@ -102,13 +102,9 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
-      ELSE                                        ;    z2 = 2.
-      ENDIF
-
       ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
       !---------------------------------------------------------------------------
 
-      CALL tra_adv_qck_i( kt, cdtype, pun, z2, ptb, ptn, pta, kjpt ) 
-      CALL tra_adv_qck_j( kt, cdtype, pvn, z2, ptb, ptn, pta, kjpt ) 
+      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
+      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
 
       ! II. The vertical fluxes are computed with the 2nd order centered scheme
@@ -119,6 +115,6 @@
    END SUBROUTINE tra_adv_qck
 
-   SUBROUTINE tra_adv_qck_i( kt, cdtype, pun, pz2,   &
-      &                                  ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,    &
+      &                                        ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!
@@ -129,8 +125,8 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun             ! zonal velocity component
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn    ! before tracer fields
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend 
-      REAL(wp)        , INTENT(in   )                               ::   pz2
       !!
       INTEGER  :: ji, jj, jk, jn           ! dummy loop indices
@@ -176,5 +172,5 @@
          !
          DO jk = 1, jpkm1  
-            zdt =  pz2 * rdttra(jk)
+            zdt =  p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.   
@@ -239,6 +235,6 @@
    END SUBROUTINE tra_adv_qck_i
 
-   SUBROUTINE tra_adv_qck_j( kt, cdtype, pvn, pz2,   &
-      &                                  ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,   &
+      &                                        ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!
@@ -250,8 +246,8 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pvn             ! meridional velocity component
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn    ! before tracer fields
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend 
-      REAL(wp)        , INTENT(in   )                               ::   pz2
       !!
       INTEGER  :: ji, jj, jk, jn           ! dummy loop indices
@@ -296,5 +292,5 @@
          !
          DO jk = 1, jpkm1  
-            zdt =  pz2 * rdttra(jk)
+            zdt =  p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.   
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 2082)
@@ -32,4 +32,5 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
 
@@ -52,6 +53,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_tvd ( kt, cdtype, pun, pvn, pwn, &
-      &                                 ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_tvd ( kt, cdtype, p2dt, pun, pvn, pwn, &
+      &                                       ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_tvd  ***
@@ -73,4 +74,5 @@
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
@@ -79,5 +81,5 @@
       INTEGER  ::   ji, jj, jk, jn          ! dummy loop indices  
       REAL(wp) ::   &
-         z2, z2dtt, zbtr, ztra,              &  ! temporary scalar
+         z2dtt, zbtr, ztra,                  &  ! temporary scalar
          zfp_ui, zfp_vj, zfp_wk,             &  !    "         "
          zfm_ui, zfm_vj, zfm_wk                 !    "         "
@@ -86,9 +88,8 @@
       !!----------------------------------------------------------------------
 
-      zwi(:,:,:) = 0.e0
-
-      IF( kt == nit000 .AND. lwp ) THEN
+
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme'
+         WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype
          WRITE(numout,*) '~~~~~~~~~~~'
          !
@@ -103,7 +104,5 @@
       END IF
       !
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
-      ELSE                                        ;    z2 = 2.
-      ENDIF
+      zwi(:,:,:) = 0.e0
       !
       !                                                          ! ===========
@@ -150,5 +149,5 @@
          ! total advective trend
          DO jk = 1, jpkm1
-            z2dtt = z2 * rdttra(jk)
+            z2dtt = p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -209,5 +208,5 @@
          ! 4. monotonicity algorithm
          ! -------------------------
-         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, z2 )
+         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
 
 
@@ -253,5 +252,5 @@
 
 
-   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, prdt )
+   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc  ***
@@ -266,5 +265,6 @@
       !!       in-space based differencing for fluid
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT( in ) ::   prdt   ! ???
+      REAL(wp), DIMENSION(jpk)         , INTENT( in    ) ::   &  
+         p2dt                               ! vertical profile of tracer time-step
       REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( in    ) ::   &
          pbef,                            & ! before field
@@ -299,5 +299,5 @@
       DO jk = 1, jpkm1
          ikm1 = MAX(jk-1,1)
-         z2dtt = prdt * rdttra(jk)
+         z2dtt = p2dt(jk)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_ubs.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 2082)
@@ -21,4 +21,5 @@
    USE diaptr          ! poleward transport diagnostics
    USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -40,6 +41,6 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_ubs ( kt, cdtype, pun, pvn, pwn, &
-      &                                 ptb, ptn, pta, kjpt   )
+   SUBROUTINE tra_adv_ubs ( kt, cdtype, p2dt, pun, pvn, pwn, &
+      &                                       ptb, ptn, pta, kjpt   )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_ubs  ***
@@ -81,4 +82,5 @@
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
       INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
+      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend 
@@ -88,5 +90,5 @@
       REAL(wp) ::   zfp_ui, zfm_ui, zcenut  !    "         "
       REAL(wp) ::   zfp_vj, zfm_vj, zcenvt  !    "         "    !    "         "
-      REAL(wp) ::   z2dtt, z2                   
+      REAL(wp) ::   z2dtt                   
       REAL(wp) ::   ztak, zfp_wk, zfm_wk    !    "         "
       REAL(wp) ::   zeeu, zeev, z_hdivn     
@@ -96,15 +98,11 @@
 
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme'
+         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
          !
          l_trd = .FALSE.
          IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
-      ENDIF
-      !
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
-      ELSE                                        ;    z2 = 2.
       ENDIF
       !
@@ -232,5 +230,5 @@
          ! update and guess with monotonic sheme
          DO jk = 1, jpkm1
-            z2dtt = z2 * rdttra(jk)
+            z2dtt = p2dt(jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -255,5 +253,5 @@
          END DO
          !
-         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, z2 )      !  monotonicity algorithm
+         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
 
          !  final trend with corrected fluxes
@@ -288,5 +286,5 @@
    END SUBROUTINE tra_adv_ubs
 
-   SUBROUTINE nonosc_z( pbef, pcc, paft, prdt )
+   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc_z  ***
@@ -301,5 +299,5 @@
       !!       in-space based differencing for fluid
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   )                          ::   prdt   ! ???
+      REAL(wp), INTENT(in   ), DIMENSION(jpk)          ::   p2dt            ! vertical profile of tracer time-step
       REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
       REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
@@ -356,5 +354,5 @@
 
       DO jk = 1, jpkm1
-         z2dtt = prdt * rdttra(jk)
+         z2dtt = p2dt(jk)
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -374,5 +372,4 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-
                za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
                zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 2082)
@@ -32,4 +32,5 @@
    USE lbclnk         ! ocean lateral boundary conditions
    USE prtctl         ! Print control
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -377,7 +378,7 @@
       
       !                                        !* bottom temperature, salinity, velocity and depth
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp)  WRITE(numout,*) ' '
-         IF(lwp)  WRITE(numout,*) ' trabbl:bbl  : Compute bbl velocities and diffusive coefficients in ', cdtype, ' at time step ', kt
+         IF(lwp)  WRITE(numout,*) ' trabbl:bbl  : Compute bbl velocities and diffusive coefficients in ', cdtype
          IF(lwp)  WRITE(numout,*) ' '
       ENDIF
@@ -391,6 +392,6 @@
 #endif
             ik = mbkt(ji,jj)                        ! bottom T-level index
-            ztb (ji,jj) = tb(ji,jj,ik)              ! bottom before T and S
-            zsb (ji,jj) = sb(ji,jj,ik)
+            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem)              ! bottom before T and S
+            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal)
             zdep(ji,jj) = fsdept_0(ji,jj,ik)        ! bottom T-level reference depth
             !
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf.F90	(revision 2082)
@@ -62,5 +62,4 @@
       !!
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ztrdt, ztrds
-      REAL(wp), DIMENSION(jpi,jpj,jpts) :: zgtsu, zgtsv
       !!----------------------------------------------------------------------
 
@@ -70,24 +69,21 @@
       ENDIF
 
-      zgtsu(:,:,jp_tem) = gtu(:,:)        ;        zgtsu(:,:,jp_sal) = gsu(:,:)
-      zgtsv(:,:,jp_tem) = gtv(:,:)        ;        zgtsv(:,:,jp_sal) = gsv(:,:)
-
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts        )  ! iso-level laplacian
-      CASE ( 1 )   ;   CALL tra_ldf_iso   ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts, ahtb0 )  ! rotated laplacian 
-      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts        )  ! iso-level bilaplacian
-      CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt , 'TRA',               tsb, tsa, jpts        )  ! s-coord. horizontal bilaplacian
+      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts        )  ! iso-level laplacian
+      CASE ( 1 )   ;   CALL tra_ldf_iso   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )  ! rotated laplacian 
+      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts        )  ! iso-level bilaplacian
+      CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt, 'TRA',             tsb, tsa, jpts        )  ! s-coord. horizontal bilaplacian
          !
       CASE ( -1 )                                     ! esopa: test all possibility with control print
-         CALL tra_ldf_lap   ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts        ) 
+         CALL tra_ldf_lap   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts        ) 
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf0 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_ldf_iso   ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts, ahtb0 )  
+         CALL tra_ldf_iso   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )  
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf1 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_ldf_bilap ( kt , 'TRA', zgtsu, zgtsv, tsb, tsa, jpts        ) 
+         CALL tra_ldf_bilap ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts        ) 
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf2 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_ldf_bilapg( kt , 'TRA',               tsb, tsa, jpts        ) 
+         CALL tra_ldf_bilapg( kt, 'TRA',             tsb, tsa, jpts        ) 
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf3 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_bilap.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 2082)
@@ -28,4 +28,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -93,7 +94,7 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator'
+         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
       ENDIF
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_bilapg.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_bilapg.F90	(revision 2082)
@@ -24,4 +24,5 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link)
    USE diaptr          ! poleward transport diagnostics 
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -80,7 +81,7 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
+         IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
       ENDIF
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 2082)
@@ -28,4 +28,5 @@
    USE ldfslp          ! iso-neutral slopes
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 #if defined key_diaar5
    USE phycst          ! physical constants
@@ -116,7 +117,7 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator'
+         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_lap.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 2082)
@@ -23,4 +23,5 @@
    USE in_out_manager  ! I/O manager
    USE diaptr          ! poleward transport diagnostics
+   USE trc_oce         ! share passive tracers/Ocean variables
 
 
@@ -84,7 +85,7 @@
       !!----------------------------------------------------------------------
       
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion'
+         IF(lwp) WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
          e1ur(:,:) = e2u(:,:) / e1u(:,:)
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/tranpc.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/tranpc.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/tranpc.F90	(revision 2082)
@@ -74,5 +74,5 @@
          inpci = 0
 
-         CALL eos( tsa(:,:,:,jp_tem), tsa(:,:,:,jp_sal), rhd, zrhop )         ! Potential density
+         CALL eos( tsa, rhd, zrhop )         ! Potential density
 
          IF( l_trdtra )   THEN                    !* Save ta and sa trends
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2082)
@@ -48,5 +48,5 @@
    PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
 
-   REAL(wp), DIMENSION(jpk) ::   r2dt_t   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
+   REAL(wp), DIMENSION(jpk) ::   r2dt   ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
 
    !! * Substitutions
@@ -120,6 +120,6 @@
  
       ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler)
-      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt(:) =     rdttra(:)      ! at nit000             (Euler)
+      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
       ENDIF
 
@@ -145,5 +145,5 @@
       IF( l_trdtra ) THEN              ! trend of the Asselin filter (tb filtered - tb)/dt     
          DO jk = 1, jpkm1
-            zfact = 1.e0 / r2dt_t(jk)             
+            zfact = 1.e0 / r2dt(jk)             
             ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
             ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
@@ -194,5 +194,5 @@
       !!----------------------------------------------------------------------
 
-      IF( kt == kit000 ) THEN
+      IF( kt == kit000 )  THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 2082)
@@ -79,11 +79,11 @@
 
       SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )    ;    CALL tra_zdf_exp( kt , 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme 
-      CASE ( 1 )    ;    CALL tra_zdf_imp( kt , 'TRA', r2dt,            tsb, tsa, jpts )  !   implicit scheme 
+      CASE ( 0 )    ;    CALL tra_zdf_exp( kt, 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme 
+      CASE ( 1 )    ;    CALL tra_zdf_imp( kt, 'TRA', r2dt,            tsb, tsa, jpts )  !   implicit scheme 
       CASE ( -1 )                                       ! esopa: test all possibility with control print
-         CALL tra_zdf_exp( kt , 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )
+         CALL tra_zdf_exp( kt, 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf0 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_zdf_imp( kt , 'TRA', r2dt,            tsb, tsa, jpts ) 
+         CALL tra_zdf_imp( kt, 'TRA', r2dt,            tsb, tsa, jpts ) 
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf1 - Ta: ', mask1=tmask,               &
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_exp.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_exp.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_exp.F90	(revision 2082)
@@ -29,4 +29,5 @@
    USE zdfddm          ! ocean vertical physics: double diffusion
    USE in_out_manager  ! I/O manager
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -87,7 +88,7 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_zdf_exp : explicit vertical mixing'
+         IF(lwp) WRITE(numout,*) 'tra_zdf_exp : explicit vertical mixing on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_imp.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 2082)
@@ -31,4 +31,5 @@
    USE domvvl          ! variable volume
    USE ldftra          ! lateral mixing type
+   USE trc_oce         ! share passive tracers/Ocean variables
 
    IMPLICIT NONE
@@ -105,7 +106,7 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
          IF(lwp)WRITE(numout,*)
-         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing'
+         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
          IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
          zavi = 0.e0      ! avoid warning at compilation phase when lk_ldfslp=F
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 2082)
@@ -14,6 +14,4 @@
    !!   zps_hde      :  Horizontal DErivative of T, S and rd at the last
    !!                   ocean level (Z-coord. with Partial Steps)
-   !!   zps_hde_trc  :  Horizontal DErivative of passive tracers at the last
-   !!                   ocean level (Z-coord. with Partial Steps)
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -31,7 +29,4 @@
    PUBLIC zps_hde          ! routine called by step.F90
    PUBLIC zps_hde_init     ! routine called by opa.F90
-#if defined key_top
-   PUBLIC zps_hde_trc 
-#endif
 
    !! * module variables
@@ -49,7 +44,7 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE zps_hde ( kt, ptem, psal, prd ,   &
-                            pgtu, pgsu, pgru,   &
-                            pgtv, pgsv, pgrv  )
+
+   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   &
+                                 prd, pgru, pgrv    )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE zps_hde  ***
@@ -91,159 +86,27 @@
       !!          di(rho) = rd~ - rd(i,j,k) or rd (i+1,j,k) - rd~
       !!
-      !! ** Action  : - pgtu, pgsu, pgru: horizontal gradient of T, S
-      !!                and rd at U-points 
-      !!              - pgtv, pgsv, pgrv: horizontal gradient of T, S
+      !! ** Action  : - pgtu, pgtv: horizontal gradient of tracer at U/V-points
+      !!              - pgru, pgrv: horizontal gradient of rd if present at U/V-points 
       !!                and rd at V-points 
       !!----------------------------------------------------------------------
       !! * Arguments
-      INTEGER, INTENT( in ) ::   kt ! ocean time-step index
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in )  :: ptem, psal, prd  ! 3D T, S and rd fields
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT( out ) :: pgtu, pgsu, pgru !  horizontal grad. of T, S and rd at u-point 
-      REAL(wp), DIMENSION(jpi,jpj)    , INTENT( out ) :: pgtv, pgsv, pgrv !  horizontal grad. of T, S and rd at v-point 
-      !! * Local declarations
-      INTEGER ::   ji , jj          ! Dummy loop indices
-      INTEGER ::   iku, ikv         ! partial step level at u- and v-points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zti, ztj, zsi, zsj   ! interpolated value of T, S
-      REAL(wp), DIMENSION(jpi,jpj) ::   zri, zrj             ! interpolated value of rd
-      REAL(wp), DIMENSION(jpi,jpj) ::   zhgi, zhgj           ! depth of interpolation for eos2d
-      REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv                 ! temporary scalars
-
-
-      ! Interpolation of T and S at the last ocean level
-# if defined key_vectopt_loop
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
-# else
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-# endif
-            ! last level
-            iku = mbatu(ji,jj)
-            ikv = mbatv(ji,jj)
-            ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
-            ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
-
-            ! i- direction
-            IF( ze3wu >= 0. ) THEN      ! case 1
-               ! interpolated values of T and S
-               zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
-               zti(ji,jj) = ptem(ji+1,jj,iku) + zmaxu * ( ptem(ji+1,jj,iku-1) - ptem(ji+1,jj,iku) )
-               zsi(ji,jj) = psal(ji+1,jj,iku) + zmaxu * ( psal(ji+1,jj,iku-1) - psal(ji+1,jj,iku) )
-               ! depth of the partial step level
-               zhgi(ji,jj) = fsdept(ji,jj,iku)
-               ! gradient of T and S
-               pgtu(ji,jj) = umask(ji,jj,1) * ( zti(ji,jj) - ptem(ji,jj,iku) )
-               pgsu(ji,jj) = umask(ji,jj,1) * ( zsi(ji,jj) - psal(ji,jj,iku) )
-
-            ELSE                        ! case 2
-               ! interpolated values of T and S
-               zmaxu =  -ze3wu / fse3w(ji,jj,iku)
-               zti(ji,jj) = ptem(ji,jj,iku) + zmaxu * ( ptem(ji,jj,iku-1) - ptem(ji,jj,iku) )
-               zsi(ji,jj) = psal(ji,jj,iku) + zmaxu * ( psal(ji,jj,iku-1) - psal(ji,jj,iku) )
-               ! depth of the partial step level
-               zhgi(ji,jj) = fsdept(ji+1,jj,iku)
-               ! gradient of T and S 
-               pgtu(ji,jj) = umask(ji,jj,1) * ( ptem(ji+1,jj,iku) - zti (ji,jj) )
-               pgsu(ji,jj) = umask(ji,jj,1) * ( psal(ji+1,jj,iku) - zsi (ji,jj) )
-            ENDIF
-
-            ! j- direction
-            IF( ze3wv >= 0. ) THEN      ! case 1
-               ! interpolated values of T and S
-               zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
-               ztj(ji,jj) = ptem(ji,jj+1,ikv) + zmaxv * ( ptem(ji,jj+1,ikv-1) - ptem(ji,jj+1,ikv) )
-               zsj(ji,jj) = psal(ji,jj+1,ikv) + zmaxv * ( psal(ji,jj+1,ikv-1) - psal(ji,jj+1,ikv) )
-               ! depth of the partial step level
-               zhgj(ji,jj) = fsdept(ji,jj,ikv) 
-               ! gradient of T and S
-               pgtv(ji,jj) = vmask(ji,jj,1) * ( ztj(ji,jj) - ptem(ji,jj,ikv) )
-               pgsv(ji,jj) = vmask(ji,jj,1) * ( zsj(ji,jj) - psal(ji,jj,ikv) )
-
-            ELSE                        ! case 2
-               ! interpolated values of T and S
-               zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
-               ztj(ji,jj) = ptem(ji,jj,ikv) + zmaxv * ( ptem(ji,jj,ikv-1) - ptem(ji,jj,ikv) )
-               zsj(ji,jj) = psal(ji,jj,ikv) + zmaxv * ( psal(ji,jj,ikv-1) - psal(ji,jj,ikv) ) 
-               ! depth of the partial step level
-               zhgj(ji,jj) = fsdept(ji,jj+1,ikv) 
-               ! gradient of T and S
-               pgtv(ji,jj) = vmask(ji,jj,1) * ( ptem(ji,jj+1,ikv) - ztj(ji,jj) )
-               pgsv(ji,jj) = vmask(ji,jj,1) * ( psal(ji,jj+1,ikv) - zsj(ji,jj) )
-            ENDIF
-# if ! defined key_vectopt_loop
-         END DO
-# endif
-      END DO
-
-      ! Compute interpolated rd from zti, zsi, ztj, zsj for the 2 cases at the depth of the partial
-      ! step and store it in  zri, zrj for each  case
-      CALL eos( zti, zsi, zhgi, zri )
-      CALL eos( ztj, zsj, zhgj, zrj )
-
-
-      ! Gradient of density at the last level 
-# if defined key_vectopt_loop
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
-# else
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-# endif
-            iku = mbatu(ji,jj)
-            ikv = mbatv(ji,jj)
-            ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
-            ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
-            IF( ze3wu >= 0. ) THEN    ! i-direction: case 1
-               pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji,jj) - prd(ji,jj,iku) )
-            ELSE                      ! i-direction: case 2
-               pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )
-            ENDIF
-            IF( ze3wv >= 0. ) THEN    ! j-direction: case 1
-               pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj) - prd(ji,jj,ikv) )  
-            ELSE                      ! j-direction: case 2
-               pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )
-            ENDIF
-# if ! defined key_vectopt_loop
-         END DO
-# endif
-      END DO
-
-      ! Lateral boundary conditions on each gradient
-      CALL lbc_lnk( pgtu , 'U', -1. )   ;   CALL lbc_lnk( pgtv , 'V', -1. )
-      CALL lbc_lnk( pgsu , 'U', -1. )   ;   CALL lbc_lnk( pgsv , 'V', -1. )
-      CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )
-
-   END SUBROUTINE zps_hde
-
-#if defined key_top
-   !!----------------------------------------------------------------------
-   !!   'key_top'                                                TOP models
-   !!----------------------------------------------------------------------
-   SUBROUTINE zps_hde_trc ( kt, kjpt, ptra, pgtru, pgtrv )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE zps_hde_trc  ***
-      !!                    
-      !! ** Purpose :   Compute the horizontal derivative of passive tracers
-      !!      TRA at u- and v-points with a linear interpolation for z-coordinate
-      !!      with partial steps.
-      !!
-      !! ** Method  :   the same for T & S
-      !!
-      !! ** Action  : - pgtru : horizontal gradient of TRA at U-points 
-      !!              - pgtrv : horizontal gradient of TRA at V-points 
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER                              , INTENT( in )  ::  kt    ! ocean time-step index
-      INTEGER                              , INTENT( in )  ::  kjpt  ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT( in )  ::  ptra  ! 4D tracers fields
-      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT( out ) ::  pgtru, pgtrv  ! horizontal grad. of TRA u- and v-points 
+      INTEGER                              , INTENT( in )           ::  kt    ! ocean time-step index
+      INTEGER                              , INTENT( in )           ::  kjpt  ! number of tracers
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT( in )           ::  pta   ! 4D active or passive tracers fields
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT( out)           ::  pgtu, pgtv  ! horizontal grad. of ptra u- and v-points 
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT( in ), OPTIONAL ::  prd   ! 3D rd fields
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT( out), OPTIONAL ::  pgru, pgrv  ! horizontal grad. of prd u- and v-points 
       !! * Local declarations
       INTEGER  ::   ji, jj, jn      ! Dummy loop indices
       INTEGER  ::   iku, ikv        ! partial step level at u- and v-points
-      REAL(wp) ::   ztrai, ztraj, ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars
-      !!----------------------------------------------------------------------
-
+      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::   zti, ztj     ! interpolated value of tracer
+      REAL(wp), DIMENSION(jpi,jpj)      ::   zri, zrj     ! interpolated value of rd
+      REAL(wp), DIMENSION(jpi,jpj)      ::   zhi, zhj     ! depth of interpolation for eos2d
+      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars
+      !!----------------------------------------------------------------------
+
+
+      ! Interpolation of tracers at the last ocean level
       DO jn = 1, kjpt
-         ! Interpolation of passive tracers at the last ocean level
 # if defined key_vectopt_loop
          jj = 1
@@ -262,14 +125,14 @@
                IF( ze3wu >= 0. ) THEN      ! case 1
                   zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
-                  ! interpolated values of passive tracers
-                  ztrai = ptra(ji+1,jj,iku,jn) + zmaxu * ( ptra(ji+1,jj,iku-1,jn) - ptra(ji+1,jj,iku,jn) )
-                  ! gradient of passive tracers
-                  pgtru(ji,jj,jn) = umask(ji,jj,1) * ( ztrai - ptra(ji,jj,iku,jn) )
+                  ! interpolated values of tracers
+                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku-1,jn) - pta(ji+1,jj,iku,jn) )
+                  ! gradient of  tracers
+                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
                ELSE                        ! case 2
                   zmaxu = -ze3wu / fse3w(ji,jj,iku)
-                  ! interpolated values of passive tracers
-                  ztrai = ptra(ji,jj,iku,jn) + zmaxu * ( ptra(ji,jj,iku-1,jn) - ptra(ji,jj,iku,jn) )
-                  ! gradient of passive tracers
-                  pgtru(ji,jj,jn) = umask(ji,jj,1) * ( ptra(ji+1,jj,iku,jn) - ztrai )
+                  ! interpolated values of tracers
+                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku-1,jn) - pta(ji,jj,iku,jn) )
+                  ! gradient of tracers
+                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
                ENDIF
 
@@ -277,14 +140,14 @@
                IF( ze3wv >= 0. ) THEN      ! case 1
                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
-                  ! interpolated values of passive tracers
-                  ztraj = ptra(ji,jj+1,ikv,jn) + zmaxv * ( ptra(ji,jj+1,ikv-1,jn) - ptra(ji,jj+1,ikv,jn) )
-                  ! gradient of passive tracers
-                  pgtrv(ji,jj,jn) = vmask(ji,jj,1) * ( ztraj - ptra(ji,jj,ikv,jn) )
+                  ! interpolated values of tracers
+                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv-1,jn) - pta(ji,jj+1,ikv,jn) )
+                  ! gradient of tracers
+                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
                ELSE                        ! case 2
                   zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
-                  ! interpolated values of passive tracers
-                  ztraj = ptra(ji,jj,ikv,jn) + zmaxv * ( ptra(ji,jj,ikv-1,jn) - ptra(ji,jj,ikv,jn) )
-                  ! gradient of passive tracers
-                  pgtrv(ji,jj,jn) = vmask(ji,jj,1) * ( ptra(ji,jj+1,ikv,jn) - ztraj )
+                  ! interpolated values of tracers
+                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv-1,jn) - pta(ji,jj,ikv,jn) )
+                  ! gradient of tracers
+                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
                ENDIF
 # if ! defined key_vectopt_loop
@@ -294,11 +157,76 @@
 
          ! Lateral boundary conditions on each gradient
-         CALL lbc_lnk( pgtru(:,:,jn) , 'U', -1. )
-         CALL lbc_lnk( pgtrv(:,:,jn) , 'V', -1. )
+         CALL lbc_lnk( pgtu(:,:,jn) , 'U', -1. )
+         CALL lbc_lnk( pgtv(:,:,jn) , 'V', -1. )
 
       END DO
 
-   END SUBROUTINE zps_hde_trc
-#endif
+      ! horizontal derivative of rd
+      IF( PRESENT( prd ) ) THEN
+         ! depth of the partial step level
+# if defined key_vectopt_loop
+         jj = 1
+         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
+# else
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+# endif
+               iku = mbatu(ji,jj)
+               ikv = mbatv(ji,jj)
+               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
+               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
+               IF( ze3wu >= 0. ) THEN   
+                  zhi(ji,jj) = fsdept(ji  ,jj,iku)
+               ELSE                   
+                  zhi(ji,jj) = fsdept(ji+1,jj,iku)
+               ENDIF
+               IF( ze3wv >= 0. ) THEN  
+                  zhj(ji,jj) = fsdept(ji,jj  ,ikv)
+               ELSE                   
+                  zhj(ji,jj) = fsdept(ji,jj+1,ikv)
+               ENDIF
+# if ! defined key_vectopt_loop
+            END DO
+# endif
+         END DO
+
+         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
+         ! step and store it in  zri, zrj for each  case
+         CALL eos( zti, zhi, zri )
+         CALL eos( ztj, zhj, zrj )
+
+         ! Gradient of density at the last level 
+# if defined key_vectopt_loop
+         jj = 1
+         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
+# else
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+# endif
+               iku = mbatu(ji,jj)
+               ikv = mbatv(ji,jj)
+               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
+               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
+               IF( ze3wu >= 0. ) THEN    ! i-direction: case 1
+                  pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji,jj) - prd(ji,jj,iku) )
+               ELSE                      ! i-direction: case 2
+                  pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )
+               ENDIF
+               IF( ze3wv >= 0. ) THEN    ! j-direction: case 1
+                  pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj) - prd(ji,jj,ikv) )
+               ELSE                      ! j-direction: case 2
+                  pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )
+               ENDIF
+# if ! defined key_vectopt_loop
+            END DO
+# endif
+         END DO
+
+         ! Lateral boundary conditions on each gradient
+         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )
+         !
+      END IF
+      !
+   END SUBROUTINE zps_hde
 
    SUBROUTINE zps_hde_init
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRD/trdicp.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRD/trdicp.F90	(revision 2082)
@@ -370,5 +370,5 @@
          zkepe(:,:,:) = 0.e0
    
-         CALL eos( tn, sn, rhd, rhop )       ! now potential and in situ densities
+         CALL eos( tsn, rhd, rhop )       ! now potential and in situ densities
 
          ! Density flux at w-point
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/istate.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/istate.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/istate.F90	(revision 2082)
@@ -42,4 +42,5 @@
    USE dynspg_exp      ! pressure gradient schemes
    USE dynspg_ts       ! pressure gradient schemes
+   USE traswp          ! Swap arrays                      (tra_swp routine)
    
    IMPLICIT NONE
@@ -122,10 +123,11 @@
          ENDIF
 
-         CALL eos( tb, sb, rhd, rhop )        ! before potential and in situ densities
+         ta(:,:,:) = 0.       ;     sa(:,:,:) = 0.
+         CALL tra_swap
+         CALL eos( tsb, rhd, rhop )        ! before potential and in situ densities
          
          IF( ln_zps .AND. .NOT. lk_c1d )   &
-            &             CALL zps_hde( nit000, tb, sb, rhd,  &  ! Partial steps: before Horizontal DErivative
-            &                                  gtu, gsu, gru, &  ! of t, s, rd at the bottom ocean level
-            &                                  gtv, gsv, grv )
+                          CALL zps_hde( nit000, 'TRA', tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
+            &                                          rhd, gru , grv , jpts )      ! of t, s, rd at the last ocean level
          
       ENDIF
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/oce.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/oce.F90	(revision 2082)
@@ -39,6 +39,6 @@
    !! interpolated gradient (only used in zps case)
    !! ---------------------
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gtu, gsu, gru   !: horizontal gradient of T, S and rd at bottom u-point
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   gtv, gsv, grv   !: horizontal gradient of T, S and rd at bottom v-point 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpts) ::   gtsu, gtsv !: horizontal gradient of T, S bottom u-point
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)      ::   gru , grv   !: horizontal gradient of rd at bottom u-point
 
    !! free surface                       !  before  !  now     !  after   !
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/opa.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/opa.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/opa.F90	(revision 2082)
@@ -37,5 +37,4 @@
    USE step_oce        ! Time stepping module definition
    USE sbc_oce         ! surface boundary condition: ocean
-   USE trdmod_oce      ! ocean variables trends
    USE domcfg          ! domain configuration               (dom_cfg routine)
    USE mppini          ! shared/distributed memory setting (mpp_init routine)
@@ -48,5 +47,4 @@
    USE zdfini
    USE phycst          ! physical constant                  (par_cst routine)
-   USE trdmod          ! momentum/tracers trends       (trd_mod_init routine)
    USE step            ! OPA time-stepping                  (stp     routine)
 #if defined key_oasis3
@@ -160,7 +158,5 @@
       !!
       !!----------------------------------------------------------------------
-#if defined key_oasis3 || defined key_oasis4 || defined key_iomput
       INTEGER :: ilocal_comm
-#endif
       CHARACTER(len=80),dimension(10) ::   cltxt = ''
       INTEGER                         ::   ji   ! local loop indices
@@ -192,4 +188,5 @@
       narea = mynode( cltxt, ilocal_comm )    ! Nodes selection (control print return in cltxt)
 # else
+      ilocal_comm = 0
       narea = mynode( cltxt )                 ! Nodes selection (control print return in cltxt)
 # endif
@@ -220,7 +217,7 @@
       !                             !--------------------------------!
 
-      CALL opa_flg                          ! Control prints & Benchmark
-
-                                            ! Domain decomposition
+      CALL opa_flg                           ! Control prints & Benchmark
+
+      !                                      ! Domain decomposition
       IF( jpni*jpnj == jpnij ) THEN   ;   CALL mpp_init      ! standard cutting out
       ELSE                            ;   CALL mpp_init2     ! eliminate land processors
@@ -231,4 +228,5 @@
      
 
+      !                                      ! General initialization
                             CALL     phy_cst    ! Physical constants
                             CALL     eos_init   ! Equation of state
@@ -245,20 +243,21 @@
 
       !                                     ! Ocean physics
-                            CALL     sbc_init   ! Read namsbc namelist : surface module (needed for iom_init)
-
-      !                                     ! Vertical physics
-                            CALL     zdf_init   ! namelist read
-                            CALL zdf_bfr_init   ! bottom friction
-      IF( lk_zdfric     )   CALL zdf_ric_init   ! Richardson number dependent Kz
-      IF( lk_zdftke_old )   CALL zdf_tke_init   ! TKE closure scheme for Kz (old scheme)
-      IF( lk_zdftke     )   CALL     tke_init   ! TKE closure scheme for Kz
-      IF( lk_zdfkpp     )   CALL zdf_kpp_init   ! KPP closure scheme for Kz
-      IF( lk_zdftmx     )   CALL zdf_tmx_init   ! tidal vertical mixing
+                            CALL     sbc_init   !  Forcings : surface module 
+
+      !                                         ! Vertical physics
+                            CALL     zdf_init      ! namelist read
+                            CALL zdf_bfr_init      ! bottom friction
+      IF( lk_zdfric     )   CALL zdf_ric_init      ! Richardson number dependent Kz
+      IF( lk_zdftke_old )   CALL zdf_tke_init      ! TKE closure scheme for Kz (old scheme)
+      IF( lk_zdftke     )   CALL     tke_init      ! TKE closure scheme for Kz
+      IF( lk_zdfkpp     )   CALL zdf_kpp_init      ! KPP closure scheme for Kz
+      IF( lk_zdftmx     )   CALL zdf_tmx_init      ! tidal vertical mixing
       IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   & 
-         &                  CALL zdf_ddm_init   ! double diffusive mixing
-      !                                     ! Lateral physics
-                            CALL ldf_tra_init   ! Lateral ocean tracer physics
-                            CALL ldf_dyn_init   ! Lateral ocean momentum physics
-      IF( lk_ldfslp )       CALL ldf_slp_init   ! slope of lateral mixing
+         &                  CALL zdf_ddm_init      ! double diffusive mixing
+      !                                         ! Lateral physics
+                            CALL ldf_tra_init      ! Lateral ocean tracer physics
+                            CALL ldf_dyn_init      ! Lateral ocean momentum physics
+      IF( lk_ldfslp )       CALL ldf_slp_init      ! slope of lateral mixing
+
       !                                     ! Active tracers
                             CALL tra_qsr_init   ! penetrative solar radiation qsr
@@ -280,7 +279,6 @@
 #if defined key_top
       !                                     ! Passive tracers
-                            CALL     trc_ini
-#endif
-
+                            CALL     trc_init
+#endif
       !                                     ! diagnostics
                             CALL     iom_init   ! iom_put initialization
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/step.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/step.F90	(revision 2082)
@@ -105,6 +105,6 @@
       ! Ocean physics update                (ua, va, ta, sa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL bn2( tb, sb, rn2b )     ! before Brunt-Vaisala frequency
-                         CALL bn2( tn, sn, rn2  )     ! now    Brunt-Vaisala frequency
+                         CALL bn2( tsb, rn2b )        ! before Brunt-Vaisala frequency
+                         CALL bn2( tsn, rn2  )        ! now    Brunt-Vaisala frequency
       !
       !  VERTICAL PHYSICS   
@@ -137,8 +137,7 @@
       !
       IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
-                         CALL eos( tb, sb, rhd )                ! before in situ density
-         IF( ln_zps )    CALL zps_hde( kstp, tb, sb, rhd,  &    ! Partial steps: before horizontal gradient
-            &                              gtu, gsu, gru,  &    ! of t, s, rd at the last ocean level
-            &                              gtv, gsv, grv )
+                         CALL eos( tsb, rhd )                ! before in situ density
+         IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
+            &                                      rhd, gru , grv  )      ! of t, s, rd at the last ocean level
                          CALL ldf_slp( kstp, rhd, rn2b )        ! before slope of the lateral mixing
       ENDIF
@@ -169,8 +168,6 @@
       ! Active tracers                              (ua, va used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                             ta(:,:,:) = 0.e0               ! set tracer trends to zero
-                             sa(:,:,:) = 0.e0
-
-                             CALL tra_swap
+                             tsa(:,:,:,:) = 0.e0               ! set tracer trends to zero
+
                              CALL tra_sbc    ( kstp )       ! surface boundary condition
       IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
@@ -192,21 +189,16 @@
          IF( ln_zdfnpc   )   CALL tra_npc    ( kstp )            ! update after fields by non-penetrative convection
                              CALL tra_nxt    ( kstp )            ! tracer fields at next time step
-                             CALL tra_unswap
-                             CALL eos( ta, sa, rhd, rhop )       ! Time-filtered in situ density for hpg computation
-         IF( ln_zps      )   CALL zps_hde( kstp, ta, sa, rhd,   &   ! Partial steps: time filtered hor. derivative
-            &                                   gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
-            &                                   gtv, gsv, grv )
+                             CALL eos( tsa, rhd, rhop )       ! Time-filtered in situ density for hpg computation
+         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps: time filtered hor. derivative
+            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
          
       ELSE                                                  ! centered hpg  (eos then time stepping)
-                             CALL tra_unswap
-                             CALL eos( tn, sn, rhd, rhop )       ! now in situ density for hpg computation
-         IF( ln_zps      )   CALL zps_hde( kstp, tn, sn, rhd,   &   ! Partial steps: now horizontal derivative
-            &                                   gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
-            &                                   gtv, gsv, grv )
-                             CALL tra_swap
+                             CALL eos( tsn, rhd, rhop )       ! now in situ density for hpg computation
+         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: now hor. derivative
+            &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level
          IF( ln_zdfnpc   )   CALL tra_npc    ( kstp )       ! update after fields by non-penetrative convection
                              CALL tra_nxt    ( kstp )       ! tracer fields at next time step
+      ENDIF 
                              CALL tra_unswap
-      ENDIF 
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/step_oce.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/step_oce.F90	(revision 2082)
@@ -81,4 +81,5 @@
    USE trdmld_rst      ! restart for mixed-layer trends
    USE trdmod_oce      ! ocean momentum/tracers trends
+   USE trdmod          ! momentum/tracers trends   
    USE trdvor          ! vorticity budget                 (trd_vor routine)
    USE diagap          ! hor. mean model-data gap         (dia_gap routine)
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/trc_oce.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/trc_oce.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/trc_oce.F90	(revision 2082)
@@ -20,5 +20,8 @@
    PUBLIC   trc_oce_rgb_read   ! routine called by traqsr.F90
    PUBLIC   trc_oce_ext_lev    ! function called by traqsr.F90 at least
-   
+ 
+   INTEGER , PUBLIC  ::  nn_dttrc      !: frequency of step on passive tracers
+   INTEGER , PUBLIC  ::  nittrc000     !: first time step of passive tracers model
+  
    REAL(wp), PUBLIC , DIMENSION(jpi,jpj,jpk) ::   etot3   !: light absortion coefficient
 
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/C14b/trcsms_c14b.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/C14b/trcsms_c14b.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/C14b/trcsms_c14b.F90	(revision 2082)
@@ -234,7 +234,7 @@
             !  Computation of solubility  
             IF (tmask(ji,jj,1) >  0.) THEN
-               ztp  = ( tn(ji,jj,1) + 273.16 ) * 0.01
+               ztp  = ( tsn(ji,jj,1,jp_tem) + 273.16 ) * 0.01
                zsk  = 0.023517 + ztp * ( -0.023656 + 0.0047036 * ztp )
-               zsol = EXP( -60.2409 + 93.4517 / ztp  + 23.3585 * LOG( ztp ) + zsk * sn(ji,jj,1) )
+               zsol = EXP( -60.2409 + 93.4517 / ztp  + 23.3585 * LOG( ztp ) + zsk * tsn(ji,jj,1,jp_sal) )
                ! convert solubilities [mol/(l * atm)] -> [mol/(m^3 * ppm)]
                zsol = zsol * 1.0e-03
@@ -247,5 +247,5 @@
 
             ! Computes the Schmidt number of CO2 in seawater
-            zt   = tn(ji,jj,1)
+            zt   = tsn(ji,jj,1,jp_tem)
             zsch = 2073.1 + zt * ( -125.62 + zt * (3.6276 - 0.043219 * zt ) )
 
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/CFC/trcsms_cfc.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/CFC/trcsms_cfc.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/CFC/trcsms_cfc.F90	(revision 2082)
@@ -129,8 +129,8 @@
                ! coefficient for solubility for CFC-11/12 in  mol/l/atm
                IF( tmask(ji,jj,1) .GE. 0.5 ) THEN
-                  ztap  = ( tn(ji,jj,1) + 273.16 ) * 0.01
+                  ztap  = ( tsn(ji,jj,1,jp_tem) + 273.16 ) * 0.01
                   zdtap = sob(1,jl) + ztap * ( sob(2,jl) + ztap * sob(3,jl) ) 
                   zsol  =  EXP( soa(1,jl) + soa(2,jl) / ztap + soa(3,jl) * LOG( ztap )   &
-                     &                    + soa(4,jl) * ztap * ztap + sn(ji,jj,1) * zdtap ) 
+                     &                    + soa(4,jl) * ztap * ztap + tsn(ji,jj,1,jp_sal) * zdtap ) 
                ELSE
                   zsol  = 0.e0
@@ -143,5 +143,5 @@
                ! Computation of speed transfert
                !    Schmidt number
-               zt1  = tn(ji,jj,1)
+               zt1  = tsn(ji,jj,1,jp_tem)
                zt2  = zt1 * zt1 
                zt3  = zt1 * zt2
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcbio.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcbio.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcbio.F90	(revision 2082)
@@ -81,5 +81,5 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( kt == nittrc000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' trc_bio: LOBSTER bio-model'
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcexp.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcexp.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcexp.F90	(revision 2082)
@@ -60,5 +60,5 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( kt == nittrc000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' trc_exp: LOBSTER export'
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcopt.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcopt.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcopt.F90	(revision 2082)
@@ -65,5 +65,5 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( kt == nittrc000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' trc_opt : LOBSTER optic-model'
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcsed.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcsed.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/LOBSTER/trcsed.F90	(revision 2082)
@@ -67,5 +67,5 @@
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
+      IF( kt == nittrc000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' trc_sed: LOBSTER sedimentation'
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zche.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zche.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zche.F90	(revision 2082)
@@ -181,8 +181,8 @@
 
             !                             ! SET ABSOLUTE TEMPERATURE
-            ztkel = tn(ji,jj,1) + 273.16
+            ztkel = tsn(ji,jj,1,jp_tem) + 273.16
             zqtt  = ztkel * 0.01
             zqtt2 = zqtt * zqtt
-            zsal  = sn(ji,jj,1) + (1.- tmask(ji,jj,1) ) * 35.
+            zsal  = tsn(ji,jj,1,jp_sal) + (1.- tmask(ji,jj,1) ) * 35.
             zlqtt = LOG( zqtt )
 
@@ -214,7 +214,7 @@
 
                ! SET ABSOLUTE TEMPERATURE
-               ztkel   = tn(ji,jj,jk) + 273.16
+               ztkel   = tsn(ji,jj,jk,jp_tem) + 273.16
                zqtt    = ztkel * 0.01
-               zsal    = sn(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
+               zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35.
                zsqrt  = SQRT( zsal )
                zsal15  = zsqrt * zsal
@@ -224,5 +224,5 @@
                zis2   = zis * zis
                zisqrt = SQRT( zis )
-               ztc     = tn(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
+               ztc     = tsn(ji,jj,jk,jp_tem) + ( 1.- tmask(ji,jj,jk) ) * 20.
 
                ! CHLORINITY (WOOSTER ET AL., 1969)
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zflx.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zflx.F90	(revision 2082)
@@ -128,5 +128,5 @@
 !CDIR NOVERRCHK
          DO ji = 1, jpi
-            ztc  = MIN( 35., tn(ji,jj,1) )
+            ztc  = MIN( 35., tsn(ji,jj,1,jp_tem) )
             ztc2 = ztc * ztc
             ztc3 = ztc * ztc2 
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zint.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zint.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zint.F90	(revision 2082)
@@ -55,6 +55,6 @@
       ! -------------------------------------------
 
-      tgfunc (:,:,:) = EXP( 0.063913 * tn(:,:,:) )
-      tgfunc2(:,:,:) = EXP( 0.07608  * tn(:,:,:) )
+      tgfunc (:,:,:) = EXP( 0.063913 * tsn(:,:,:,jp_tem) )
+      tgfunc2(:,:,:) = EXP( 0.07608  * tsn(:,:,:,jp_tem) )
 
       ! Computation of the silicon dependant half saturation
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zlim.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zlim.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zlim.F90	(revision 2082)
@@ -161,5 +161,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ztemp = MAX( 0., tn(ji,jj,jk) )
+               ztemp = MAX( 0., tsn(ji,jj,jk,jp_tem) )
                xfracal(ji,jj,jk) = caco3r * xlimphy(ji,jj,jk)   &
                   &                       * MAX( 0.0001, ztemp / ( 2.+ ztemp ) )   &
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zprod.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zprod.F90	(revision 2082)
@@ -147,5 +147,5 @@
                ! Computation of the P-I slope for nanos and diatoms
                IF( etot(ji,jj,jk) > 1.E-3 ) THEN
-                   ztn    = MAX( 0., tn(ji,jj,jk) - 15. )
+                   ztn    = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
                    zadap  = 0.+ 1.* ztn / ( 2.+ ztn )
                    zadap2 = 0.e0
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zrem.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zrem.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zrem.F90	(revision 2082)
@@ -287,6 +287,6 @@
                zsatur  = ( sio3eq(ji,jj,jk) - trn(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn )
                zsatur  = MAX( rtrn, zsatur )
-               zsatur2 = zsatur * ( 1. + tn(ji,jj,jk) / 400.)**4
-               znusil  = 0.225  * ( 1. + tn(ji,jj,jk) / 15.) * zsatur + 0.775 * zsatur2**9
+               zsatur2 = zsatur * ( 1. + tsn(ji,jj,jk,jp_tem) / 400.)**4
+               znusil  = 0.225  * ( 1. + tsn(ji,jj,jk,jp_tem) / 15.) * zsatur + 0.775 * zsatur2**9
 #    if defined key_degrad
                zsiremin = xsirem * xstep * znusil * facvol(ji,jj,jk)
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zsed.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zsed.F90	(revision 2082)
@@ -383,8 +383,8 @@
       imois2 = nmonth
 
-      ! 1. first call kt=nit000
+      ! 1. first call kt=nittrc000
       ! -----------------------
 
-      IF( kt == nit000 ) THEN
+      IF( kt == nittrc000 ) THEN
          ! initializations
          nflx1  = 0
@@ -402,5 +402,5 @@
       ! ----------------
 
-      IF( kt == nit000 .OR. imois /= nflx1 ) THEN
+      IF( kt == nittrc000 .OR. imois /= nflx1 ) THEN
 
          ! Calendar computation
@@ -423,13 +423,4 @@
          CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,1), nflx1 )
          CALL iom_get ( numdust, jpdom_data, 'dust', dustmo(:,:,2), nflx2 )
-
-         IF(lwp .AND. nitend-nit000 <= 100 ) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) ' read clio flx ok'
-            WRITE(numout,*)
-               WRITE(numout,*)
-               WRITE(numout,*) 'Clio month: ',nflx1,'  field: dust'
-               CALL prihre( dustmo(:,:,1),jpi,jpj,1,jpi,20,1,jpj,10,1e9,numout )
-         ENDIF
 
       ENDIF
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/par_sed.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/par_sed.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/par_sed.F90	(revision 2082)
@@ -16,4 +16,6 @@
       jpjm1    =>   jpjm1 ,  & !: jpj - 1
       jpij     =>   jpij       !: jpi x jpj
+      jp_tem   =>   jp_tem     !: indice of temperature
+      jp_sal   =>   jp_sal     !: indice of salintity
 
 #if ! defined key_sed_off
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/sed.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/sed.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/sed.F90	(revision 2082)
@@ -34,6 +34,5 @@
 
    USE oce , ONLY :            &
-      tn      =>    tn    ,  & !: pot. temperature (celsius)
-      sn      =>    sn         !: salinity (psu)
+      tsn      =>    tsn       & !: pot. temperature (celsius) and salinity (psu)
 
    USE trc, ONLY :  &
@@ -215,5 +214,4 @@
    INTEGER, PUBLIC :: &
      numsed = 27
-	
 #else
    !!======================================================================
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/seddta.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/seddta.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/SED/seddta.F90	(revision 2082)
@@ -131,12 +131,12 @@
                   trc_data(ji,jj,9 ) = sinking2(ji,jj,ikt)
                   trc_data(ji,jj,10) = sinkcal (ji,jj,ikt)
-                  trc_data(ji,jj,11) = tn      (ji,jj,ikt)
-                  trc_data(ji,jj,12) = sn      (ji,jj,ikt)
+                  trc_data(ji,jj,11) = tsn     (ji,jj,ikt,jp_tem)
+                  trc_data(ji,jj,12) = tsn     (ji,jj,ikt,jp_sal)
 #   else
                   trc_data(ji,jj,7 ) = sinksil (ji,jj,ikt)
                   trc_data(ji,jj,8 ) = sinking (ji,jj,ikt)
                   trc_data(ji,jj,9 ) = sinkcal (ji,jj,ikt)
-                  trc_data(ji,jj,10) = tn      (ji,jj,ikt)
-                  trc_data(ji,jj,11) = sn      (ji,jj,ikt)       
+                  trc_data(ji,jj,10) = tsn     (ji,jj,ikt,jp_tem)
+                  trc_data(ji,jj,11) = tsn     (ji,jj,ikt,jp_sal)       
 #   endif
                ENDIF
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcadv.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcadv.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcadv.F90	(revision 2082)
@@ -33,5 +33,8 @@
 
    PUBLIC   trc_adv    ! routine called by step module
+
    INTEGER ::   nadv   ! choice of the type of advection scheme
+   REAL(wp), DIMENSION(jpk) ::   r2dt  ! vertical profile time-step, = 2 rdttra
+      !                                ! except at nit000 (=rdttra) if neuler=0
 
    !! * Substitutions
@@ -64,4 +67,10 @@
       IF( kt == nittrc000 )   CALL trc_adv_ctl          ! initialisation & control of options
 
+      IF( neuler == 0 .AND. kt == nittrc000 ) THEN     ! at nit000
+         r2dt(:) =  rdttra(:) * FLOAT(nn_dttrc)          ! = rdtra (restarting with Euler time stepping)
+      ELSEIF( kt <= nittrc000 + nn_dttrc ) THEN          ! at nit000 or nit000+1
+         r2dt(:) = 2. * rdttra(:) * FLOAT(nn_dttrc)      ! = 2 rdttra (leapfrog)
+      ENDIF
+
       !                                                   ! effective transport
       DO jk = 1, jpkm1
@@ -78,28 +87,28 @@
       !
       SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==!
-      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )   !  2nd order centered
-      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )   !  TVD 
-      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, 'TRC', zun, zvn, zwn, trb,      tra, jptra )   !  MUSCL 
-      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )   !  MUSCL2 
-      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )   !  UBS 
-      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )   !  QUICKEST 
+      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )   !  2nd order centered
+      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  TVD 
+      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, 'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra )   !  MUSCL 
+      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  MUSCL2 
+      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  UBS 
+      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  QUICKEST 
       !
       CASE (-1 )                                      !==  esopa: test all possibility with control print  ==!
-         CALL tra_adv_cen2  ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_cen2  ( kt, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )          
          WRITE(charout, FMT="('adv1')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
-         CALL tra_adv_tvd   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_tvd   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
          WRITE(charout, FMT="('adv2')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
-         CALL tra_adv_muscl ( kt, 'TRC', zun, zvn, zwn, trb,      tra, jptra )          
+         CALL tra_adv_muscl ( kt, 'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra )          
          WRITE(charout, FMT="('adv3')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
-         CALL tra_adv_muscl2( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_muscl2( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
          WRITE(charout, FMT="('adv4')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
-         CALL tra_adv_ubs   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_ubs   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
          WRITE(charout, FMT="('adv5')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
-         CALL tra_adv_qck   ( kt, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_qck   ( kt, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
          WRITE(charout, FMT="('adv6')")  ; CALL prt_ctl_trc_info(charout)
                                            CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcldf.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 2082)
@@ -70,20 +70,20 @@
 
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt , 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level laplacian
-      CASE ( 1 )   ;   CALL tra_ldf_iso   ( kt , 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )  ! rotated laplacian 
-      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt , 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level bilaplacian
-      CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt , 'TRC',             trb, tra, jptra            )  ! s-coord. horizontal bilaplacian
+      CASE ( 0 )   ;   CALL tra_ldf_lap   ( kt, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level laplacian
+      CASE ( 1 )   ;   CALL tra_ldf_iso   ( kt, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )  ! rotated laplacian 
+      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level bilaplacian
+      CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt, 'TRC',             trb, tra, jptra            )  ! s-coord. horizontal bilaplacian
          !
       CASE ( -1 )                                     ! esopa: test all possibility with control print
-         CALL tra_ldf_lap   ( kt , 'TRC', gtru, gtrv, trb, tra, jptra            )
+         CALL tra_ldf_lap   ( kt, 'TRC', gtru, gtrv, trb, tra, jptra            )
          WRITE(charout, FMT="('ldf0 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-         CALL tra_ldf_iso   ( kt , 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
+         CALL tra_ldf_iso   ( kt, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
          WRITE(charout, FMT="('ldf1 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-         CALL tra_ldf_bilap ( kt , 'TRC', gtru, gtrv, trb, tra, jptra            )
+         CALL tra_ldf_bilap ( kt, 'TRC', gtru, gtrv, trb, tra, jptra            )
          WRITE(charout, FMT="('ldf2 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-         CALL tra_ldf_bilapg( kt , 'TRC',             trb, tra, jptra            )
+         CALL tra_ldf_bilapg( kt, 'TRC',             trb, tra, jptra            )
          WRITE(charout, FMT="('ldf3 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcnxt.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 2082)
@@ -44,5 +44,5 @@
    PUBLIC trc_nxt          ! routine called by step.F90
 
-  REAL(wp), DIMENSION(jpk) ::   r2dt_t
+  REAL(wp), DIMENSION(jpk) ::   r2dt
    !!----------------------------------------------------------------------
    !!   TOP 1.0 , LOCEAN-IPSL (2005) 
@@ -109,6 +109,6 @@
 
       ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt ==  nittrc000) THEN  ;  r2dt_t(:) =     rdttra(:) * FLOAT( nn_dttrc )  ! at nit000             (Euler)
-      ELSEIF( kt <= nittrc000 + 1 )           THEN  ;  r2dt_t(:) = 2.* rdttra(:) * FLOAT( nn_dttrc )  ! at nit000 or nit000+1 (Leapfrog)
+      IF( neuler == 0 .AND. kt ==  nittrc000) THEN  ;  r2dt(:) =     rdttra(:) * FLOAT( nn_dttrc )  ! at nit000             (Euler)
+      ELSEIF( kt <= nittrc000 + 1 )           THEN  ;  r2dt(:) = 2.* rdttra(:) * FLOAT( nn_dttrc )  ! at nit000 or nit000+1 (Leapfrog)
       ENDIF
 
@@ -133,5 +133,5 @@
          DO jn = 1, jptra
             DO jk = 1, jpkm1
-               zfact = 1.e0 / r2dt_t(jk)  
+               zfact = 1.e0 / r2dt(jk)  
                ztrdt(:,:,jk,jn) = ( trb(:,:,jk,jn) - ztrdt(:,:,jk,jn) ) * zfact 
                CALL trd_tra( kt, 'TRC', jn, jptra_trd_atf, ztrdt )
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trctrp.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 2082)
@@ -28,5 +28,5 @@
    USE trcrad          ! positivity                          (trc_rad routine)
    USE trcsbc          ! surface boundary condition          (trc_sbc routine)
-   USE zpshde          ! partial step: hor. derivative   (zps_hde_trc routine)
+   USE zpshde          ! partial step: hor. derivative       (zps_hde routine)
 
 #if defined key_agrif
@@ -74,5 +74,5 @@
 #endif
                                 CALL trc_zdf( kstp )            ! vertical mixing and after tracer fields
-         IF( ln_zps    )        CALL zps_hde_trc( kstp, jptra, trn, gtru, gtrv )  ! Partial steps: now horizontal gradient of passive
+         IF( ln_zps    )        CALL zps_hde( kstp, jptra, trn, gtru, gtrv )  ! Partial steps: now horizontal gradient of passive
                                                                 ! tracers at the bottom ocean level
                                 CALL trc_nxt( kstp )            ! tracer fields at next time step     
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trczdf.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 2082)
@@ -80,12 +80,12 @@
       SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
       CASE ( -1 )                                       ! esopa: test all possibility with control print
-         CALL tra_zdf_exp( kt , 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra ) 
+         CALL tra_zdf_exp( kt, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra ) 
          WRITE(charout, FMT="('zdf1 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-         CALL tra_zdf_imp( kt , 'TRC', r2dt,                trb, tra, jptra ) 
+         CALL tra_zdf_imp( kt, 'TRC', r2dt,                trb, tra, jptra ) 
          WRITE(charout, FMT="('zdf2 ')") ;  CALL prt_ctl_trc_info(charout)
                                             CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
-      CASE ( 0 ) ;  CALL tra_zdf_exp( kt , 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra )    !   explicit scheme 
-      CASE ( 1 ) ;  CALL tra_zdf_imp( kt , 'TRC', r2dt,                trb, tra, jptra )    !   implicit scheme          
+      CASE ( 0 ) ;  CALL tra_zdf_exp( kt, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra )    !   explicit scheme 
+      CASE ( 1 ) ;  CALL tra_zdf_imp( kt, 'TRC', r2dt,                trb, tra, jptra )    !   implicit scheme          
 
       END SELECT
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/oce_trc.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/oce_trc.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/oce_trc.F90	(revision 2082)
@@ -173,9 +173,7 @@
    USE oce , ONLY :   l_traldf_rot => l_traldf_rot  !: rotated laplacian operator for lateral diffusion
 #if defined key_offline
-   USE oce , ONLY :   gtu     =>    gtu     !: t-, s- and rd horizontal gradient at u- and
-   USE oce , ONLY :   gsu     =>    gsu     !: v-points at bottom ocean level
+   USE oce , ONLY :   gtsu    =>    gtsu    !: t-, s- and rd horizontal gradient at u- and
+   USE oce , ONLY :   gtsv    =>    gtsv    !:
    USE oce , ONLY :   gru     =>    gru     !:
-   USE oce , ONLY :   gtv     =>    gtv     !:
-   USE oce , ONLY :   gsv     =>    gsv     !:
    USE oce , ONLY :   grv     =>    grv     !: 
 # if defined key_degrad
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/trc.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/trc.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/trc.F90	(revision 2082)
@@ -55,6 +55,4 @@
    !! passive tracers restart (input and output)
    !! ------------------------------------------  
-   INTEGER , PUBLIC  ::  nn_dttrc      !: frequency of step on passive tracers
-   INTEGER , PUBLIC  ::  nittrc000     !: first time step of passive tracers model
    LOGICAL , PUBLIC  ::  ln_rsttr      !: boolean term for restart i/o for passive tracers (namelist)
    LOGICAL , PUBLIC  ::  lrst_trc      !: logical to control the trc restart write
Index: /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/trcini.F90
===================================================================
--- /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/trcini.F90	(revision 2081)
+++ /branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/trcini.F90	(revision 2082)
@@ -29,5 +29,5 @@
    USE daymod
 #endif
-   USE zpshde          ! partial step: hor. derivative   (zps_hde_trc routine)
+   USE zpshde          ! partial step: hor. derivative   (zps_hde routine)
    USE in_out_manager  ! I/O manager
    USE prtctl_trc      ! Print control passive tracers (prt_ctl_trc_init routine)
@@ -49,7 +49,7 @@
 CONTAINS
    
-   SUBROUTINE trc_ini
+   SUBROUTINE trc_init
       !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE trc_ini  ***
+      !!                     ***  ROUTINE trc_init  ***
       !!
       !! ** Purpose :   Initialization of the passive tracer fields 
@@ -67,5 +67,5 @@
 
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'trc_ini : initial set up of the passive tracers'
+      IF(lwp) WRITE(numout,*) 'trc_init : initial set up of the passive tracers'
       IF(lwp) WRITE(numout,*) '~~~~~~~'
 
@@ -138,5 +138,5 @@
       
       IF( ln_zps .AND. .NOT. lk_trc_c1d )   &              ! Partial steps: before horizontal gradient of passive
-      &                     CALL zps_hde_trc( nittrc000, jptra, trb, gtru, gtrv )       ! tracers at the bottom ocean level
+      &                     CALL zps_hde( nittrc000, jptra, trb, gtru, gtrv )       ! tracers at the bottom ocean level
 
 
@@ -170,5 +170,5 @@
       ENDIF
 
-   END SUBROUTINE trc_ini
+   END SUBROUTINE trc_init
 
 #else
@@ -177,6 +177,6 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_ini                      ! Dummy routine   
-   END SUBROUTINE trc_ini
+   SUBROUTINE trc_init                      ! Dummy routine   
+   END SUBROUTINE trc_init
 #endif
 
