Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/C1D_SRC/step_c1d.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/C1D_SRC/step_c1d.F90	(revision 2287)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/C1D_SRC/step_c1d.F90	(revision 2382)
@@ -13,48 +13,11 @@
    !!   stp_c1d        : NEMO system time-stepping in c1d case
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
-   USE zdf_oce         ! ocean vertical physics variables
-   USE in_out_manager  ! I/O manager
-   USE iom             !
-   USE lbclnk
-
-   USE daymod          ! calendar                         (day     routine)
-
-   USE dtatem          ! ocean temperature data           (dta_tem routine)
-   USE dtasal          ! ocean salinity    data           (dta_sal routine)
-   USE sbcmod          ! surface boundary condition       (sbc     routine)
-   USE sbcrnf          ! surface boundary condition: runoff variables
-
-   USE trcstp          ! passive tracer time-stepping      (trc_stp routine)
-
-   USE traqsr          ! solar radiation penetration      (tra_qsr routine)
-   USE trasbc          ! surface boundary condition       (tra_sbc routine)
-   !   zdfkpp          ! KPP non-local tracer fluxes      (tra_kpp routine)
-   USE trazdf          ! vertical mixing                  (tra_zdf routine)
-   USE tranxt          ! time-stepping                    (tra_nxt routine)
-   USE tranpc          ! non-penetrative convection       (tra_npc routine)
-
-   USE eosbn2          ! equation of state                (eos_bn2 routine)
-
+   USE step_oce         ! time stepping definition modules 
+#if defined key_top
+   USE trcstp           ! passive tracer time-stepping      (trc_stp routine)
+#endif
    USE dyncor_c1d      ! Coriolis term (c1d case)         (dyn_cor_1d     )
-   USE dynzdf          ! vertical diffusion               (dyn_zdf routine)
    USE dynnxt_c1d      ! time-stepping                    (dyn_nxt routine)
    USE diawri_c1d      ! write outputs                (dia_wri_c1d routine)
-
-   USE zdfbfr          ! bottom friction                  (zdf_bfr routine)
-   USE zdftke          ! TKE vertical mixing              (zdf_tke routine)
-   USE zdfkpp          ! KPP vertical mixing              (zdf_kpp routine)
-   USE zdfddm          ! double diffusion mixing          (zdf_ddm routine)
-   USE zdfevd          ! enhanced vertical diffusion      (zdf_evd routine)
-   USE zdfric          ! Richardson vertical mixing       (zdf_ric routine)
-   USE zdfmxl          ! Mixed-layer depth                (zdf_mxl routine)
-
-
-   USE diawri          ! Standard run outputs             (dia_wri routine)
-
-   USE stpctl          ! time stepping control            (stp_ctl routine)
-   USE restart         ! ocean restart                    (rst_wri routine)
-   USE prtctl          ! Print control                    (prt_ctl routine)
 
    IMPLICIT NONE
@@ -69,7 +32,6 @@
    !! NEMO/C1D 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -96,53 +58,54 @@
       !! ---------------------------------------------------------------------
 
-      indic = 1                    ! reset to no error condition
-
-      CALL day( kstp )             ! Calendar
-
-      CALL rst_opn( kstp )         ! Open the restart file
+                             indic = 0                ! reset to no error condition
+      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
+                             CALL iom_setkt( kstp )   ! say to iom that we are at time step kstp
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update data, open boundaries, surface boundary condition (including sea-ice)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
       IF( lk_dtatem  )   CALL dta_tem( kstp )         ! update 3D temperature data
       IF( lk_dtasal  )   CALL dta_sal( kstp )         ! update 3D salinity data
-
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
 
-      IF( ninist == 1 ) THEN                          ! Output the initial state and forcings
-                        CALL dia_wri_state( 'output.init', kstp )   ;   ninist = 0
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Ocean physics update                (ua, va, ta, sa used as workspace)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL bn2( tsb, rn2b )        ! before Brunt-Vaisala frequency
+                         CALL bn2( tsn, rn2  )        ! now    Brunt-Vaisala frequency
+      !  VERTICAL PHYSICS   
+                         CALL zdf_bfr( kstp )         ! bottom friction
+      !                                               ! Vertical eddy viscosity and diffusivity coefficients
+      IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz
+      IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz
+      IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
+      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz
+      IF( lk_zdfcst  )   THEN                            ! Constant Kz (reset avt, avm[uv] to the background value)
+         avt (:,:,:) = rn_avt0 * tmask(:,:,:)
+         avmu(:,:,:) = rn_avm0 * umask(:,:,:)
+         avmv(:,:,:) = rn_avm0 * vmask(:,:,:)
       ENDIF
 
+      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
+         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:)   ;   END DO
+      ENDIF
+      IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
+
+      IF( lk_zdftmx  )   CALL zdf_tmx( kstp )         ! tidal vertical mixing
+
+      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
+         &               CALL zdf_ddm( kstp )         ! double diffusive mixing
+         
+                         CALL zdf_mxl( kstp )         ! mixed layer depth
+
+                                                      ! write tke information in the restart file
+      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
+                                                      ! write gls information in the restart file
+      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Ocean physics update
+      ! diagnostics and outputs             (ua, va, ta, sa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !-----------------------------------------------------------------------
-      !  VERTICAL PHYSICS
-      !-----------------------------------------------------------------------
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
-
-                        CALL bn2( tb, sb, rn2 )             ! before Brunt-Vaisala frequency
-      
-      !                                                     ! Vertical eddy viscosity and diffusivity coefficients
-      IF( lk_zdfric )   CALL zdf_ric( kstp )                     ! Richardson number dependent Kz
-      IF( lk_zdftke )   CALL zdf_tke( kstp )                     ! TKE closure scheme for Kz
-      IF( lk_zdfkpp )   CALL zdf_kpp( kstp )                     ! KPP closure scheme for Kz
-      IF( lk_zdfcst )   THEN                                     ! Constant Kz (reset avt, avm to the background value)
-         avt (:,:,:) = avt0 * tmask(:,:,:)
-         avmu(:,:,:) = avm0 * umask(:,:,:)
-         avmv(:,:,:) = avm0 * vmask(:,:,:)
-      ENDIF
-
-      IF( ln_rnf_mouth ) THEN                                ! increase diffusivity at rivers mouths
-         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:)   ;   END DO
-      ENDIF
-      IF( ln_zdfevd )   CALL zdf_evd( kstp )                 ! enhanced vertical eddy diffusivity
-      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
-         &              CALL zdf_ddm( kstp )                 ! double diffusive mixing
-                        CALL zdf_bfr( kstp )                 ! bottom friction
-                        CALL zdf_mxl( kstp )                 ! mixed layer depth
+                         CALL dia_wri_c1d( kstp, indic )       ! ocean model: outputs
 
 #if defined key_top
@@ -150,29 +113,24 @@
       ! Passive Tracer Model
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
-                             CALL trc_stp( kstp, indic )            ! time-stepping
+                         CALL trc_stp( kstp, indic )            ! time-stepping
 #endif
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Active tracers
+      ! Active tracers                              (ua, va used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ua, va arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
                              tsa(:,:,:,:) = 0.e0                ! set tracer trends to zero
 
                              CALL tra_sbc    ( kstp )        ! surface boundary condition
       IF( ln_traqsr      )   CALL tra_qsr    ( kstp )        ! penetrative solar radiation qsr
-      IF( lk_zdfkpp )        CALL tra_kpp    ( kstp )        ! KPP non-local tracer fluxes
+      IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )        ! KPP non-local tracer fluxes
                              CALL tra_zdf    ( kstp )        ! vertical mixing
                              CALL tra_nxt    ( kstp )        ! tracer fields at next time step
       IF( ln_zdfnpc      )   CALL tra_npc    ( kstp )        ! applied non penetrative convective adjustment on (t,s)
-                             CALL eos( tsb, rhd, rhop )   ! now (swap=before) in situ density for dynhpg module
+                             CALL eos( tsb, rhd, rhop )      ! now (swap=before) in situ density for dynhpg module
+                             CALL tra_unswap                 ! udate T & S 3D arrays  (to be suppressed)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Dynamics
+      ! Dynamics                                    (ta, sa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ta, sa arrays are used as workspace in this section 
-      !-----------------------------------------------------------------------
                                ua(:,:,:) = 0.e0               ! set dynamics trends to zero
                                va(:,:,:) = 0.e0
@@ -188,12 +146,4 @@
       IF( kstp == nit000     )   CALL iom_close( numror )             ! close input  ocean restart file
       IF( lrst_oce           )   CALL rst_write  ( kstp )             ! write output ocean restart file
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! diagnostics and outputs
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
-
-      IF( nstop == 0 )           CALL dia_wri_c1d( kstp, indic )       ! ocean model: outputs
       !
    END SUBROUTINE stp_c1d
