Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 2892)
@@ -93,5 +93,5 @@
       !
       IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   &
-         &              CALL tra_adv_eiv( kt, zun, zvn, zwn, 'TRA' )          ! add the eiv transport (if necessary)
+         &              CALL tra_adv_eiv( kt, nit000, zun, zvn, zwn, 'TRA' )    ! add the eiv transport (if necessary)
       !
       CALL iom_put( "uocetr_eff", zun )                                         ! output effective transport      
@@ -100,28 +100,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', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD 
-      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts )   !  MUSCL 
-      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  MUSCL2 
-      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS 
-      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST 
+      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt,nit000, 'TRA',         zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  2nd order centered
+      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt,nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD 
+      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt,nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts )   !  MUSCL 
+      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt,nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  MUSCL2 
+      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt,nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS 
+      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt,nit000, 'TRA', r2dtra, 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,nit000, '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', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_tvd   ( kt,nit000, 'TRA', r2dtra, 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', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts )          
+         CALL tra_adv_muscl ( kt,nit000, 'TRA', r2dtra, 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', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_muscl2( kt,nit000, 'TRA', r2dtra, 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', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_ubs   ( kt,nit000, 'TRA', r2dtra, 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', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )          
+         CALL tra_adv_qck   ( kt,nit000, 'TRA', r2dtra, 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' )
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90	(revision 2892)
@@ -32,7 +32,4 @@
    USE trc_oce         ! share passive tracers/Ocean variables
    USE lib_mpp         ! MPP library
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
    IMPLICIT NONE
@@ -56,5 +53,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_cen2( kt, cdtype, pun, pvn, pwn,        &
+   SUBROUTINE tra_adv_cen2( kt,nit000, cdtype, pun, pvn, pwn,     &
       &                                 ptb, ptn, pta, kjpt   ) 
       !!----------------------------------------------------------------------
@@ -119,4 +116,5 @@
       !
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
@@ -138,9 +136,5 @@
       ENDIF
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme on ', cdtype
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90	(revision 2892)
@@ -45,5 +45,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_eiv( kt, pun, pvn, pwn, cdtype )
+   SUBROUTINE tra_adv_eiv( kt,nit000, pun, pvn, pwn, cdtype )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_eiv  ***
@@ -69,4 +69,5 @@
 #endif
       INTEGER                         , INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   nit000   ! first time step index
       CHARACTER(len=3)                , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun      ! in : 3 ocean velocity components 
@@ -203,6 +204,7 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE tra_adv_eiv( kt, pun, pvn, pwn, cdtype )              ! Empty routine
+   SUBROUTINE tra_adv_eiv( kt,nit000, pun, pvn, pwn, cdtype )              ! Empty routine
       INTEGER  ::   kt    
+      INTEGER  ::   nit000    
       CHARACTER(len=3) ::   cdtype
       REAL, DIMENSION(:,:,:) ::   pun, pvn, pwn
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 2892)
@@ -25,7 +25,4 @@
    USE diaptr          ! poleward transport diagnostics
    USE trc_oce         ! share passive tracers/Ocean variables
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
 
@@ -47,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_muscl( kt, cdtype, p2dt, pun, pvn, pwn, &
+   SUBROUTINE tra_adv_muscl( kt,nit000,cdtype, p2dt, pun, pvn, pwn, &
       &                                        ptb, pta, kjpt )
       !!----------------------------------------------------------------------
@@ -69,4 +66,5 @@
       !
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
@@ -86,9 +84,5 @@
       ENDIF
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 2892)
@@ -23,7 +23,4 @@
    USE diaptr          ! poleward transport diagnostics
    USE trc_oce         ! share passive tracers/Ocean variables
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
 
@@ -45,5 +42,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_muscl2( kt, cdtype, p2dt, pun, pvn, pwn,      &
+   SUBROUTINE tra_adv_muscl2( kt,nit000, cdtype, p2dt, pun, pvn, pwn,      &
       &                                         ptb, ptn, pta, kjpt )
       !!----------------------------------------------------------------------
@@ -67,4 +64,5 @@
       !!
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
@@ -84,9 +82,5 @@
       ENDIF
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme on ', cdtype
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90	(revision 2892)
@@ -26,7 +26,4 @@
    USE diaptr          ! poleward transport diagnostics
    USE trc_oce         ! share passive tracers/Ocean variables
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
    IMPLICIT NONE
@@ -48,5 +45,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_qck ( kt, cdtype, p2dt, pun, pvn, pwn,      &
+   SUBROUTINE tra_adv_qck ( kt, nit000, cdtype, p2dt, pun, pvn, pwn,      &
       &                                       ptb, ptn, pta, kjpt )
       !!----------------------------------------------------------------------
@@ -85,4 +82,5 @@
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
@@ -93,9 +91,5 @@
       !!----------------------------------------------------------------------
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 2892)
@@ -32,7 +32,4 @@
    USE diaptr          ! poleward transport diagnostics
    USE trc_oce         ! share passive tracers/Ocean variables
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
 
@@ -54,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_tvd ( kt, cdtype, p2dt, pun, pvn, pwn,      &
+   SUBROUTINE tra_adv_tvd ( kt,nit000, cdtype, p2dt, pun, pvn, pwn,      &
       &                                       ptb, ptn, pta, kjpt )
       !!----------------------------------------------------------------------
@@ -74,4 +71,5 @@
       !
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
@@ -92,9 +90,5 @@
       ENDIF
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90	(revision 2892)
@@ -22,7 +22,4 @@
    USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
    USE trc_oce         ! share passive tracers/Ocean variables
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
    IMPLICIT NONE
@@ -43,5 +40,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_ubs ( kt, cdtype, p2dt, pun, pvn, pwn,      &
+   SUBROUTINE tra_adv_ubs ( kt, nit000, cdtype, p2dt, pun, pvn, pwn,      &
       &                                       ptb, ptn, pta, kjpt )
       !!----------------------------------------------------------------------
@@ -83,4 +80,5 @@
       !
       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
@@ -100,9 +98,5 @@
       ENDIF
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 2892)
@@ -34,7 +34,4 @@
    USE lbclnk         ! ocean lateral boundary conditions
    USE prtctl         ! Print control
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
    IMPLICIT NONE
@@ -314,5 +311,5 @@
 
 
-   SUBROUTINE bbl( kt, cdtype )
+   SUBROUTINE bbl( kt, nit000, cdtype )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE bbl  ***
@@ -346,4 +343,5 @@
       !
       INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER         , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
       !!
@@ -392,9 +390,5 @@
       ENDIF
      
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp)  WRITE(numout,*)
          IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 2892)
@@ -68,29 +68,29 @@
 
       SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )   ;   CALL tra_ldf_lap     ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts        )  ! iso-level laplacian
+      CASE ( 0 )   ;   CALL tra_ldf_lap     ( kt,nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts        )  ! iso-level laplacian
       CASE ( 1 )                                                                              ! rotated laplacian
          IF( ln_traldf_grif ) THEN                                                          
-                       CALL tra_ldf_iso_grif( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )      ! Griffies operator
+                       CALL tra_ldf_iso_grif( kt, nit000,'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )      ! Griffies operator
          ELSE                                                                                
-                       CALL tra_ldf_iso     ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )      ! Madec operator
-         ENDIF
-      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. geopot. bilap.
+                       CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )      ! Madec operator
+         ENDIF
+      CASE ( 2 )   ;   CALL tra_ldf_bilap   ( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts        )  ! iso-level bilaplacian
+      CASE ( 3 )   ;   CALL tra_ldf_bilapg  ( kt, nit000, 'TRA',             tsb, tsa, jpts        )  ! s-coord. geopot. bilap.
          !
       CASE ( -1 )                                ! esopa: test all possibility with control print
-         CALL tra_ldf_lap   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts        ) 
+         CALL tra_ldf_lap   ( kt, nit000, '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' )
          IF( ln_traldf_grif ) THEN
-            CALL tra_ldf_iso_grif( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )
+            CALL tra_ldf_iso_grif( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )
          ELSE
-            CALL tra_ldf_iso     ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )  
+            CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )  
          ENDIF
          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', gtsu, gtsv, tsb, tsa, jpts        ) 
+         CALL tra_ldf_bilap ( kt, nit000, '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, nit000, '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/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 2892)
@@ -29,7 +29,4 @@
    USE trc_oce         ! share passive tracers/Ocean variables
    USE lib_mpp         ! MPP library
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
    IMPLICIT NONE
@@ -50,5 +47,5 @@
 CONTAINS
  
-   SUBROUTINE tra_ldf_bilap( kt, cdtype, pgu, pgv,      &
+   SUBROUTINE tra_ldf_bilap( kt, nit000, cdtype, pgu, pgv,      &
       &                                  ptb, pta, kjpt )  
       !!----------------------------------------------------------------------
@@ -82,4 +79,5 @@
       !!
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
@@ -96,9 +94,5 @@
       ENDIF
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator on ', cdtype
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90	(revision 2892)
@@ -25,7 +25,4 @@
    USE trc_oce         ! share passive tracers/Ocean variables
    USE lib_mpp         ! MPP library
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
    IMPLICIT NONE
@@ -45,5 +42,5 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf_bilapg( kt, cdtype, ptb, pta, kjpt )
+   SUBROUTINE tra_ldf_bilapg( kt, nit000, cdtype, ptb, pta, kjpt )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE tra_ldf_bilapg  ***
@@ -73,4 +70,5 @@
       !
       INTEGER         , INTENT(in   )                      ::   kt       ! ocean time-step index
+      INTEGER         , INTENT(in   )                      ::   nit000   ! first time step index
       CHARACTER(len=3), INTENT(in   )                      ::   cdtype   ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                      ::   kjpt     ! number of tracers
@@ -85,9 +83,5 @@
       ENDIF
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate on ', cdtype
@@ -352,5 +346,6 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE tra_ldf_bilapg( kt, cdtype, ptb, pta, kjpt )      ! Empty routine
+   SUBROUTINE tra_ldf_bilapg( kt,nit000, cdtype, ptb, pta, kjpt )      ! Empty routine
+      INTEGER :: kt, nit000
       CHARACTER(len=3) ::   cdtype
       REAL, DIMENSION(:,:,:,:) ::   ptb, pta
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 2892)
@@ -32,7 +32,4 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
 #endif
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
    IMPLICIT NONE
@@ -52,5 +49,5 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf_iso( kt, cdtype, pgu, pgv,              &
+   SUBROUTINE tra_ldf_iso( kt,nit000, cdtype, pgu, pgv,              &
       &                                ptb, pta, kjpt, pahtb0 )
       !!----------------------------------------------------------------------
@@ -99,4 +96,5 @@
       !
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
@@ -124,11 +122,4 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
       ENDIF
-#if defined key_top
-      IF( kt == nittrc000 .AND. cdtype == 'TRC' )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-#endif
       !
       !                                                          ! ===========
@@ -311,5 +302,6 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE tra_ldf_iso( kt, cdtype, pgu, pgv, ptb, pta, kjpt, pahtb0 )      ! Empty routine
+   SUBROUTINE tra_ldf_iso( kt, nit000,cdtype, pgu, pgv, ptb, pta, kjpt, pahtb0 )      ! Empty routine
+      INTEGER:: kt, nit000
       CHARACTER(len=3) ::   cdtype
       REAL, DIMENSION(:,:,:) ::   pgu, pgv   ! tracer gradient at pstep levels
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90	(revision 2892)
@@ -26,7 +26,4 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! MPP library
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
    IMPLICIT NONE
@@ -51,5 +48,5 @@
 CONTAINS
 
-  SUBROUTINE tra_ldf_iso_grif( kt, cdtype, pgu, pgv,              &
+  SUBROUTINE tra_ldf_iso_grif( kt, nit000, cdtype, pgu, pgv,              &
        &                                   ptb, pta, kjpt, pahtb0 )
       !!----------------------------------------------------------------------
@@ -98,4 +95,5 @@
       !
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
@@ -130,9 +128,5 @@
       !zdkt(1:jpi,1:jpj,0:1) => wrk_3d_9(:,:,1:2)
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 2892)
@@ -44,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf_lap( kt, cdtype, pgu, pgv,      &
+   SUBROUTINE tra_ldf_lap( kt, nit000, cdtype, pgu, pgv,      &
       &                                ptb, pta, kjpt ) 
       !!----------------------------------------------------------------------
@@ -66,4 +66,5 @@
       !
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 2892)
@@ -48,7 +48,4 @@
    USE agrif_opa_interp
 #endif
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
    IMPLICIT NONE
@@ -151,6 +148,6 @@
       ELSE                                            ! Leap-Frog + Asselin filter time stepping
          !
-         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
-         ELSE                 ;   CALL tra_nxt_fix( kt, 'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level 
+         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt,nit000, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)     
+         ELSE                 ;   CALL tra_nxt_fix( kt,nit000, 'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level 
          ENDIF
       ENDIF 
@@ -182,5 +179,5 @@
 
 
-   SUBROUTINE tra_nxt_fix( kt, cdtype, ptb, ptn, pta, kjpt )
+   SUBROUTINE tra_nxt_fix( kt, nit000, cdtype, ptb, ptn, pta, kjpt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_nxt_fix  ***
@@ -206,4 +203,5 @@
       !!----------------------------------------------------------------------
       INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
+      INTEGER         , INTENT(in   )                               ::   nit000   ! first time step index
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
@@ -217,9 +215,5 @@
       !!----------------------------------------------------------------------
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping', cdtype
@@ -252,5 +246,5 @@
 
 
-   SUBROUTINE tra_nxt_vvl( kt, cdtype, ptb, ptn, pta, kjpt )
+   SUBROUTINE tra_nxt_vvl( kt, nit000, cdtype, ptb, ptn, pta, kjpt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_nxt_vvl  ***
@@ -277,4 +271,5 @@
       !!----------------------------------------------------------------------
       INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
+      INTEGER         , INTENT(in   )                               ::   nit000   ! first time step index
       CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
       INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
@@ -289,9 +284,5 @@
       !!----------------------------------------------------------------------
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping', cdtype
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 2892)
@@ -76,11 +76,11 @@
 
       SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )    ;    CALL tra_zdf_exp( kt, 'TRA', r2dtra, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme 
-      CASE ( 1 )    ;    CALL tra_zdf_imp( kt, 'TRA', r2dtra,            tsb, tsa, jpts )  !   implicit scheme 
+      CASE ( 0 )    ;    CALL tra_zdf_exp( kt, nit000, 'TRA', r2dtra, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme 
+      CASE ( 1 )    ;    CALL tra_zdf_imp( kt, nit000, 'TRA', r2dtra,            tsb, tsa, jpts )  !   implicit scheme 
       CASE ( -1 )                                       ! esopa: test all possibility with control print
-         CALL tra_zdf_exp( kt, 'TRA', r2dtra, nn_zdfexp, tsb, tsa, jpts )
+         CALL tra_zdf_exp( kt, nit000, 'TRA', r2dtra, 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', r2dtra,            tsb, tsa, jpts ) 
+         CALL tra_zdf_imp( kt, nit000, 'TRA', r2dtra,            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/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_exp.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_exp.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_exp.F90	(revision 2892)
@@ -31,7 +31,4 @@
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
    IMPLICIT NONE
@@ -51,5 +48,5 @@
 CONTAINS
 
-   SUBROUTINE tra_zdf_exp( kt, cdtype, p2dt, kn_zdfexp,   &
+   SUBROUTINE tra_zdf_exp( kt,nit000, cdtype, p2dt, kn_zdfexp,   &
       &                                ptb , pta      , kjpt )
       !!----------------------------------------------------------------------
@@ -80,4 +77,5 @@
       !
       INTEGER                              , INTENT(in   ) ::   kt          ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000      ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype      ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt        ! number of tracers
@@ -96,9 +94,5 @@
       ENDIF
 
-#if defined key_top
-      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
-#else
       IF( kt == nit000 )  THEN
-#endif
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_zdf_exp : explicit vertical mixing on ', cdtype
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 2892)
@@ -35,7 +35,4 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! MPP library
-#if defined key_top
-   USE trc, ONLY: nittrc000  !get first time step for passive tracers
-#endif
 
    IMPLICIT NONE
@@ -58,5 +55,5 @@
 CONTAINS
  
-   SUBROUTINE tra_zdf_imp( kt, cdtype, p2dt, ptb, pta, kjpt ) 
+   SUBROUTINE tra_zdf_imp( kt,nit000, cdtype, p2dt, ptb, pta, kjpt ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_zdf_imp  ***
@@ -82,4 +79,5 @@
       !
       INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   nit000          ! first time step index
       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
       INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
@@ -105,15 +103,4 @@
          ENDIF
       ENDIF
-#if defined key_top
-      IF( kt == nittrc000 .AND. cdtype == 'TRC' )  THEN
-         IF(lwp)WRITE(numout,*)
-         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
-         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
-         !
-         IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator
-         ELSE                ;    r_vvl = 0._wp       
-         ENDIF
-      ENDIF
-#endif
       !
       !                                               ! ============= !
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90	(revision 2892)
@@ -103,31 +103,31 @@
 
       !                                                   ! add the eiv transport (if necessary)
-      IF( lk_traldf_eiv )   CALL tra_adv_eiv( kt, zun, zvn, zwn, 'TRC' )
+      IF( lk_traldf_eiv )   CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' )
       !
       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', 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 )   ;    CALL tra_adv_cen2  ( kt, nittrc000, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )   !  2nd order centered
+      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  TVD 
+      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra )   !  MUSCL 
+      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  MUSCL2 
+      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  UBS 
+      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, nittrc000, '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, nittrc000, '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', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_tvd   ( kt, nittrc000, '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', r2dt, zun, zvn, zwn, trb,      tra, jptra )          
+         CALL tra_adv_muscl ( kt, nittrc000, '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', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_muscl2( kt, nittrc000, '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', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_ubs   ( kt, nittrc000, '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', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
+         CALL tra_adv_qck   ( kt, nittrc000, '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/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90	(revision 2892)
@@ -67,20 +67,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, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level laplacian
+      CASE ( 1 )   ;   CALL tra_ldf_iso   ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )  ! rotated laplacian 
+      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level bilaplacian
+      CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt, nittrc000, '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, nittrc000, '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, nittrc000, '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, nittrc000, '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, nittrc000, '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/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90	(revision 2892)
@@ -139,6 +139,6 @@
       ELSE
          ! Leap-Frog + Asselin filter time stepping
-         IF( lk_vvl ) THEN   ;   CALL tra_nxt_vvl( kt, 'TRC', trb, trn, tra, jptra )      ! variable volume level (vvl) 
-         ELSE                ;   CALL tra_nxt_fix( kt, 'TRC', trb, trn, tra, jptra )      ! fixed    volume level 
+         IF( lk_vvl ) THEN   ;   CALL tra_nxt_vvl( kt,nittrc000, 'TRC', trb, trn, tra, jptra )      ! variable volume level (vvl) 
+         ELSE                ;   CALL tra_nxt_fix( kt,nittrc000, 'TRC', trb, trn, tra, jptra )      ! fixed    volume level 
          ENDIF
       ENDIF
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 2892)
@@ -88,12 +88,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,nittrc000, '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,nittrc000, '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,nittrc000, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra )    !   explicit scheme 
+      CASE ( 1 ) ;  CALL tra_zdf_imp( kt,nittrc000, 'TRC', r2dt,                trb, tra, jptra )    !   implicit scheme          
 
       END SELECT
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 2892)
@@ -108,4 +108,5 @@
    USE dom_oce , ONLY :   e3w_0      =>   e3w_0      !: reference depth of w-points (m)
    USE dom_oce , ONLY :   gdepw_0    =>   gdepw_0    !: reference depth of w-points (m)
+# if ! defined key_zco
    USE dom_oce , ONLY :   gdep3w     =>  gdep3w      !: ???
    USE dom_oce , ONLY :   gdept      =>  gdept       !: depth of t-points (m)
@@ -118,5 +119,5 @@
    USE dom_oce , ONLY :   e3uw       =>  e3uw        !: uw-points (m)
    USE dom_oce , ONLY :   e3vw       =>  e3vw        !: vw-points (m)
-
+# endif
    USE dom_oce , ONLY :   ln_zps     =>  ln_zps      !: partial steps flag
    USE dom_oce , ONLY :   ln_sco     =>  ln_sco      !: s-coordinate flag
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/sshwzv_trc.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/sshwzv_trc.F90	(revision 2891)
+++ 	(revision )
@@ -1,193 +1,0 @@
-MODULE sshwzv_trc   
-   !!==============================================================================
-   !!                       ***  MODULE  sshwzv  ***
-   !! Ocean dynamics : sea surface height and vertical velocity
-   !!==============================================================================
-   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
-   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA 
-   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
-   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   ssh_wzv        : after ssh & now vertical velocity
-   !!   ssh_nxt        : filter ans swap the ssh arrays
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE domvvl          ! Variable volume
-   USE divcur          ! hor. divergence and curl      (div & cur routines)
-   USE iom             ! I/O library
-   USE restart         ! only for lrst_oce
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE phycst
-   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
-   USE lib_mpp         ! MPP library
-   USE obc_par         ! open boundary cond. parameter
-   USE obc_oce
-   USE bdy_oce
-   USE diaar5, ONLY:   lk_diaar5
-   USE iom
-   USE sbcrnf, ONLY: h_rnf, nk_rnf   ! River runoff 
-   USE trc, ONLY:  nn_dttrc, nittrc000 !TOP_TRC substepping
-#if defined key_agrif
-   USE agrif_opa_update
-   USE agrif_opa_interp
-#endif
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ssh_wzv_trc    ! called by step.F90
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id: sshwzv.F90 2715 2011-03-30 15:58:35Z rblod $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE ssh_wzv_trc( kt ) 
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE ssh_wzv  ***
-      !!                   
-      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
-      !!              and update the now vertical coordinate (lk_vvl=T).
-      !!
-      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
-      !!      velocity is computed by integrating the horizontal divergence  
-      !!      from the bottom to the surface minus the scale factor evolution.
-      !!        The boundary conditions are w=0 at the bottom (no flux) and.
-      !!
-      !! ** action  :   ssha    : after sea surface height
-      !!                wn      : now vertical velocity
-      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
-      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
-      !!
-      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
-      !!----------------------------------------------------------------------
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE oce     , ONLY:   z3d   => ta                           ! ta used as 3D workspace
-      USE wrk_nemo, ONLY:   zhdiv => wrk_2d_1 , z2d => wrk_2d_2   ! 2D workspace
-      !
-      INTEGER, INTENT(in) ::   kt   ! time step
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
-      !!----------------------------------------------------------------------
-
-      IF( wrk_in_use(2, 1,2) ) THEN
-         CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable')   ;   RETURN
-      ENDIF
-
-      IF( kt == nittrc000 ) THEN
-         !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'ssh_wzv_trc : after sea surface height and now vertical velocity '
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
-         !
-         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
-         !
-         !
-      ENDIF
-
-      !                                           !------------------------------------------!
-      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
-         !                                        !------------------------------------------!
-         DO jk = 1, jpkm1
-            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
-            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
-            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
-            !
-            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
-            fse3u (:,:,jk) = fse3u_n (:,:,jk)
-            fse3v (:,:,jk) = fse3v_n (:,:,jk)
-            fse3f (:,:,jk) = fse3f_n (:,:,jk)
-            fse3w (:,:,jk) = fse3w_n (:,:,jk)
-            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
-            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
-         END DO
-         !
-         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
-         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
-         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
-         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
-         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )
-         ! 
-      ENDIF
-      !
-      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
-      !
-      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nittrc000 )   z2dt = rdt
-
-      !                                           !------------------------------!
-      !                                           !   After Sea Surface Height   !
-      !                                           !------------------------------!
-      zhdiv(:,:) = 0._wp
-      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
-        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
-      END DO
-      !                                                ! Sea surface elevation time stepping
-      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
-      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
-      z1_rau0 = 0.5 / rau0
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
-
-#if defined key_agrif
-      CALL agrif_ssh( kt )
-#endif
-#if defined key_obc
-      IF( Agrif_Root() ) THEN 
-         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
-         CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm)
-      ENDIF
-#endif
-#if defined key_bdy
-      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
-      CALL lbc_lnk( ssha, 'T', 1. ) 
-#endif
-
-      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
-      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1      ! NO Vector Opt.
-               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
-                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
-                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
-               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
-                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
-                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
-            END DO
-         END DO
-         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
-      ENDIF
-      
-
-      !                                           !------------------------------!
-      !                                           !     Now Vertical Velocity    !
-      !                                           !------------------------------!
-      z1_2dt = 1.e0 / z2dt
-      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
-         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
-         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
-            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
-            &                         * tmask(:,:,jk) * z1_2dt
-#if defined key_bdy
-         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
-#endif
-      END DO
-
-      !
-      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('ssh_wzv: failed to release workspace arrays')
-      !
-   END SUBROUTINE ssh_wzv_trc
-
-
-   !!======================================================================
-END MODULE sshwzv_trc
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 2892)
@@ -119,4 +119,37 @@
   REAL(wp), PUBLIC,ALLOCATABLE, SAVE, DIMENSION (:,:,:) ::   vslp_tm !: j-direction slope at u-, w-points
 #endif
+
+  REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:,:,:) :: tsn_tm      !: t/s average     [m/s]
+  REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) ::     fr_i_tm    !: average ice fraction     [m/s]
+  REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) ::     hmld_tm    !: mixed layer depth average [m]
+  REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) ::     qsr_tm     !: solar radiation average [m]
+  REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) ::     wndm_tm    !: 10m wind average [m]
+  REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) ::     emp_tm     !: freshwater budget: volume flux [Kg/m2/s]
+  REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) ::     emp_b_hold !:hold emp from the beginning of each sub-stepping[m]  
+  REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) ::    emps_tm    !: freshwater budget:concentration/dilution [Kg/m2/s]
+# if defined key_zdfddm
+  REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:,:) ::  avs_tm      !: salinity vertical diffusivity coeff. at w-point   [m/s]
+# endif
+#if defined key_traldf_c3d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahtt_tm, ahtu_tm, ahtv_tm, ahtw_tm   !: ** 3D coefficients ** at T-,U-,V-,W-points
+#elif defined key_traldf_c2d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ahtt_tm, ahtu_tm, ahtv_tm, ahtw_tm   !: ** 2D coefficients ** at T-,U-,V-,W-points
+#elif defined key_traldf_c1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::   ahtt_tm, ahtu_tm, ahtv_tm, ahtw_tm   !: ** 1D coefficients ** at T-,U-,V-,W-points
+#else
+   REAL(wp), PUBLIC                                      ::   ahtt_tm, ahtu_tm, ahtv_tm, ahtw_tm   !: ** 0D coefficients ** at T-,U-,V-,W-points
+#endif
+#if defined key_traldf_eiv
+# if defined key_traldf_c3d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aeiu_tm , aeiv_tm , aeiw_tm   !: ** 3D coefficients **
+# elif defined key_traldf_c2d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aeiu_tm , aeiv_tm , aeiw_tm   !: ** 2D coefficients **
+# elif defined key_traldf_c1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::   aeiu_tm , aeiv_tm, aeiw_tm   !: ** 1D coefficients **
+# else
+   REAL(wp), PUBLIC                                      ::   aeiu_tm , aeiv_tm , aeiw_tm   !: ** 0D coefficients **
+# endif
+# endif
+
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)  ::   rnf_tm !: river runoff
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)  ::   h_rnf_tm !: depth in metres to the bottom of the relevant grid box
@@ -133,4 +166,29 @@
   REAL(wp), PUBLIC,ALLOCATABLE, SAVE, DIMENSION(:,:,:)::     hdivn_temp, rotn_temp
   REAL(wp), PUBLIC,ALLOCATABLE, SAVE, DIMENSION(:,:,:)::     hdivb_temp, rotb_temp
+  REAL(wp), PUBLIC,ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  tsn_temp
+  REAL(wp), PUBLIC,ALLOCATABLE, SAVE, DIMENSION(:,:)::       hmld_temp, qsr_temp, emp_temp, emps_temp,fr_i_temp,wndm_temp,emp_b_temp
+# if defined key_zdfddm
+  REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:,:) ::  avs_temp      !: salinity vertical diffusivity coeff. at w-point   [m/s]
+# endif
+#if defined key_traldf_c3d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  ahtt_temp, ahtu_temp, ahtv_temp, ahtw_temp   !: ** 3D coefficients ** at T-,U-,V-,W-points
+#elif defined key_traldf_c2d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  ahtt_temp, ahtu_temp, ahtv_temp, ahtw_temp   !: ** 2D coefficients ** at T-,U-,V-,W-points
+#elif defined key_traldf_c1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::  ahtt_temp, ahtu_temp, ahtv_temp, ahtw_temp   !: ** 1D coefficients ** at T-,U-,V-,W-points
+#else
+   REAL(wp), PUBLIC                                      ::  ahtt_temp, ahtu_temp, ahtv_temp, ahtw_temp   !: ** 0D coefficients ** at T-,U-,V-,W-points
+#endif
+#if defined key_traldf_eiv
+# if defined key_traldf_c3d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  aeiu_temp , aeiv_temp , aeiw_temp   !: ** 3D coefficients **
+# elif defined key_traldf_c2d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  aeiu_temp , aeiv_temp , aeiw_temp   !: ** 2D coefficients **
+# elif defined key_traldf_c1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::  aeiu_temp , aeiv_temp, aeiw_temp   !: ** 1D coefficients **
+# else
+   REAL(wp), PUBLIC                                      ::  aeiu_temp , aeiv_temp , aeiw_temp   !: ** 0D coefficients **
+# endif
+# endif
 
    !!----------------------------------------------------------------------
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 2892)
@@ -29,5 +29,5 @@
    USE zpshde          ! partial step: hor. derivative   (zps_hde routine)
    USE prtctl_trc      ! Print control passive tracers (prt_ctl_trc_init routine)
-   USE trcsubstp       ! variables to substep passive tracers
+   USE trcsub       ! variables to substep passive tracers
    
    IMPLICIT NONE
@@ -174,5 +174,5 @@
 
       IF(nn_dttrc .NE. 1) THEN
-          CALL trc_sub_stp_ini      !initialize variables for substepping passive tracers
+          CALL trc_sub_ini      !initialize variables for substepping passive tracers
       ENDIF
       !
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 2892)
@@ -77,5 +77,5 @@
       ! we open and define the tracer restart file one tracer time step before writing the data (-> at nitrst - 2*nn_dttrc + 1)
       ! except if we write tracer restart files every tracer time step or if a tracer restart file was writen at nitend - 2*nn_dttrc + 1
-      IF( kt == nitrst - 2*nn_dttrc + 1 .OR. nstock == nn_dttrc .OR. ( kt == nitend - nn_dttrc + 1 .AND. .NOT. lrst_trc ) ) THEN
+      IF( kt == nitrst - 2*nn_dttrc .OR. nstock == nn_dttrc .OR. ( kt == nitend - nn_dttrc .AND. .NOT. lrst_trc ) ) THEN
          ! beware of the format used to write kt (default is i8.8, that should be large enough)
          IF( nitrst > 1.0e9 ) THEN   ;   WRITE(clkt,*       ) nitrst
@@ -151,5 +151,4 @@
       REAL(wp) :: zarak0
       !!----------------------------------------------------------------------
-
 
       CALL trc_rst_cal( kt, 'WRITE' )   ! calendar
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 2891)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 2892)
@@ -22,5 +22,5 @@
    USE iom
    USE in_out_manager
-   USE trcsubstp
+   USE trcsub
 
    IMPLICIT NONE
@@ -77,5 +77,5 @@
          !
          IF(nn_dttrc .NE. 1) THEN
-            CALL trc_sub_stp_reset(kt)
+            CALL trc_sub_reset(kt)
          ENDIF
       ENDIF
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcsub.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 2892)
+++ branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 2892)
@@ -0,0 +1,1151 @@
+MODULE trcsub
+   !!======================================================================
+   !!                       ***  MODULE trcsubstp  ***
+   !!TOP :   Averages physics variables for TOP substepping. 
+   !!======================================================================
+   !! History :  1.0  !  2011-10  (K. Edwards)  Original
+   !!----------------------------------------------------------------------
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   trc_sub    : passive tracer system sub-stepping 
+   !!----------------------------------------------------------------------
+   USE oce_trc          ! ocean dynamics and active tracers variables
+   USE trc
+   USE prtctl_trc       ! Print control for debbuging
+   USE iom
+   USE in_out_manager
+   USE lbclnk
+#if defined key_zdftke
+   USE zdftke          ! twice TKE (en)
+#endif
+#if defined key_zdfgls
+   USE zdfgls, ONLY: en
+#endif
+   USE zdf_oce
+   USE domvvl
+   USE divcur          ! hor. divergence and curl      (div & cur routines)
+   USE sbcrnf, ONLY: h_rnf, nk_rnf   ! River runoff 
+   USE trc, ONLY:  nn_dttrc, nittrc000 !TOP_TRC substepping
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE bdy_oce
+#if defined key_agrif
+   USE agrif_opa_update
+   USE agrif_opa_interp
+#endif
+
+   IMPLICIT NONE
+
+   PUBLIC   trc_sub_stp    ! called by trc_stp
+   PUBLIC   trc_sub_ini    ! called by trc_ini to initialize substepping arrays.
+   PUBLIC   trc_sub_reset  ! called by trc_stp to reset physics variables
+   PUBLIC   trc_sub_ssh    ! called by trc_stp to reset physics variables
+
+#  include "top_substitute.h90"
+   
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id: trcstp.F90 2528 2010-12-27 17:33:53Z rblod $ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE trc_sub_stp( kt )
+      !!-------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_stp  ***
+      !!                      
+      !! ** Purpose : Average variables needed for sub-stepping passive tracers
+      !! 
+      !! ** Method  : Called every timestep to increment _tm (time mean) variables
+      !!              on TOP steps, calculate averages.
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::  kt        !ocean time-step index
+      INTEGER ::	        ji,jj,jk  !dummy loop indices
+      !!-------------------------------------------------------------------
+
+       IF (kt == nit000) THEN
+           IF(lwp) WRITE(numout,*)
+           IF(lwp) WRITE(numout,*) 'trc_sub_stp : substepping of the passive tracers'
+           IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+
+           sshb_hold(:,:)=sshn(:,:)
+           sshu_b_hold(:,:)=sshu_n(:,:)
+           sshv_b_hold(:,:)=sshv_n(:,:)
+           emp_b_hold(:,:)=emp_b(:,:)
+       ENDIF  
+       IF( MOD( kt , nn_dttrc ) /= 0 ) THEN
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  un_tm(ji,jj,jk)= un_tm(ji,jj,jk)+un(ji,jj,jk)*fse3u(ji,jj,jk)
+                  vn_tm(ji,jj,jk)= vn_tm(ji,jj,jk)+vn(ji,jj,jk)*fse3v(ji,jj,jk)
+		  wn_tm(ji,jj,jk)= wn_tm(ji,jj,jk)+wn(ji,jj,jk)*fse3w(ji,jj,jk)
+                  avt_tm(ji,jj,jk)= avt_tm(ji,jj,jk)+avt(ji,jj,jk)*fse3t(ji,jj,jk)
+                  tsn_tm(ji,jj,jk,jp_sal)= tsn_tm(ji,jj,jk,jp_sal)+tsn(ji,jj,jk,jp_sal)*fse3t(ji,jj,jk)
+                  tsn_tm(ji,jj,jk,jp_tem)= tsn_tm(ji,jj,jk,jp_tem)+tsn(ji,jj,jk,jp_tem)*fse3t(ji,jj,jk)
+# if defined key_zdfddm
+                  avs_tm(ji,jj,jk)= avs_tm(ji,jj,jk)+avs(ji,jj,jk)*fse3w(ji,jj,jk)
+# endif
+#if defined key_ldfslp
+		  wslpi_tm(ji,jj,jk)= wslpi_tm(ji,jj,jk)+wslpi(ji,jj,jk)*fse3w(ji,jj,jk)
+		  wslpj_tm(ji,jj,jk)= wslpj_tm(ji,jj,jk)+wslpj(ji,jj,jk)*fse3w(ji,jj,jk)
+		  uslp_tm(ji,jj,jk)= uslp_tm(ji,jj,jk)+uslp(ji,jj,jk)*fse3u(ji,jj,jk)
+		  vslp_tm(ji,jj,jk)= vslp_tm(ji,jj,jk)+vslp(ji,jj,jk)*fse3v(ji,jj,jk)
+#endif
+               END DO
+            END DO
+         END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+         	sshn_tm(ji,jj)=sshn_tm(ji,jj)+sshn(ji,jj)
+                sshu_n_tm(ji,jj)=sshu_n_tm(ji,jj)+sshu_n(ji,jj)
+                sshv_n_tm(ji,jj)=sshv_n_tm(ji,jj)+sshv_n(ji,jj)
+	        rnf_tm(ji,jj)= rnf_tm(ji,jj)+rnf(ji,jj)
+	        h_rnf_tm(ji,jj)= h_rnf_tm(ji,jj)+h_rnf(ji,jj)
+         	hmld_tm(ji,jj)=hmld_tm(ji,jj)+hmld(ji,jj)
+         	qsr_tm(ji,jj)=qsr_tm(ji,jj)+qsr(ji,jj)
+         	fr_i_tm(ji,jj)=fr_i_tm(ji,jj)+fr_i(ji,jj)
+         	wndm_tm(ji,jj)=wndm_tm(ji,jj)+wndm(ji,jj)
+         	emp_tm(ji,jj)=emp_tm(ji,jj)+emp(ji,jj)
+         	emps_tm(ji,jj)=emps_tm(ji,jj)+emps(ji,jj)
+             END DO
+         END DO
+#if defined key_traldf_c3d
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ahtt_tm(ji,jj,jk)= ahtt_tm(ji,jj,jk)+ahtt(ji,jj,jk)*fse3t(ji,jj,jk)
+                  ahtu_tm(ji,jj,jk)= ahtu_tm(ji,jj,jk)+ahtu(ji,jj,jk)*fse3u(ji,jj,jk)
+                  ahtv_tm(ji,jj,jk)= ahtv_tm(ji,jj,jk)+ahtv(ji,jj,jk)*fse3v(ji,jj,jk)
+                  ahtw_tm(ji,jj,jk)= ahtw_tm(ji,jj,jk)+ahtw(ji,jj,jk)*fse3w(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+#elif defined key_traldf_c2d
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ahtt_tm(ji,jj)= ahtt_tm(ji,jj)+ahtt(ji,jj)
+                  ahtu_tm(ji,jj)= ahtu_tm(ji,jj)+ahtu(ji,jj)
+                  ahtv_tm(ji,jj)= ahtv_tm(ji,jj)+ahtv(ji,jj)
+                  ahtw_tm(ji,jj)= ahtw_tm(ji,jj)+ahtw(ji,jj)
+               END DO
+            END DO
+#elif defined key_traldf_c1d
+          DO jk = 1, jpk
+             ahtt_tm(jk)= ahtt_tm(jk)+ahtt(jk)
+             ahtu_tm(jk)= ahtu_tm(jk)+ahtu(jk)
+             ahtv_tm(jk)= ahtv_tm(jk)+ahtv(jk)
+             ahtw_tm(jk)= ahtw_tm(jk)+ahtw(jk)
+          END DO
+#else
+             ahtt_tm= ahtt_tm+ahtt
+             ahtu_tm= ahtu_tm+ahtu
+             ahtv_tm= ahtv_tm+ahtv
+             ahtw_tm= ahtw_tm+ahtw
+#endif
+#if defined key_traldf_eiv
+# if defined key_traldf_c3d
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  aeiu_tm(ji,jj,jk)= aeiu_tm(ji,jj,jk)+aeiu(ji,jj,jk)*fse3u(ji,jj,jk)
+                  aeiv_tm(ji,jj,jk)= aeiv_tm(ji,jj,jk)+aeiv(ji,jj,jk)*fse3v(ji,jj,jk)
+                  aeiw_tm(ji,jj,jk)= aeiw_tm(ji,jj,jk)+aeiw(ji,jj,jk)*fse3w(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+# elif defined key_traldf_c2d
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  aeiu_tm(ji,jj)= aeiu_tm(ji,jj)+aeiu(ji,jj)
+                  aeiv_tm(ji,jj)= aeiv_tm(ji,jj)+aeiv(ji,jj)
+                  aeiw_tm(ji,jj)= aeiw_tm(ji,jj)+aeiw(ji,jj)
+               END DO
+            END DO
+# elif defined key_traldf_c1d
+          DO jk = 1, jpk
+             aeiu_tm(jk)= aeiu_tm(jk)+aeiu(jk)
+             aeiv_tm(jk)= aeiv_tm(jk)+aeiv(jk)
+             aeiw_tm(jk)= aeiw_tm(jk)+aeiw(jk)
+          END DO
+# else
+             aeiu_tm= aeiu_tm+aeiu
+             aeiv_tm= aeiv_tm+aeiv
+             aeiw_tm= aeiw_tm+aeiw
+# endif
+# endif
+
+      ELSE
+!   it is time to substep
+!   1. set temporary arrays to hold physics variables
+            sshn_temp(:,:)=sshn(:,:)
+            sshu_n_temp(:,:)=sshu_n(:,:)
+            sshv_n_temp(:,:)=sshv_n(:,:)
+            sshf_n_temp(:,:)=sshf_n(:,:)
+            sshb_temp(:,:)=sshb(:,:)
+            sshu_b_temp(:,:)=sshu_b(:,:)
+            sshv_b_temp(:,:)=sshv_b(:,:)
+            ssha_temp(:,:)=ssha(:,:)
+            sshu_a_temp(:,:)=sshu_a(:,:)
+            sshv_a_temp(:,:)=sshv_a(:,:)
+            rnf_temp(:,:)=rnf(:,:)
+            h_rnf_temp(:,:)=h_rnf(:,:)
+            rotn_temp(:,:,:)=rotn(:,:,:)
+            hdivn_temp(:,:,:)=hdivn(:,:,:)
+            rotb_temp(:,:,:)=rotb(:,:,:)
+            hdivb_temp(:,:,:)=hdivb(:,:,:)
+            un_temp(:,:,:)=un(:,:,:)
+            vn_temp(:,:,:)=vn(:,:,:)
+            wn_temp(:,:,:)=wn(:,:,:)
+            tsn_temp(:,:,:,:)=tsn(:,:,:,:)
+            avt_temp(:,:,:)=avt(:,:,:)
+            hmld_temp(:,:)=hmld(:,:)
+            qsr_temp(:,:)=qsr(:,:)
+            wndm_temp(:,:)=wndm(:,:)
+            emp_temp(:,:)=emp(:,:)
+            emp_b_temp(:,:)=emp_b(:,:)
+            emps_temp(:,:)=emps(:,:)
+#if defined key_traldf_c3d
+            ahtu_temp(:,:,:)=ahtu(:,:,:)
+            ahtv_temp(:,:,:)=ahtv(:,:,:)
+            ahtw_temp(:,:,:)=ahtw(:,:,:)
+            ahtt_temp(:,:,:)=ahtt(:,:,:)
+#elif defined key_traldf_c2d
+            ahtu_temp(:,:)=ahtu(:,:)
+            ahtv_temp(:,:)=ahtv(:,:)
+            ahtw_temp(:,:)=ahtw(:,:)
+            ahtt_temp(:,:)=ahtt(:,:)
+#elif defined key_traldf_c1d
+            ahtu_temp(:)=ahtu(:)
+            ahtv_temp(:)=ahtv(:)
+            ahtw_temp(:)=ahtw(:)
+            ahtt_temp(:)=ahtt(:)
+#else
+            ahtu_temp=ahtu
+            ahtv_temp=ahtv
+            ahtw_temp=ahtw
+            ahtt_temp=ahtt
+#endif
+# if defined key_zdfddm
+            avs_temp(:,:,:)=avs(:,:,:)
+# endif
+#if defined key_traldf_eiv
+#if defined key_traldf_c3d
+            aeiu_temp(:,:,:)=aeiu(:,:,:)
+            aeiv_temp(:,:,:)=aeiv(:,:,:)
+            aeiw_temp(:,:,:)=aeiw(:,:,:)
+#elif defined key_traldf_c2d
+            aeiu_temp(:,:)=aeiu(:,:)
+            aeiv_temp(:,:)=aeiv(:,:)
+            aeiw_temp(:,:)=aeiw(:,:)
+#elif defined key_traldf_c1d
+            aeiu_temp(:)=aeiu(:)
+            aeiv_temp(:)=aeiv(:)
+            aeiw_temp(:)=aeiw(:)
+#else
+            aeiu_temp=aeiu
+            aeiv_temp=aeiv
+            aeiw_temp=aeiw
+#endif
+#endif
+#if defined key_ldfslp
+            wslpi_temp(:,:,:)=wslpi(:,:,:)
+            wslpj_temp(:,:,:)=wslpj(:,:,:)
+            uslp_temp(:,:,:)=uslp(:,:,:)
+            vslp_temp(:,:,:)=vslp(:,:,:)
+#endif
+            hu_temp(:,:)=hu(:,:)
+            hv_temp(:,:)=hv(:,:)
+            hur_temp(:,:)=hur(:,:)
+            hvr_temp(:,:)=hvr(:,:)
+
+!   2. create averages and reassign variables
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  sshn_tm(ji,jj)= sshn_tm(ji,jj)+sshn(ji,jj)
+                  sshn(ji,jj)=sshn_tm(ji,jj)/(nn_dttrc+1)
+		  
+                  sshu_n_tm(ji,jj)= sshu_n_tm(ji,jj)+sshu_n(ji,jj)
+                  sshu_n(ji,jj)=sshu_n_tm(ji,jj)/(nn_dttrc+1)
+		  
+                  sshv_n_tm(ji,jj)= sshv_n_tm(ji,jj)+sshv_n(ji,jj)
+                  sshv_n(ji,jj)=sshv_n_tm(ji,jj)/(nn_dttrc+1)
+
+		  sshb(ji,jj)=sshb_hold(ji,jj)
+
+		  sshu_b(ji,jj)=sshu_b_hold(ji,jj)
+
+		  sshv_b(ji,jj)=sshv_b_hold(ji,jj)
+
+                  rnf_tm(ji,jj)= rnf_tm(ji,jj)+rnf(ji,jj)
+                  rnf(ji,jj)=rnf_tm(ji,jj)/(nn_dttrc+1)
+ 
+                  h_rnf_tm(ji,jj)= h_rnf_tm(ji,jj)+h_rnf(ji,jj)
+                  h_rnf(ji,jj)=h_rnf_tm(ji,jj)/(nn_dttrc+1)
+
+                  hmld_tm(ji,jj)= hmld_tm(ji,jj)+hmld(ji,jj)
+                  hmld(ji,jj)=hmld_tm(ji,jj)/(nn_dttrc+1)
+ 
+                  fr_i_tm(ji,jj)= fr_i_tm(ji,jj)+fr_i(ji,jj)
+                  fr_i(ji,jj)=fr_i_tm(ji,jj)/(nn_dttrc+1)
+
+!  variables that are initialized after averages initialized
+                  wndm_tm(ji,jj)= wndm_tm(ji,jj)+wndm(ji,jj)
+                  IF(kt == nittrc000)THEN
+                      wndm(ji,jj)=wndm_tm(ji,jj)/(nn_dttrc)
+                  ELSE
+                      wndm(ji,jj)=wndm_tm(ji,jj)/(nn_dttrc+1)
+                  ENDIF 
+
+                  qsr_tm(ji,jj)= qsr_tm(ji,jj)+qsr(ji,jj)
+                  IF(kt == nittrc000)THEN
+                      qsr(ji,jj)=qsr_tm(ji,jj)/(nn_dttrc)
+                  ELSE
+                      qsr(ji,jj)=qsr_tm(ji,jj)/(nn_dttrc+1)
+                  ENDIF 
+
+		  emp_b(ji,jj)=emp_b_hold(ji,jj)
+
+                  emp_tm(ji,jj)= emp_tm(ji,jj)+emp(ji,jj)
+                  IF(kt == nittrc000)THEN
+                     emp(ji,jj)=emp_tm(ji,jj)/(nn_dttrc)
+                  ELSE
+                     emp(ji,jj)=emp_tm(ji,jj)/(nn_dttrc+1)
+                  ENDIF 
+
+                  emps_tm(ji,jj)= emps_tm(ji,jj)+emps(ji,jj)
+                  IF(kt == nittrc000)THEN
+                     emps(ji,jj)=emps_tm(ji,jj)/(nn_dttrc)
+                  ELSE
+                     emps(ji,jj)=emps_tm(ji,jj)/(nn_dttrc+1)
+                  ENDIF
+
+               END DO
+            END DO
+         DO jk = 1, jpk
+             IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
+             !                                           !------------------------------!
+                 fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
+                 fse3u (:,:,jk) = fse3u_n (:,:,jk)
+                 fse3v (:,:,jk) = fse3v_n (:,:,jk)
+                 fse3f (:,:,jk) = fse3f_n (:,:,jk)
+                 fse3w (:,:,jk) = fse3w_n (:,:,jk)
+                 fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
+                 fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
+             ENDIF
+         END DO
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+
+                  un_tm(ji,jj,jk)= un_tm(ji,jj,jk)+un(ji,jj,jk)*e3u_temp(ji,jj,jk)
+                  un(ji,jj,jk)=(un_tm(ji,jj,jk)/(nn_dttrc+1))/fse3u(ji,jj,jk)
+		  
+                  vn_tm(ji,jj,jk)= vn_tm(ji,jj,jk)+vn(ji,jj,jk)*e3v_temp(ji,jj,jk)
+                  vn(ji,jj,jk)=(vn_tm(ji,jj,jk)/(nn_dttrc+1))/fse3v(ji,jj,jk)
+		  
+                  wn_tm(ji,jj,jk)= wn_tm(ji,jj,jk)+wn(ji,jj,jk)*e3w_temp(ji,jj,jk)
+                  wn(ji,jj,jk)=(wn_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)
+		  
+                  avt_tm(ji,jj,jk)= avt_tm(ji,jj,jk)+avt(ji,jj,jk)*e3t_temp(ji,jj,jk)
+                  avt(ji,jj,jk)=(avt_tm(ji,jj,jk)/(nn_dttrc+1))/fse3t(ji,jj,jk)	
+	  	  
+                  tsn_tm(ji,jj,jk,jp_tem)= tsn_tm(ji,jj,jk,jp_tem)+tsn(ji,jj,jk,jp_tem)*e3t_temp(ji,jj,jk)
+                  tsn(ji,jj,jk,jp_tem)=(tsn_tm(ji,jj,jk,jp_tem)/(nn_dttrc+1))/fse3t(ji,jj,jk)
+	  
+                  tsn_tm(ji,jj,jk,jp_sal)= tsn_tm(ji,jj,jk,jp_sal)+tsn(ji,jj,jk,jp_sal)*e3t_temp(ji,jj,jk)
+                  tsn(ji,jj,jk,jp_sal)=(tsn_tm(ji,jj,jk,jp_sal)/(nn_dttrc+1))/fse3t(ji,jj,jk)
+	
+# if defined key_zdfddm
+                  avs_tm(ji,jj,jk)= avs_tm(ji,jj,jk)+avs(ji,jj,jk)*e3w_temp(ji,jj,jk)
+                  avs(ji,jj,jk)=(avs_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)
+
+# endif
+	  
+#if defined key_ldfslp
+                  wslpi_tm(ji,jj,jk)= wslpi_tm(ji,jj,jk)+wslpi(ji,jj,jk)*e3w_temp(ji,jj,jk)
+                  wslpi(ji,jj,jk)=(wslpi_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)
+
+                  wslpj_tm(ji,jj,jk)= wslpj_tm(ji,jj,jk)+wslpj(ji,jj,jk)*e3w_temp(ji,jj,jk)
+                  wslpj(ji,jj,jk)=(wslpj_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)
+
+                  uslp_tm(ji,jj,jk)= uslp_tm(ji,jj,jk)+uslp(ji,jj,jk)*e3u_temp(ji,jj,jk)
+                  uslp(ji,jj,jk)=(uslp_tm(ji,jj,jk)/(nn_dttrc+1))/fse3u(ji,jj,jk)
+
+                  vslp_tm(ji,jj,jk)= vslp_tm(ji,jj,jk)+vslp(ji,jj,jk)*e3v_temp(ji,jj,jk)
+                  vslp(ji,jj,jk)=(vslp_tm(ji,jj,jk)/(nn_dttrc+1))/fse3v(ji,jj,jk)
+#endif
+	       END DO
+            END DO
+         END DO
+#if defined key_traldf_c3d
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ahtt_tm(ji,jj,jk)= ahtt_tm(ji,jj,jk)+ahtt(ji,jj,jk)*fse3t(ji,jj,jk)
+                  ahtt(ji,jj,jk)=(ahtt_tm(ji,jj,jk)/(nn_dttrc+1))/fse3t(ji,jj,jk)	
+
+                  ahtu_tm(ji,jj,jk)= ahtu_tm(ji,jj,jk)+ahtu(ji,jj,jk)*fse3u(ji,jj,jk)
+                  ahtu(ji,jj,jk)=(ahtu_tm(ji,jj,jk)/(nn_dttrc+1))/fse3u(ji,jj,jk)	
+
+                  ahtv_tm(ji,jj,jk)= ahtv_tm(ji,jj,jk)+ahtv(ji,jj,jk)*fse3v(ji,jj,jk)
+                  ahtv(ji,jj,jk)=(ahtv_tm(ji,jj,jk)/(nn_dttrc+1))/fse3v(ji,jj,jk)	
+
+                  ahtw_tm(ji,jj,jk)= ahtw_tm(ji,jj,jk)+ahtw(ji,jj,jk)*fse3w(ji,jj,jk)
+                  ahtw(ji,jj,jk)=(ahtw_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)	
+
+              END DO
+            END DO
+         END DO
+#elif defined key_traldf_c2d
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ahtt_tm(ji,jj)= ahtt_tm(ji,jj)+ahtt(ji,jj)
+                  ahtt(ji,jj)=ahtt_tm(ji,jj)/(nn_dttrc+1)
+
+                  ahtu_tm(ji,jj)= ahtu_tm(ji,jj)+ahtu(ji,jj)
+                  ahtu(ji,jj)=ahtu_tm(ji,jj)/(nn_dttrc+1)
+
+                  ahtv_tm(ji,jj)= ahtv_tm(ji,jj)+ahtv(ji,jj)
+                  ahtv(ji,jj)=ahtv_tm(ji,jj)/(nn_dttrc+1)
+
+                  ahtw_tm(ji,jj)= ahtw_tm(ji,jj)+ahtw(ji,jj)
+                  ahtw(ji,jj)=ahtw_tm(ji,jj)/(nn_dttrc+1)
+
+             END DO
+            END DO
+#elif defined key_traldf_c1d
+          DO jk = 1, jpk
+             ahtt_tm(jk)= ahtt_tm(jk)+ahtt(jk)
+             ahtt(jk)=ahtt_tm(jk)/(nn_dttrc+1)
+
+             ahtu_tm(jk)= ahtu_tm(jk)+ahtu(jk)
+             ahtu(jk)=ahtu_tm(jk)/(nn_dttrc+1)
+
+             ahtv_tm(jk)= ahtv_tm(jk)+ahtv(jk)
+             ahtv(jk)=ahtv_tm(jk)/(nn_dttrc+1)
+
+             ahtw_tm(jk)= ahtw_tm(jk)+ahtw(jk)
+             ahtw(jk)=ahtw_tm(jk)/(nn_dttrc+1)
+          END DO
+#else
+             ahtt_tm= ahtt_tm+ahtt
+             ahtt=ahtt_tm/(nn_dttrc+1)
+
+             ahtu_tm= ahtu_tm+ahtu
+             ahtu=ahtu_tm/(nn_dttrc+1)
+
+             ahtv_tm= ahtv_tm+ahtv
+             ahtv=ahtv_tm/(nn_dttrc+1)
+
+             ahtw_tm= ahtw_tm+ahtw
+             ahtw=ahtw_tm/(nn_dttrc+1)
+
+#endif
+#if defined key_traldf_eiv
+# if defined key_traldf_c3d
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  aeiu_tm(ji,jj,jk)= aeiu_tm(ji,jj,jk)+aeiu(ji,jj,jk)*fse3u(ji,jj,jk)
+                  aeiu(ji,jj,jk)=(aeiu_tm(ji,jj,jk)/(nn_dttrc+1))/fse3u(ji,jj,jk)	
+
+                  aeiv_tm(ji,jj,jk)= aeiv_tm(ji,jj,jk)+aeiv(ji,jj,jk)*fse3v(ji,jj,jk)
+                  aeiv(ji,jj,jk)=(aeiv_tm(ji,jj,jk)/(nn_dttrc+1))/fse3v(ji,jj,jk)	
+
+                  aeiw_tm(ji,jj,jk)= aeiw_tm(ji,jj,jk)+aeiw(ji,jj,jk)*fse3w(ji,jj,jk)
+                  aeiw(ji,jj,jk)=(aeiw_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)	
+
+             END DO
+            END DO
+         END DO
+# elif defined key_traldf_c2d
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  aeiu_tm(ji,jj)= aeiu_tm(ji,jj)+aeiu(ji,jj)
+                  aeiu(ji,jj)=aeiu_tm(ji,jj)/(nn_dttrc+1)
+
+                  aeiv_tm(ji,jj)= aeiv_tm(ji,jj)+aeiv(ji,jj)
+                  aeiv(ji,jj)=aeiv_tm(ji,jj)/(nn_dttrc+1)
+
+                  aeiw_tm(ji,jj)= aeiw_tm(ji,jj)+aeiw(ji,jj)
+                  aeiw(ji,jj)=aeiw_tm(ji,jj)/(nn_dttrc+1)
+             END DO
+            END DO
+# elif defined key_traldf_c1d
+          DO jk = 1, jpk
+             aeiu_tm(jk)= aeiu_tm(jk)+aeiu(jk)
+             aeiu(jk)=aeiu_tm(jk)/(nn_dttrc+1)
+
+             aeiv_tm(jk)= aeiv_tm(jk)+aeiv(jk)
+             aeiv(jk)=aeiv_tm(jk)/(nn_dttrc+1)
+
+             aeiw_tm(jk)= aeiw_tm(jk)+aeiw(jk)
+             aeiw(jk)=aeiw_tm(jk)/(nn_dttrc+1)
+          END DO
+# else
+             aeiu_tm= aeiu_tm+aeiu
+             aeiu=aeiu_tm/(nn_dttrc+1)
+
+             aeiv_tm= aeiv_tm+aeiv
+             aeiv=aeiv_tm/(nn_dttrc+1)
+
+             aeiw_tm= aeiw_tm+aeiw
+             aeiw=aeiw_tm/(nn_dttrc+1)
+# endif
+# endif
+
+         CALL lbc_lnk( sshn(:,:), 'T', 1. ) 
+         CALL lbc_lnk( sshu_n(:,:), 'U', 1. ) 
+         CALL lbc_lnk( sshv_n(:,:), 'V', 1. ) 
+         CALL lbc_lnk( sshf_n(:,:), 'F', 1. ) 
+         CALL lbc_lnk( sshb(:,:), 'T', 1. ) 
+         CALL lbc_lnk( sshu_b(:,:), 'U', 1. ) 
+         CALL lbc_lnk( sshv_b(:,:), 'V', 1. ) 
+         CALL lbc_lnk( ssha(:,:), 'T', 1. ) 
+         CALL lbc_lnk( sshu_a(:,:), 'U', 1. ) 
+         CALL lbc_lnk( sshv_a(:,:), 'V', 1. ) 
+         CALL lbc_lnk( rnf(:,:), 'T', 1. ) 
+         CALL lbc_lnk( h_rnf(:,:), 'T', 1. ) 
+         CALL lbc_lnk( rotn(:,:,:), 'F', 1. ) 
+         CALL lbc_lnk( hdivn(:,:,:), 'T', 1. ) 
+         CALL lbc_lnk( rotb(:,:,:), 'F', 1. ) 
+         CALL lbc_lnk( hdivb(:,:,:), 'T', 1. ) 
+         CALL lbc_lnk( un(:,:,:), 'U', 1. ) 
+         CALL lbc_lnk( vn(:,:,:), 'V', 1. ) 
+         CALL lbc_lnk( wn(:,:,:), 'V', 1. ) 
+         CALL lbc_lnk( tsn(:,:,:,jp_tem), 'T', 1. ) 
+         CALL lbc_lnk( tsn(:,:,:,jp_sal), 'T', 1. ) 
+         CALL lbc_lnk( avt(:,:,:), 'T', 1. ) 
+         CALL lbc_lnk( hmld(:,:), 'T', 1. ) 
+         CALL lbc_lnk( qsr(:,:), 'T', 1. ) 
+         CALL lbc_lnk( wndm(:,:), 'T', 1. ) 
+         CALL lbc_lnk( emp(:,:), 'T', 1. ) 
+         CALL lbc_lnk( emp_b(:,:), 'T', 1. ) 
+         CALL lbc_lnk( emps(:,:), 'T', 1. ) 
+#if defined key_traldf_c3d
+         CALL lbc_lnk( ahtt(:,:,:), 'T', 1. ) 
+         CALL lbc_lnk( ahtu(:,:,:), 'U', 1. ) 
+         CALL lbc_lnk( ahtv(:,:,:), 'V', 1. ) 
+         CALL lbc_lnk( ahtw(:,:,:), 'W', 1. ) 
+#elif defined key_traldf_c2d
+         CALL lbc_lnk( ahtt(:,:), 'T', 1. ) 
+         CALL lbc_lnk( ahtu(:,:), 'U', 1. ) 
+         CALL lbc_lnk( ahtv(:,:), 'V', 1. ) 
+         CALL lbc_lnk( ahtw(:,:), 'W', 1. ) 
+#endif
+# if defined key_zdfddm
+          CALL lbc_lnk( avs, 'W', 1. ) 
+# endif
+#if defined key_traldf_eiv
+#if defined key_traldf_c3d
+         CALL lbc_lnk( aeiu(:,:,:), 'U', 1. ) 
+         CALL lbc_lnk( aeiv(:,:,:), 'V', 1. ) 
+         CALL lbc_lnk( aeiw(:,:,:), 'W', 1. ) 
+#elif defined key_traldf_c2d
+         CALL lbc_lnk( aeiu(:,:), 'U', 1. ) 
+         CALL lbc_lnk( aeiv(:,:), 'V', 1. ) 
+         CALL lbc_lnk( aeiw(:,:), 'W', 1. ) 
+#endif
+#endif
+#if defined key_ldfslp
+         CALL lbc_lnk( uslp(:,:,:), 'U', 1. ) 
+         CALL lbc_lnk( vslp(:,:,:), 'V', 1. ) 
+         CALL lbc_lnk( wslpi(:,:,:), 'W', 1. ) 
+         CALL lbc_lnk( wslpj(:,:,:), 'W', 1. ) 
+#endif
+         CALL lbc_lnk( hu(:,:), 'U', 1. ) 
+         CALL lbc_lnk( hv(:,:), 'V', 1. ) 
+         CALL lbc_lnk( hur(:,:), 'U', 1. ) 
+         CALL lbc_lnk( hvr(:,:), 'V', 1. ) 
+
+         CALL trc_sub_ssh( kt )         ! after ssh & vertical velocity
+      ENDIF
+
+
+   END SUBROUTINE trc_sub_stp
+  !!-------------------------------------------------------------------
+
+  !!-------------------------------------------------------------------
+   SUBROUTINE trc_sub_ini
+      !!-------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_sub_ini  ***
+      !!                      
+      !! ** Purpose : Initialize variables needed for sub-stepping passive tracers
+      !! 
+      !! ** Method  : 
+      !!              Compute the averages for sub-stepping
+      !!-------------------------------------------------------------------
+      INTEGER ::   jk, jn ,ji,jj   ! dummy loop indices
+      INTEGER ::   ierr
+      !!-------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'trc_sub_ini : initial set up of the passive tracers substepping'
+      IF(lwp) WRITE(numout,*) '~~~~~~~'
+
+      ierr =  trc_sub_alloc    ()
+      IF( lk_mpp    )   CALL mpp_sum( ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'top_sub_alloc : unable to allocate standard ocean arrays' )
+
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               un_tm(ji,jj,jk)= un(ji,jj,jk)*fse3u(ji,jj,jk)
+               vn_tm(ji,jj,jk)= vn(ji,jj,jk)*fse3v(ji,jj,jk)
+	       wn_tm(ji,jj,jk)= wn(ji,jj,jk)*fse3w(ji,jj,jk)
+               avt_tm(ji,jj,jk)=avt(ji,jj,jk)*fse3t(ji,jj,jk)
+               tsn_tm(ji,jj,jk,jp_tem)= tsn_tm(ji,jj,jk,jp_tem)+tsn(ji,jj,jk,jp_tem)*fse3t(ji,jj,jk)
+               tsn_tm(ji,jj,jk,jp_sal)= tsn_tm(ji,jj,jk,jp_sal)+tsn(ji,jj,jk,jp_sal)*fse3t(ji,jj,jk)
+#if defined key_zdfddm
+               avs_tm(ji,jj,jk)= avs_tm(ji,jj,jk)+avs(ji,jj,jk)*fse3w(ji,jj,jk)
+#endif
+#if defined key_ldfslp
+               wslpi_tm(ji,jj,jk)= wslpi(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
+               wslpj_tm(ji,jj,jk)= wslpj(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
+               uslp_tm(ji,jj,jk)= uslp(ji,jj,jk)*fse3u(ji,jj,jk) !0. 
+               vslp_tm(ji,jj,jk)= vslp(ji,jj,jk)*fse3v(ji,jj,jk) !0. 
+#endif
+           END DO
+         END DO
+      END DO
+
+!   physics variables that are set after initialization:
+      wndm_tm(:,:)=0.  !wndm(:,:)
+      qsr_tm(:,:)= 0.  !qsr(:,:)
+      emp_tm(:,:)= 0.  !emp(:,:)
+      emps_tm(:,:)=0.  !emps(:,:)
+ 
+      sshn_tm(:,:) = sshn(:,:) !0.
+      sshu_n_tm(:,:) = sshu_n(:,:) !0.
+      sshv_n_tm(:,:) = sshv_n(:,:) !0.
+      rnf_tm(:,:)= rnf(:,:) !0. 
+      h_rnf_tm(:,:)= h_rnf(:,:) !0.
+      hmld_tm(:,:)=hmld(:,:)
+      fr_i_tm(:,:)=fr_i(:,:)
+#if defined key_traldf_c3d
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ahtt_tm(ji,jj,jk)= ahtt(ji,jj,jk)*fse3t(ji,jj,jk)
+                  ahtu_tm(ji,jj,jk)= ahtu(ji,jj,jk)*fse3u(ji,jj,jk)
+                  ahtv_tm(ji,jj,jk)= ahtv(ji,jj,jk)*fse3v(ji,jj,jk)
+                  ahtw_tm(ji,jj,jk)= ahtw(ji,jj,jk)*fse3w(ji,jj,jk)
+             END DO
+            END DO
+         END DO
+#elif defined key_traldf_c2d
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ahtt_tm(ji,jj)= ahtt(ji,jj)
+                  ahtu_tm(ji,jj)= ahtu(ji,jj)
+                  ahtv_tm(ji,jj)= ahtv(ji,jj)
+                  ahtw_tm(ji,jj)= ahtw(ji,jj)
+             END DO
+            END DO
+#elif defined key_traldf_c1d
+          DO jk = 1, jpk
+             ahtt_tm(jk)= ahtt(jk)
+             ahtu_tm(jk)= ahtu(jk)
+             ahtv_tm(jk)= ahtv(jk)
+             ahtw_tm(jk)= ahtw(jk)
+          END DO
+#else
+             ahtt_tm= ahtt
+             ahtu_tm= ahtu
+             ahtv_tm= ahtv
+             ahtw_tm= ahtw
+#endif
+#if defined key_traldf_eiv
+# if defined key_traldf_c3d
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  aeiu_tm(ji,jj,jk)= aeiu(ji,jj,jk)*fse3u(ji,jj,jk)
+                  aeiv_tm(ji,jj,jk)= aeiv(ji,jj,jk)*fse3v(ji,jj,jk)
+                  aeiw_tm(ji,jj,jk)= aeiw(ji,jj,jk)*fse3w(ji,jj,jk)
+             END DO
+            END DO
+         END DO
+# elif defined key_traldf_c2d
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  aeiu_tm(ji,jj)= aeiu(ji,jj)
+                  aeiv_tm(ji,jj)= aeiv(ji,jj)
+                  aeiw_tm(ji,jj)= aeiw(ji,jj)
+             END DO
+            END DO
+# elif defined key_traldf_c1d
+          DO jk = 1, jpk
+             aeiu_tm(jk)= aeiu(jk)
+             aeiv_tm(jk)= aeiv(jk)
+             aeiw_tm(jk)= aeiw(jk)
+          END DO
+# else
+             aeiu_tm= aeiu
+             aeiv_tm= aeiv
+             aeiw_tm= aeiw
+# endif
+# endif
+
+   END SUBROUTINE trc_sub_ini
+
+  !!-------------------------------------------------------------------
+   SUBROUTINE trc_sub_reset(kt)
+      !!-------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_sub_reset  ***
+      !!                      
+      !! ** Purpose : Reset physics variables averaged for substepping
+      !! 
+      !! ** Method  : 
+      !!              Compute the averages for sub-stepping
+      !!-------------------------------------------------------------------
+      INTEGER ::   jk, jn ,ji,jj    ! dummy loop indices
+      INTEGER ::   ierr
+      INTEGER, INTENT( in ) ::  kt  ! ocean time-step index
+      !!-------------------------------------------------------------------
+
+!   restore physics variables
+        avt(:,:,:)=avt_temp(:,:,:)
+        un(:,:,:)=un_temp(:,:,:)
+        vn(:,:,:)=vn_temp(:,:,:)
+        wn(:,:,:)=wn_temp(:,:,:)
+        sshn(:,:)=sshn_temp(:,:)
+        sshb(:,:)=sshb_temp(:,:)
+        ssha(:,:)=ssha_temp(:,:)
+        sshu_n(:,:)=sshu_n_temp(:,:)
+        sshu_b(:,:)=sshu_b_temp(:,:)
+        sshu_a(:,:)=sshu_a_temp(:,:)
+        sshv_n(:,:)=sshv_n_temp(:,:)
+        sshv_b(:,:)=sshv_b_temp(:,:)
+        sshv_a(:,:)=sshv_a_temp(:,:)
+        sshf_n(:,:)=sshf_n_temp(:,:)
+#if defined key_ldfslp
+        wslpi(:,:,:)=wslpi_temp(:,:,:)
+        wslpj(:,:,:)=wslpj_temp(:,:,:)
+        uslp(:,:,:)=uslp_temp(:,:,:)
+        vslp(:,:,:)=vslp_temp(:,:,:)
+#endif
+        rnf(:,:)=rnf_temp(:,:)
+        h_rnf(:,:)=h_rnf_temp(:,:)
+        hdivn(:,:,:)=hdivn_temp(:,:,:)
+        rotn(:,:,:)=rotn_temp(:,:,:)
+        hdivb(:,:,:)=hdivb_temp(:,:,:)
+        rotb(:,:,:)=rotb_temp(:,:,:)
+        hu(:,:)=hu_temp(:,:)
+        hv(:,:)=hv_temp(:,:)
+        hur(:,:)=hur_temp(:,:)
+        hvr(:,:)=hvr_temp(:,:)
+        tsn(:,:,:,:)=tsn_temp(:,:,:,:)
+# if defined key_zdfddm
+        avs(:,:,:)=avs_temp(:,:,:)
+# endif
+       hmld(:,:)=hmld_temp(:,:)
+       qsr(:,:)=qsr_temp(:,:)
+       fr_i(:,:)=fr_i_temp(:,:)
+       wndm(:,:)=wndm_temp(:,:)
+       emp(:,:)=emp_temp(:,:)
+       emps(:,:)=emps_temp(:,:)
+       emp_b(:,:)=emp_b_temp(:,:)
+#if defined key_traldf_c3d
+            ahtu_temp(:,:,:)=ahtu(:,:,:)
+            ahtv_temp(:,:,:)=ahtv(:,:,:)
+            ahtw_temp(:,:,:)=ahtw(:,:,:)
+            ahtt_temp(:,:,:)=ahtt(:,:,:)
+#elif defined key_traldf_c2d
+            ahtu_temp(:,:)=ahtu(:,:)
+            ahtv_temp(:,:)=ahtv(:,:)
+            ahtw_temp(:,:)=ahtw(:,:)
+            ahtt_temp(:,:)=ahtt(:,:)
+#elif defined key_traldf_c1d
+            ahtu_temp(:)=ahtu(:)
+            ahtv_temp(:)=ahtv(:)
+            ahtw_temp(:)=ahtw(:)
+            ahtt_temp(:)=ahtt(:)
+#else
+            ahtu_temp=ahtu
+            ahtv_temp=ahtv
+            ahtw_temp=ahtw
+            ahtt_temp=ahtt
+#endif
+#if defined key_traldf_eiv
+#if defined key_traldf_c3d
+            aeiu_temp(:,:,:)=aeiu(:,:,:)
+            aeiv_temp(:,:,:)=aeiv(:,:,:)
+            aeiw_temp(:,:,:)=aeiw(:,:,:)
+#elif defined key_traldf_c2d
+            aeiu_temp(:,:)=aeiu(:,:)
+            aeiv_temp(:,:)=aeiv(:,:)
+            aeiw_temp(:,:)=aeiw(:,:)
+#elif defined key_traldf_c1d
+            aeiu_temp(:)=aeiu(:)
+            aeiv_temp(:)=aeiv(:)
+            aeiw_temp(:)=aeiw(:)
+#else
+            aeiu_temp=aeiu
+            aeiv_temp=aeiv
+            aeiw_temp=aeiw
+#endif
+#endif        
+
+
+       !                                           !------------------------------!
+       IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
+       !                                           !------------------------------!
+         DO jk = 1, jpkm1
+            fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays
+            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
+            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
+            !
+            fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
+            fse3u (:,:,jk) = fse3u_n (:,:,jk)
+            fse3v (:,:,jk) = fse3v_n (:,:,jk)
+            fse3f (:,:,jk) = fse3f_n (:,:,jk)
+            fse3w (:,:,jk) = fse3w_n (:,:,jk)
+            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
+            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
+         END DO
+      ENDIF
+
+!   start new averages
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  un_tm(ji,jj,jk)=un(ji,jj,jk)*fse3u(ji,jj,jk)  !0.0
+                  vn_tm(ji,jj,jk)=vn(ji,jj,jk)*fse3v(ji,jj,jk)  !0.0
+                  wn_tm(ji,jj,jk)=wn(ji,jj,jk)*fse3w(ji,jj,jk)  !0.0
+                  avt_tm(ji,jj,jk)=avt(ji,jj,jk)*fse3t(ji,jj,jk)  !0.0
+                  tsn_tm(ji,jj,jk,jp_tem)=tsn(ji,jj,jk,jp_tem)*fse3t(ji,jj,jk)  !0.0
+                  tsn_tm(ji,jj,jk,jp_sal)=tsn(ji,jj,jk,jp_sal)*fse3t(ji,jj,jk)  !0.0
+# if defined key_zdfddm
+                  avs_tm==avs(ji,jj,jk)*fse3w(ji,jj,jk)  !0.0
+# endif
+#if defined key_ldfslp
+                  wslpi_tm(ji,jj,jk)= wslpi(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
+                  wslpj_tm(ji,jj,jk)=wslpj(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
+                  uslp_tm(ji,jj,jk)=uslp(ji,jj,jk)*fse3u(ji,jj,jk) !0. 
+                  vslp_tm(ji,jj,jk)=vslp(ji,jj,jk)*fse3v(ji,jj,jk) !0. 
+#endif
+ 	       END DO
+            END DO
+         END DO
+      sshn_tm(:,:) = sshn(:,:) !0.
+      sshu_n_tm(:,:) = sshu_n(:,:) !0.
+      sshv_n_tm(:,:) = sshv_n(:,:) !0.
+      rnf_tm(:,:)= rnf(:,:) !0. 
+      h_rnf_tm(:,:)= h_rnf(:,:) !0.
+      hmld_tm(:,:)=hmld(:,:)
+      fr_i_tm(:,:)=fr_i(:,:)
+      sshb_hold(:,:)=sshn(:,:)
+      sshu_b_hold(:,:)=sshu_n(:,:)
+      sshv_b_hold(:,:)=sshv_n(:,:)
+      emp_b_hold(:,:)=emp(:,:)
+      emp_tm(:,:)=emp(:,:)
+      emps_tm(:,:)=emps(:,:)
+      qsr_tm(:,:)=qsr(:,:)
+      wndm_tm(:,:)=wndm(:,:)
+#if defined key_traldf_c3d
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ahtt_tm(ji,jj,jk)= ahtt(ji,jj,jk)*fse3t(ji,jj,jk)
+                  ahtu_tm(ji,jj,jk)= ahtu(ji,jj,jk)*fse3u(ji,jj,jk)
+                  ahtv_tm(ji,jj,jk)= ahtv(ji,jj,jk)*fse3v(ji,jj,jk)
+                  ahtw_tm(ji,jj,jk)= ahtw(ji,jj,jk)*fse3w(ji,jj,jk)
+             END DO
+            END DO
+         END DO
+#elif defined key_traldf_c2d
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ahtt_tm(ji,jj)= ahtt(ji,jj)
+                  ahtu_tm(ji,jj)= ahtu(ji,jj)
+                  ahtv_tm(ji,jj)= ahtv(ji,jj)
+                  ahtw_tm(ji,jj)= ahtw(ji,jj)
+             END DO
+            END DO
+#elif defined key_traldf_c1d
+          DO jk = 1, jpk
+             ahtt_tm(jk)= ahtt(jk)
+             ahtu_tm(jk)= ahtu(jk)
+             ahtv_tm(jk)= ahtv(jk)
+             ahtw_tm(jk)= ahtw(jk)
+          END DO
+#else
+             ahtt_tm= ahtt
+             ahtu_tm= ahtu
+             ahtv_tm= ahtv
+             ahtw_tm= ahtw
+#endif
+#if defined key_traldf_eiv
+# if defined key_traldf_c3d
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  aeiu_tm(ji,jj,jk)= aeiu(ji,jj,jk)*fse3u(ji,jj,jk)
+                  aeiv_tm(ji,jj,jk)= aeiv(ji,jj,jk)*fse3v(ji,jj,jk)
+                  aeiw_tm(ji,jj,jk)= aeiw(ji,jj,jk)*fse3w(ji,jj,jk)
+             END DO
+            END DO
+         END DO
+# elif defined key_traldf_c2d
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  aeiu_tm(ji,jj)= aeiu(ji,jj)
+                  aeiv_tm(ji,jj)= aeiv(ji,jj)
+                  aeiw_tm(ji,jj)= aeiw(ji,jj)
+             END DO
+            END DO
+# elif defined key_traldf_c1d
+          DO jk = 1, jpk
+             aeiu_tm(jk)= aeiu(jk)
+             aeiv_tm(jk)= aeiv(jk)
+             aeiw_tm(jk)= aeiw(jk)
+          END DO
+# else
+             aeiu_tm= aeiu
+             aeiv_tm= aeiv
+             aeiw_tm= aeiw
+# endif
+# endif
+
+
+   END SUBROUTINE trc_sub_reset
+  !!-------------------------------------------------------------------
+   SUBROUTINE trc_sub_ssh( kt ) 
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE trc_sub_ssh  ***
+      !!                   
+      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
+      !!              and update the now vertical coordinate (lk_vvl=T).
+      !!
+      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
+      !!      velocity is computed by integrating the horizontal divergence  
+      !!      from the bottom to the surface minus the scale factor evolution.
+      !!        The boundary conditions are w=0 at the bottom (no flux) and.
+      !!
+      !! ** action  :   ssha    : after sea surface height
+      !!                wn      : now vertical velocity
+      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
+      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE oce     , ONLY:   z3d   => ta                           ! ta used as 3D workspace
+      USE wrk_nemo, ONLY:   zhdiv => wrk_2d_1 , z2d => wrk_2d_2   ! 2D workspace
+      !
+      INTEGER, INTENT(in) ::   kt   ! time step
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(2, 1,2) ) THEN
+         CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable')   ;   RETURN
+      ENDIF
+
+      IF( kt == nittrc000 ) THEN
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'trc_sub_ssh : after sea surface height and now vertical velocity '
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
+         !
+         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+         !
+         !
+      ENDIF
+
+      !                                           !------------------------------------------!
+      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
+         !                                        !------------------------------------------!
+         DO jk = 1, jpkm1
+            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
+            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
+            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
+            !
+            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
+            fse3u (:,:,jk) = fse3u_n (:,:,jk)
+            fse3v (:,:,jk) = fse3v_n (:,:,jk)
+            fse3f (:,:,jk) = fse3f_n (:,:,jk)
+            fse3w (:,:,jk) = fse3w_n (:,:,jk)
+            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
+            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
+         END DO
+         !
+         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
+         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
+         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
+         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
+         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )
+         ! 
+      ENDIF
+      !
+      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
+      !
+      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
+      IF( neuler == 0 .AND. kt == nittrc000 )   z2dt = rdt
+
+      !                                           !------------------------------!
+      !                                           !   After Sea Surface Height   !
+      !                                           !------------------------------!
+      zhdiv(:,:) = 0._wp
+      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
+        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
+      END DO
+      !                                                ! Sea surface elevation time stepping
+      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
+      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
+      z1_rau0 = 0.5 / rau0
+      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
+
+#if defined key_agrif
+      CALL agrif_ssh( kt )
+#endif
+#if defined key_obc
+      IF( Agrif_Root() ) THEN 
+         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
+         CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm)
+      ENDIF
+#endif
+#if defined key_bdy
+      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
+      CALL lbc_lnk( ssha, 'T', 1. ) 
+#endif
+
+      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
+      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1      ! NO Vector Opt.
+               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
+                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
+                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
+               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
+                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
+                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
+            END DO
+         END DO
+         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
+      ENDIF
+      
+
+      !                                           !------------------------------!
+      !                                           !     Now Vertical Velocity    !
+      !                                           !------------------------------!
+      z1_2dt = 1.e0 / z2dt
+      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
+         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
+         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
+            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
+            &                         * tmask(:,:,jk) * z1_2dt
+#if defined key_bdy
+         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
+#endif
+      END DO
+
+      !
+      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('ssh_wzv: failed to release workspace arrays')
+      !
+   END SUBROUTINE trc_sub_ssh
+  !!-------------------------------------------------------------------
+
+  !!-------------------------------------------------------------------
+   INTEGER FUNCTION trc_sub_alloc()
+      !!-------------------------------------------------------------------
+      !!                    *** ROUTINE trc_sub_alloc ***
+      !!-------------------------------------------------------------------
+      USE lib_mpp, ONLY: ctl_warn
+      INTEGER ::         ierr
+      !!-------------------------------------------------------------------
+      !
+      ALLOCATE( un_temp(jpi,jpj,jpk)        ,  vn_temp(jpi,jpj,jpk)  ,   &
+         &      wn_temp(jpi,jpj,jpk)        ,  avt_temp(jpi,jpj,jpk) ,   &
+         &      sshn_temp(jpi,jpj)          ,  sshb_temp(jpi,jpj) ,      &
+         &      ssha_temp(jpi,jpj)          ,  sshu_a_temp(jpi,jpj),     &
+         &      sshu_n_temp(jpi,jpj)        ,  sshu_b_temp(jpi,jpj),     &
+         &      sshv_n_temp(jpi,jpj)        ,  sshv_b_temp(jpi,jpj),     &
+         &      sshv_a_temp(jpi,jpj)        ,  sshf_a_temp(jpi,jpj),     &
+         &      sshf_n_temp(jpi,jpj)        ,  sshf_b_temp(jpi,jpj),     &
+#if defined key_ldfslp
+         &      wslpi_temp(jpi,jpj,jpk)     ,  wslpj_temp(jpi,jpj,jpk),  &
+         &      uslp_temp(jpi,jpj,jpk)      ,  vslp_temp(jpi,jpj,jpk),   &
+#endif
+         &      rnf_temp(jpi,jpj)           ,  h_rnf_temp(jpi,jpj) ,     &
+         &      tsn_temp(jpi,jpj,jpk,2)     ,                            &
+         &      emp_temp(jpi,jpj)           ,  emps_temp(jpi,jpj) ,      &
+         &      hmld_temp(jpi,jpj)          ,  qsr_temp(jpi,jpj) ,       &
+         &      fr_i_temp(jpi,jpj)          ,  fr_i_tm(jpi,jpj) ,        &
+         &      wndm_temp(jpi,jpj)          ,  wndm_tm(jpi,jpj) ,        &
+# if defined key_zdfddm
+         &      avt_tm(jpi,jpj,jpk)         ,  avt_temp(jpi,jpj,jpk) ,   &
+# endif
+#if defined key_traldf_c3d
+         &      ahtt_tm(jpi,jpj,jpk)        ,  ahtt_temp(jpi,jpj,jpk),   &
+         &      ahtu_tm(jpi,jpj,jpk)        ,  ahtu_temp(jpi,jpj,jpk),   &
+         &      ahtv_tm(jpi,jpj,jpk)        ,  ahtv_temp(jpi,jpj,jpk),   &
+         &      ahtw_tm(jpi,jpj,jpk)        ,  ahtw_temp(jpi,jpj,jpk),   &
+#elif defined key_traldf_c2d
+         &      ahtt_tm(jpi,jpj)            ,  ahtt_temp(jpi,jpj),       &
+         &      ahtu_tm(jpi,jpj)            ,  ahtu_temp(jpi,jpj),       &
+         &      ahtv_tm(jpi,jpj)            ,  ahtv_temp(jpi,jpj),       &
+         &      ahtw_tm(jpi,jpj)            ,  ahtw_temp(jpi,jpj),       &
+#elif defined key_traldf_c1d
+         &      ahtt_tm(jpk)                ,  ahtt_temp(jpk),           &
+         &      ahtu_tm(jpk)                ,  ahtu_temp(jpk),           &
+         &      ahtv_tm(jpk)                ,  ahtv_temp(jpk),           &
+         &      ahtw_tm(jpk)                ,  ahtw_temp(jpk),           &
+#endif
+#if defined key_traldf_eiv
+# if defined key_traldf_c3d
+         &      aeiu_tm(jpi,jpj,jpk)        ,  aeiu_temp(jpi,jpj,jpk),   &
+         &      aeiv_tm(jpi,jpj,jpk)        ,  aeiv_temp(jpi,jpj,jpk),   &
+         &      aeiw_tm(jpi,jpj,jpk)        ,  aeiw_temp(jpi,jpj,jpk),   &
+# elif defined key_traldf_c2d
+         &      aeiu_tm(jpi,jpj)            ,  aeiu_temp(jpi,jpj),       &
+         &      aeiv_tm(jpi,jpj)            ,  aeiv_temp(jpi,jpj),       &
+         &      aeiw_tm(jpi,jpj)            ,  aeiw_temp(jpi,jpj),       &
+# elif defined key_traldf_c1d
+         &      aeiu_tm(jpk)                ,  aeiu_temp(jpk),           &
+         &      aeiv_tm(jpk)                ,  aeiv_temp(jpk),           &
+         &      aeiw_tm(jpk)                ,  aeiw_temp(jpk),           &
+# endif
+# endif
+         &      hdivn_temp(jpi,jpj,jpk)     ,  hdivb_temp(jpi,jpj,jpk),  &
+         &      rotn_temp(jpi,jpj,jpk)      ,  rotb_temp(jpi,jpj,jpk),   &
+         &      hu_temp(jpi,jpj)            ,  hv_temp(jpi,jpj),         &
+         &      hur_temp(jpi,jpj)           ,  hvr_temp(jpi,jpj),        &
+         &      un_tm(jpi,jpj,jpk)          ,  vn_tm(jpi,jpj,jpk)  ,     &
+         &      wn_tm(jpi,jpj,jpk)          ,  avt_tm(jpi,jpj,jpk) ,     &
+         &      sshn_tm(jpi,jpj)            ,  sshb_hold(jpi,jpj) ,      &
+         &      sshu_n_tm(jpi,jpj)          ,  sshu_b_hold(jpi,jpj),     &
+         &      sshv_n_tm(jpi,jpj)          ,  sshv_b_hold(jpi,jpj),     &
+         &      tsn_tm(jpi,jpj,jpk,2)       ,                            &
+         &      emp_tm(jpi,jpj)             ,  emps_tm(jpi,jpj) ,        &
+         &      emp_b_hold(jpi,jpj)         ,                            &
+         &      hmld_tm(jpi,jpj)            ,  qsr_tm(jpi,jpj) ,         &
+#if defined key_ldfslp
+         &      wslpi_tm(jpi,jpj,jpk)       ,  wslpj_tm(jpi,jpj,jpk),    &
+         &      uslp_tm(jpi,jpj,jpk)        ,  vslp_tm(jpi,jpj,jpk),     &
+#endif
+         &      rnf_tm(jpi,jpj)             ,  h_rnf_tm(jpi,jpj) ,       &
+         &                                    STAT=trc_sub_alloc )  
+      IF( trc_sub_alloc /= 0 )   CALL ctl_warn('trc_sub_alloc: failed to allocate arrays')
+
+      !
+   END FUNCTION trc_sub_alloc
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default key                                     NO passive tracers
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_sub_stp( kt )        ! Empty routine
+      WRITE(*,*) 'trc_sub_stp: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_sub_stp
+   SUBROUTINE trc_sub_ini        ! Empty routine
+      WRITE(*,*) 'trc_sub_ini: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_sub_ini
+
+#endif
+
+   !!======================================================================
+END MODULE trcsub
Index: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcsubstp.F90
===================================================================
--- branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/trcsubstp.F90	(revision 2891)
+++ 	(revision )
@@ -1,399 +1,0 @@
-MODULE trcsubstp
-   !!======================================================================
-   !!                       ***  MODULE trcsubstp  ***
-   !! Time-stepping    : time loop of opa for passive tracer
-   !!======================================================================
-   !! History :  1.0  !  2004-03  (C. Ethe)  Original
-   !!----------------------------------------------------------------------
-#if defined key_top
-   !!----------------------------------------------------------------------
-   !!   trc_sub_stp     : passive tracer system sub-stepping 
-   !!----------------------------------------------------------------------
-   USE oce_trc          ! ocean dynamics and active tracers variables
-   USE trc
-   USE prtctl_trc       ! Print control for debbuging
-   USE iom
-   USE in_out_manager
-#if defined key_zdftke
-   USE zdftke          ! twice TKE (en)
-#endif
-#if defined key_zdfgls
-   USE zdfgls, ONLY: en
-#endif
-   USE zdf_oce
-   USE domvvl
-   USE sshwzv_trc
-   USE sbcrnf, ONLY:h_rnf
-
-   IMPLICIT NONE
-
-   PUBLIC   trc_sub_stp    ! called by trc_stp
-   PUBLIC   trc_sub_stp_ini ! called by trc_ini to initialize substepping arrays.
-   PUBLIC   trc_sub_stp_reset  !called by trc_stp to reset physics variables
-#  include "domzgr_substitute.h90"
-   
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Id: trcstp.F90 2528 2010-12-27 17:33:53Z rblod $ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE trc_sub_stp( kt )
-      !!-------------------------------------------------------------------
-      !!                     ***  ROUTINE trc_stp  ***
-      !!                      
-      !! ** Purpose : Average variables needed for sub-stepping passive tracers
-      !! 
-      !! ** Method  : 
-      !!              Compute the averages for sub-stepping
-      !!-------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::  kt  ! ocean time-step index
-      INTEGER::				    ji,jj,jk
-      !!-------------------------------------------------------------------
-
-       IF (kt .eq. nit000) THEN
-           sshb_hold=sshn
-           sshu_b_hold=sshu_n
-           sshv_b_hold=sshv_n
-           IF(lwp) WRITE(numout,*)
-           IF(lwp) WRITE(numout,*) 'trc_sub_stp : substepping of the passive tracers'
-           IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-       ENDIF  
-       IF( MOD( kt , nn_dttrc ) /= 0 ) THEN
-         if(lwp)write(numout,*)'trc_sub_stp at kt=',kt
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  un_tm(ji,jj,jk)= un_tm(ji,jj,jk)+un(ji,jj,jk)*fse3u(ji,jj,jk)
-                  vn_tm(ji,jj,jk)= vn_tm(ji,jj,jk)+vn(ji,jj,jk)*fse3v(ji,jj,jk)
-		  wn_tm(ji,jj,jk)= wn_tm(ji,jj,jk)+wn(ji,jj,jk)*fse3w(ji,jj,jk)
-                  avt_tm(ji,jj,jk)= avt_tm(ji,jj,jk)+avt(ji,jj,jk)*fse3t(ji,jj,jk)
-#if defined key_ldfslp
-		  wslpi_tm(ji,jj,jk)= wslpi_tm(ji,jj,jk)+wslpi(ji,jj,jk)*fse3w(ji,jj,jk)
-		  wslpj_tm(ji,jj,jk)= wslpj_tm(ji,jj,jk)+wslpj(ji,jj,jk)*fse3w(ji,jj,jk)
-		  uslp_tm(ji,jj,jk)= uslp_tm(ji,jj,jk)+uslp(ji,jj,jk)*fse3u(ji,jj,jk)
-		  vslp_tm(ji,jj,jk)= vslp_tm(ji,jj,jk)+vslp(ji,jj,jk)*fse3v(ji,jj,jk)
-#endif
-             END DO
-            END DO
-         END DO
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-         	sshn_tm(ji,jj)=sshn_tm(ji,jj)+sshn(ji,jj)
-                sshu_n_tm(ji,jj)=sshu_n_tm(ji,jj)+sshu_n(ji,jj)
-                sshv_n_tm(ji,jj)=sshv_n_tm(ji,jj)+sshv_n(ji,jj)
-	        rnf_tm(ji,jj)= rnf_tm(ji,jj)+rnf(ji,jj)
-	        h_rnf_tm(ji,jj)= h_rnf_tm(ji,jj)+h_rnf(ji,jj)
-             END DO
-         END DO
-      ELSE
-         if(lwp)write(numout,*)'trc_sub_stp at kt=',kt
-!it is time to substep
-! 1. set temporary arrays to hold physics variables
-            sshn_temp=sshn
-            sshu_n_temp=sshu_n
-            sshv_n_temp=sshv_n
-            sshf_n_temp=sshf_n
-            sshb_temp=sshb
-            sshu_b_temp=sshu_b
-            sshv_b_temp=sshv_b
-            ssha_temp=ssha
-            sshu_a_temp=sshu_a
-            sshv_a_temp=sshv_a
-            rnf_temp=rnf
-            h_rnf_temp=h_rnf
-            rotn_temp=rotn
-            hdivn_temp=hdivn
-            rotb_temp=rotb
-            hdivb_temp=hdivb
-            un_temp=un
-            vn_temp=vn
-            wn_temp=wn
-            avt_temp=avt
-#if defined key_ldfslp
-            wslpi_temp=wslpi
-            wslpj_temp=wslpj
-            uslp_temp=uslp
-            vslp_temp=vslp
-#endif
-            hu_temp=hu
-            hv_temp=hv
-            hur_temp=hur
-            hvr_temp=hvr
-
-!2. create averages and reassign variables
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  sshn_tm(ji,jj)= sshn_tm(ji,jj)+sshn(ji,jj)
-                  sshn(ji,jj)=sshn_tm(ji,jj)/(nn_dttrc+1)
-		  
-                  sshu_n_tm(ji,jj)= sshu_n_tm(ji,jj)+sshu_n(ji,jj)
-                  sshu_n(ji,jj)=sshu_n_tm(ji,jj)/(nn_dttrc+1)
-		  
-                  sshv_n_tm(ji,jj)= sshv_n_tm(ji,jj)+sshv_n(ji,jj)
-                  sshv_n(ji,jj)=sshv_n_tm(ji,jj)/(nn_dttrc+1)
-
-		  sshb(ji,jj)=sshb_hold(ji,jj)
-
-		  sshu_b(ji,jj)=sshu_b_hold(ji,jj)
-
-		  sshv_b(ji,jj)=sshv_b_hold(ji,jj)
-
-                  rnf_tm(ji,jj)= rnf_tm(ji,jj)+rnf(ji,jj)
-                  rnf(ji,jj)=rnf_tm(ji,jj)/(nn_dttrc+1)
- 
-                  h_rnf_tm(ji,jj)= h_rnf_tm(ji,jj)+h_rnf(ji,jj)
-                  h_rnf(ji,jj)=h_rnf_tm(ji,jj)/(nn_dttrc+1)
-               END DO
-            END DO
-         DO jk = 1, jpk
-             IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
-             !                                           !------------------------------!
-                 fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
-                 fse3u (:,:,jk) = fse3u_n (:,:,jk)
-                 fse3v (:,:,jk) = fse3v_n (:,:,jk)
-                 fse3f (:,:,jk) = fse3f_n (:,:,jk)
-                 fse3w (:,:,jk) = fse3w_n (:,:,jk)
-                 fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
-                 fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
-             ENDIF
-         ENDDO
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-
-                  un_tm(ji,jj,jk)= un_tm(ji,jj,jk)+un(ji,jj,jk)*e3u_temp(ji,jj,jk)
-                  un(ji,jj,jk)=(un_tm(ji,jj,jk)/(nn_dttrc+1))/fse3u(ji,jj,jk)
-		  
-                  vn_tm(ji,jj,jk)= vn_tm(ji,jj,jk)+vn(ji,jj,jk)*e3v_temp(ji,jj,jk)
-                  vn(ji,jj,jk)=(vn_tm(ji,jj,jk)/(nn_dttrc+1))/fse3v(ji,jj,jk)
-		  
-                  wn_tm(ji,jj,jk)= wn_tm(ji,jj,jk)+wn(ji,jj,jk)*e3w_temp(ji,jj,jk)
-                  wn(ji,jj,jk)=(wn_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)
-		  
-                  avt_tm(ji,jj,jk)= avt_tm(ji,jj,jk)+avt(ji,jj,jk)*e3t_temp(ji,jj,jk)
-                  avt(ji,jj,jk)=(avt_tm(ji,jj,jk)/(nn_dttrc+1))/fse3t(ji,jj,jk)	
-	  
-#if defined key_ldfslp
-                  wslpi_tm(ji,jj,jk)= wslpi_tm(ji,jj,jk)+wslpi(ji,jj,jk)*e3w_temp(ji,jj,jk)
-                  wslpi(ji,jj,jk)=(wslpi_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)
-
-                  wslpj_tm(ji,jj,jk)= wslpj_tm(ji,jj,jk)+wslpj(ji,jj,jk)*e3w_temp(ji,jj,jk)
-                  wslpj(ji,jj,jk)=(wslpj_tm(ji,jj,jk)/(nn_dttrc+1))/fse3w(ji,jj,jk)
-
-                  uslp_tm(ji,jj,jk)= uslp_tm(ji,jj,jk)+uslp(ji,jj,jk)*e3u_temp(ji,jj,jk)
-                  uslp(ji,jj,jk)=(uslp_tm(ji,jj,jk)/(nn_dttrc+1))/fse3u(ji,jj,jk)
-
-                  vslp_tm(ji,jj,jk)= vslp_tm(ji,jj,jk)+vslp(ji,jj,jk)*e3v_temp(ji,jj,jk)
-                  vslp(ji,jj,jk)=(vslp_tm(ji,jj,jk)/(nn_dttrc+1))/fse3v(ji,jj,jk)
-#endif
-	       END DO
-            END DO
-         END DO
-         CALL ssh_wzv_trc( kt )         ! after ssh & vertical velocity
-      ENDIF
-
-
-   END SUBROUTINE trc_sub_stp
-  !!-------------------------------------------------------------------
-
-  !!-------------------------------------------------------------------
-   SUBROUTINE trc_sub_stp_ini
-      !!-------------------------------------------------------------------
-      !!                     ***  ROUTINE trc_sub_stp_ini  ***
-      !!                      
-      !! ** Purpose : Initialize variables needed for sub-stepping passive tracers
-      !! 
-      !! ** Method  : 
-      !!              Compute the averages for sub-stepping
-      !!-------------------------------------------------------------------
-      INTEGER ::   jk, jn ,ji,jj   ! dummy loop indices
-      INTEGER :: ierr
-      !!-------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'trc_sub_stp_ini : initial set up of the passive tracers substepping'
-      IF(lwp) WRITE(numout,*) '~~~~~~~'
-
-      ierr =  trc_sub_alloc    ()
-      IF( lk_mpp    )   CALL mpp_sum( ierr )
-      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'top_sub_alloc : unable to allocate standard ocean arrays' )
-
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               un_tm(ji,jj,jk)= un(ji,jj,jk)*fse3u(ji,jj,jk)
-               vn_tm(ji,jj,jk)= vn(ji,jj,jk)*fse3v(ji,jj,jk)
-	       wn_tm(ji,jj,jk)= wn(ji,jj,jk)*fse3w(ji,jj,jk)
-               avt_tm(ji,jj,jk)=avt(ji,jj,jk)*fse3t(ji,jj,jk)
-#if defined key_ldfslp
-               wslpi_tm(ji,jj,jk)= wslpi(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
-               wslpj_tm(ji,jj,jk)= wslpj(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
-               uslp_tm(ji,jj,jk)= uslp(ji,jj,jk)*fse3u(ji,jj,jk) !0. 
-               vslp_tm(ji,jj,jk)= vslp(ji,jj,jk)*fse3v(ji,jj,jk) !0. 
-#endif
-           END DO
-         END DO
-      END DO
- 
-      sshn_tm = sshn !0.
-      sshu_n_tm = sshu_n !0.
-      sshv_n_tm = sshv_n !0.
-      rnf_tm= rnf !0. 
-      h_rnf_tm= h_rnf !0.
-
-   END SUBROUTINE trc_sub_stp_ini
-
-  !!-------------------------------------------------------------------
-   SUBROUTINE trc_sub_stp_reset(kt)
-      !!-------------------------------------------------------------------
-      !!                     ***  ROUTINE trc_sub_stp_reset  ***
-      !!                      
-      !! ** Purpose : Reset physics variables averaged for substepping
-      !! 
-      !! ** Method  : 
-      !!              Compute the averages for sub-stepping
-      !!-------------------------------------------------------------------
-      INTEGER ::   jk, jn ,ji,jj   ! dummy loop indices
-      INTEGER :: ierr
-      INTEGER, INTENT( in ) ::  kt  ! ocean time-step index
-      !!-------------------------------------------------------------------
-
-      IF(lwp)write(numout,*)'trc_sub_stp_reset at kt=',kt
-!restore physics variables
-        avt=avt_temp
-        un=un_temp
-        vn=vn_temp
-        wn=wn_temp
-        sshn=sshn_temp
-        sshb=sshb_temp
-        ssha=ssha_temp
-        sshu_n=sshu_n_temp
-        sshu_b=sshu_b_temp
-        sshu_a=sshu_a_temp
-        sshv_n=sshv_n_temp
-        sshv_b=sshv_b_temp
-        sshv_a=sshv_a_temp
-        sshf_n=sshf_n_temp
-#if defined key_ldfslp
-        wslpi=wslpi_temp
-        wslpj=wslpj_temp
-        uslp=uslp_temp
-        vslp=vslp_temp
-#endif
-        rnf=rnf_temp
-        h_rnf=h_rnf_temp
-        hdivn=hdivn_temp
-        rotn=rotn_temp
-        hdivb=hdivb_temp
-        rotb=rotb_temp
-        hu=hu_temp
-        hv=hv_temp
-        hur=hur_temp
-        hvr=hvr_temp
-
-       !                                           !------------------------------!
-       IF( lk_vvl ) THEN                           !  Update Now Vertical coord.  !   (only in vvl case)
-       !                                           !------------------------------!
-         DO jk = 1, jpkm1
-            fsdept(:,:,jk) = fsdept_n(:,:,jk)          ! now local depths stored in fsdep. arrays
-            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
-            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
-            !
-            fse3t (:,:,jk) = fse3t_n (:,:,jk)          ! vertical scale factors stored in fse3. arrays
-            fse3u (:,:,jk) = fse3u_n (:,:,jk)
-            fse3v (:,:,jk) = fse3v_n (:,:,jk)
-            fse3f (:,:,jk) = fse3f_n (:,:,jk)
-            fse3w (:,:,jk) = fse3w_n (:,:,jk)
-            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
-            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
-         END DO
-      ENDIF
-
-!start new averages
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  un_tm(ji,jj,jk)=un(ji,jj,jk)*fse3u(ji,jj,jk)  !0.0
-                  vn_tm(ji,jj,jk)=vn(ji,jj,jk)*fse3v(ji,jj,jk)  !0.0
-                  wn_tm(ji,jj,jk)=wn(ji,jj,jk)*fse3w(ji,jj,jk)  !0.0
-                  avt_tm(ji,jj,jk)=avt(ji,jj,jk)*fse3t(ji,jj,jk)  !0.0
-#if defined key_ldfslp
-                  wslpi_tm(ji,jj,jk)= wslpi(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
-                  wslpj_tm(ji,jj,jk)=wslpj(ji,jj,jk)*fse3w(ji,jj,jk) !0. 
-                  uslp_tm(ji,jj,jk)=uslp(ji,jj,jk)*fse3u(ji,jj,jk) !0. 
-                  vslp_tm(ji,jj,jk)=vslp(ji,jj,jk)*fse3v(ji,jj,jk) !0. 
-#endif
- 	       END DO
-            END DO
-         END DO
-        sshn_tm= sshn !0.0
-        sshu_n_tm= sshu_n !0.0
-        sshv_n_tm= sshv_n !0.0
-        rnf_tm=rnf !0. 
-        h_rnf_tm=h_rnf !0
-        sshb_hold=sshn
-        sshu_b_hold=sshu_n
-        sshv_b_hold=sshv_n
-
-   END SUBROUTINE trc_sub_stp_reset
-
-  !!-------------------------------------------------------------------
-   INTEGER FUNCTION trc_sub_alloc()
-      !!-------------------------------------------------------------------
-      !!                    *** ROUTINE trc_sub_alloc ***
-      !!-------------------------------------------------------------------
-      USE lib_mpp, ONLY: ctl_warn
-      INTEGER :: ierr
-      !!-------------------------------------------------------------------
-      !
-      ALLOCATE( un_temp(jpi,jpj,jpk)        ,  vn_temp(jpi,jpj,jpk)  ,   &
-         &      wn_temp(jpi,jpj,jpk)        ,  avt_temp(jpi,jpj,jpk) ,   &
-         &      sshn_temp(jpi,jpj)          ,  sshb_temp(jpi,jpj) ,    &
-         &      ssha_temp(jpi,jpj)          ,  sshu_a_temp(jpi,jpj),    &
-         &      sshu_n_temp(jpi,jpj)        ,  sshu_b_temp(jpi,jpj),   &
-         &      sshv_n_temp(jpi,jpj)        ,  sshv_b_temp(jpi,jpj),   &
-         &      sshv_a_temp(jpi,jpj)        ,  sshf_a_temp(jpi,jpj),   &
-         &      sshf_n_temp(jpi,jpj)        ,  sshf_b_temp(jpi,jpj),   &
-#if defined key_ldfslp
-         &      wslpi_temp(jpi,jpj,jpk) ,  wslpj_temp(jpi,jpj,jpk),   &
-         &      uslp_temp(jpi,jpj,jpk)  ,  vslp_temp(jpi,jpj,jpk),   &
-#endif
-         &      rnf_temp(jpi,jpj)           ,  h_rnf_temp(jpi,jpj) ,   &
-         &      hdivn_temp(jpi,jpj,jpk)     ,  hdivb_temp(jpi,jpj,jpk),    &
-         &      rotn_temp(jpi,jpj,jpk)      ,  rotb_temp(jpi,jpj,jpk),    &
-         &      hu_temp(jpi,jpj)            ,  hv_temp(jpi,jpj),    &
-         &      hur_temp(jpi,jpj)           ,  hvr_temp(jpi,jpj),    &
-         &      un_tm(jpi,jpj,jpk)          ,  vn_tm(jpi,jpj,jpk)  ,   &
-         &      wn_tm(jpi,jpj,jpk)          ,  avt_tm(jpi,jpj,jpk) ,   &
-         &      sshn_tm(jpi,jpj)            ,  sshb_hold(jpi,jpj) ,    &
-         &      sshu_n_tm(jpi,jpj)          ,  sshu_b_hold(jpi,jpj),   &
-         &      sshv_n_tm(jpi,jpj)          ,  sshv_b_hold(jpi,jpj),   &
-#if defined key_ldfslp
-         &      wslpi_tm(jpi,jpj,jpk)   ,  wslpj_tm(jpi,jpj,jpk),   &
-         &      uslp_tm(jpi,jpj,jpk)    ,  vslp_tm(jpi,jpj,jpk),   &
-#endif
-         &      rnf_tm(jpi,jpj)             ,  h_rnf_tm(jpi,jpj) ,   &
-         &                                    STAT=trc_sub_alloc )  
-      IF( trc_sub_alloc /= 0 )   CALL ctl_warn('trc_sub_alloc: failed to allocate arrays')
-
-      !
-   END FUNCTION trc_sub_alloc
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default key                                     NO passive tracers
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE trc_sub_stp( kt )        ! Empty routine
-      WRITE(*,*) 'trc_sub_stp: You should not have seen this print! error?', kt
-   END SUBROUTINE trc_sub_stp
-   SUBROUTINE trc_sub_stp_ini        ! Empty routine
-      WRITE(*,*) 'trc_sub_stp_ini: You should not have seen this print! error?', kt
-   END SUBROUTINE trc_sub_stp_ini
-
-#endif
-
-   !!======================================================================
-END MODULE trcsubstp
