Index: trunk/NEMO/OPA_SRC/IOM/in_out_manager.F90
===================================================================
--- trunk/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 1480)
+++ trunk/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 1481)
@@ -44,7 +44,4 @@
    !!----------------------------------------------------------------------
    INTEGER            ::   nitrst                 !: time step at which restart file should be written
-#if defined key_zdftke2
-   INTEGER            ::   nitrst_tke2            !: time step at which restart file should be written
-#endif
    !!----------------------------------------------------------------------
    !!                    output monitoring
Index: trunk/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 1480)
+++ trunk/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 1481)
@@ -232,16 +232,4 @@
          !
          !                                           ! Diagnostics and outputs 
-#if defined key_zdftke2
-         IF( kt .NE. nitend+1 ) THEN
-            IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   &
-               &             CALL lim_dia 
-                             CALL lim_wri( 1  )              ! Ice outputs 
-            IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file 
-                             CALL lim_var_glo2eqv            ! ???
-            !
-            IF( ln_nicep )   CALL lim_ctl               ! alerts in case of model crash
-            !
-         ENDIF
-#else
          IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   &
             &             CALL lim_dia 
@@ -252,5 +240,4 @@
          IF( ln_nicep )   CALL lim_ctl               ! alerts in case of model crash
          !
-#endif
       ENDIF                                          ! End sea-ice time step only
 
Index: trunk/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 1480)
+++ trunk/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 1481)
@@ -171,11 +171,5 @@
          ! ---------------- !
 
-#if defined key_zdftke2
-         IF ( kt .NE. nitend+1 ) THEN
                                         CALL lim_rst_opn_2  ( kt )      ! Open Ice restart file
-         ENDIF
-#else
-                                        CALL lim_rst_opn_2  ( kt )      ! Open Ice restart file
-#endif
          IF( .NOT. lk_c1d ) THEN                                        ! Ice dynamics & transport (not in 1D case)
                                         CALL lim_dyn_2      ( kt )           ! Ice dynamics    ( rheology/dynamics )
@@ -192,13 +186,5 @@
                                         CALL lim_thd_2      ( kt )      ! Ice thermodynamics 
                                         CALL lim_sbc_2      ( kt )      ! Ice/Ocean Mass & Heat fluxes 
-#if defined key_zdftke2
-         IF( kt .NE. nitend+1 ) THEN
-            IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   &
-               &                        CALL lim_dia_2      ( kt )      ! Ice Diagnostics 
-                                        CALL lim_wri_2      ( kt )      ! Ice outputs 
-            IF( lrst_ice )              CALL lim_rst_write_2( kt )      ! Ice restart file 
-            !
-         ENDIF
-#else
+
          IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   &
             &                           CALL lim_dia_2      ( kt )      ! Ice Diagnostics
@@ -206,5 +192,4 @@
          IF( lrst_ice )                 CALL lim_rst_write_2( kt )      ! Ice restart file
          !
-#endif
       ENDIF
       !
Index: trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90	(revision 1480)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90	(revision 1481)
@@ -50,9 +50,10 @@
 
    PUBLIC   zdf_tke2    ! routine called in step module
+   PUBLIC   tke2_rst    ! routine called in step module
 
    LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke2 = .TRUE.  !: TKE vertical mixing flag
    REAL(wp), PUBLIC                         ::   eboost              !: multiplicative coeff of the shear product.
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   en                  !: now turbulent kinetic energy
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   avmt                !: now mixing coefficient of diffusion for TKE
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   avm                 !: now mixing coefficient of diffusion for TKE
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   dissl               !: now mixing lenght of dissipation
 
@@ -166,17 +167,12 @@
       INTEGER, INTENT(in) ::   kt ! ocean time step
 
-      IF( kt == nit000  ) THEN
-         CALL zdf_tke2_init        ! Initialization (first time-step only)
-         IF( .NOT. ln_rstart )  CALL tke_tke( kt )
-         CALL tke_mlmc
-      ELSE
-         CALL tke_tke( kt )
-         IF( kt /= nitend+1 )   CALL tke_mlmc
-         IF( lrst_oce ) CALL tke2_rst( kt, 'WRITE' )      ! write en in restart file 
-      ENDIF
+      IF( kt == nit000 )   CALL zdf_tke2_init             ! Initialization (first time-step only)
+      !
+      CALL tke_tke                                        ! now tke (en)
+      CALL tke_mlmc                                       ! now avmu, avmv, avt
 
    END SUBROUTINE zdf_tke2
 
-   SUBROUTINE tke_tke( kt )
+   SUBROUTINE tke_tke
       !!----------------------------------------------------------------------
       !! Now Turbulent kinetic energy (output in en)
@@ -185,5 +181,5 @@
       !! Resolution of a tridiagonal linear system by a "methode de chasse"
       !! computation from level 2 to jpkm1  (e(1) computed and e(jpk)=0 ).
-      !! avmt represents the turbulent coefficient of the TKE
+      !! avm represents the turbulent coefficient of the TKE
       !!----------------------------------------------------------------------
       USE oce,     zwd     =>   ua   ! use ua as workspace
@@ -191,6 +187,4 @@
       USE oce,     zdiag2  =>   ta   ! use ta as workspace
       USE oce,     ztkelc  =>   sa   ! use sa as workspace
-      !
-      INTEGER, INTENT(in) ::   kt ! ocean time step
       !
       INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
@@ -209,5 +203,19 @@
       zfact1 = -.5 * rdt * rn_efave
       zfact2 = 1.5 * rdt * rn_ediss
-      zfact3 = 0.5 * rdt * rn_ediss
+      zfact3 = 0.5       * rn_ediss
+
+      ! Surface boundary condition on tke
+      ! -------------------------------------------------
+      ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0)
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
+            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
+            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
+            en(ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1)
+         END DO
+      END DO
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -217,7 +225,7 @@
          !
          ! Computation of total energy produce by LC : cumulative sum over jk
-         zpelc(:,:,1) =  MAX( rn2(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)
+         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)
          DO jk = 2, jpk
-            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
+            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
          END DO
          !
@@ -246,5 +254,5 @@
 # endif
          !
-         ! TKE Langmuir circulation source term
+         ! TKE Langmuir circulation source term added to en
 !CDIR NOVERRCHK
          DO jk = 2, jpkm1
@@ -259,5 +267,5 @@
                   zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
                   ! TKE Langmuir circulation source term
-                  ztkelc(ji,jj,jk) = ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
+                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) ) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -266,26 +274,15 @@
       ENDIF
 
-      ! Surface boundary condition on tke
-      ! -------------------------------------------------
-      ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0)
-!CDIR NOVERRCHK
-      DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
-            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
-            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
-            en (ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1)
-         END DO
-      END DO
-
-      ! Preliminary computing
+
+      ! Shear production at uw- and vw-points (energy conserving form)
       DO jk = 2, jpkm1
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
                avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) * &
-                  &                              ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )
+                  &                              ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) / & 
+                  &                              ( fse3uw_n(ji,jj,jk) * fse3uw_b(ji,jj,jk) )
                avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) * &
-                  &                              ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )
+                  &                              ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) / &
+                  &                              ( fse3vw_n(ji,jj,jk) * fse3vw_b(ji,jj,jk) )
             ENDDO
          ENDDO
@@ -304,18 +301,14 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               !                                                             ! shear prod. - stratif. destruction
-               zcoef = 0.5 / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )
-               ! shear
-               zdku = zcoef * ( avmu(ji-1,jj  ,jk) + avmu(ji,jj,jk) )        ! shear
-               zdkv = zcoef * ( avmv(ji  ,jj-1,jk) + avmv(ji,jj,jk) )
-               zesh2 =  eboost * ( zdku + zdkv ) - avt(ji,jj,jk) * rn2(ji,jj,jk)        ! coefficient (zesh2)
-               !
+               !                                                             ! shear prod. at w-point weightened by mask
+               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
                !                                                             ! Matrix
                zcof = zfact1 * tmask(ji,jj,jk)
                !                                                                  ! lower diagonal
-               zdiag2(ji,jj,jk) = zcof * ( avmt (ji,jj,jk  ) + avmt (ji,jj,jk-1) )   &
+               zdiag2(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
                   &                    / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
                !                                                                  ! upper diagonal
-               zdiag1(ji,jj,jk) = zcof * ( avmt (ji,jj,jk+1) + avmt (ji,jj,jk  ) )   &
+               zdiag1(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
                   &                    / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
                !                                                                  ! diagonal
@@ -323,12 +316,6 @@
                !
                !                                                             ! right hand side in en
-               IF( .NOT. ln_lc ) THEN                                             ! No Langmuir cells
-                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * dissl (ji,jj,jk) * en(ji,jj,jk) * tmask(ji,jj,jk) &
-                     &                        +   rdt  * zesh2 * tmask(ji,jj,jk)
-               ELSE                                                               ! Langmuir cell source term
-                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * dissl (ji,jj,jk) * en(ji,jj,jk) * tmask(ji,jj,jk) &
-                     &                        +   rdt  * zesh2 * tmask(ji,jj,jk)                           &
-                     &                        +   rdt  * ztkelc(ji,jj,jk) * tmask(ji,jj,jk)
-               ENDIF
+               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
+                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) * tmask(ji,jj,jk)
             END DO
          END DO
@@ -434,11 +421,10 @@
       USE oce,     zmxlm  =>   va   ! use va as workspace
       USE oce,     zmxld  =>   ta   ! use ta as workspace
-      USE oce,     visco  =>   sa   ! use sa as workspace
       !
       INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
-      REAL(wp) ::   zbbrau, zrn2, zraug             ! temporary scalars
-      REAL(wp) ::   zfact1, ztx2, zdku              !    -         -
-      REAL(wp) ::   zfact2, zty2, zdkv              !    -         -
-      REAL(wp) ::   zfact3, zcoef, zav              !    -         -
+      REAL(wp) ::   zrn2, zraug                     ! temporary scalars
+      REAL(wp) ::   ztx2, zdku                      !    -         -
+      REAL(wp) ::   zty2, zdkv                      !    -         -
+      REAL(wp) ::   zcoef, zav                      !    -         -
       REAL(wp) ::   zpdl, zri, zsqen                !    -         -
       REAL(wp) ::   zemxl, zemlm, zemlp             !    -         -
@@ -569,8 +555,8 @@
       ! II  Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt !
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<!
-      ! Surface boundary condition on avt and avmt : jk = 1
+      ! Surface boundary condition on avt and avm : jk = 1
       ! -------------------------------------------------
-      ! avt (1) = avmb(1) and avt (jpk) = 0.
-      ! avmt(1) = avmb(1) and avmt(jpk) = 0.
+      ! avt(1) = avmb(1) and avt(jpk) = 0.
+      ! avm(1) = avmb(1) and avm(jpk) = 0.
       !
 !CDIR NOVERRCHK
@@ -582,7 +568,6 @@
                zsqen = SQRT( en(ji,jj,jk) )
                zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
-               avmt (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
-               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
-               visco(ji,jj,jk) = avmt (ji,jj,jk) 
+               avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk)
+               avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
                dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
             END DO
@@ -591,23 +576,16 @@
 
       ! Lateral boundary conditions (sign unchanged)
-      CALL lbc_lnk( avmt, 'W', 1. )   ;   CALL lbc_lnk( visco, 'W', 1. )
+      CALL lbc_lnk( avm, 'W', 1. )
 
       DO jk = 2, jpkm1                                 ! only vertical eddy viscosity
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               avmu(ji,jj,jk) = ( visco(ji,jj,jk) + visco(ji+1,jj  ,jk) ) * umask(ji,jj,jk)   &
-               &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
-               avmv(ji,jj,jk) = ( visco(ji,jj,jk) + visco(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)   &
-               &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
+               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
+               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
       CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
-      !
-      !
-      DO jk = 2, jpkm1                                  ! Minimum value on the eddy viscosity
-         avmu(:,:,jk) = MAX( avmu(:,:,jk), avmb(jk) ) * umask(:,:,jk)
-         avmv(:,:,jk) = MAX( avmv(:,:,jk), avmb(jk) ) * vmask(:,:,jk)
-      END DO
+
 
       ! Prandtl number
@@ -623,5 +601,5 @@
                   zdkv = avmv(ji,jj-1,jk) * (vn(ji,jj-1,jk-1)-vn(ji,jj-1,jk)) * (vb(ji,jj-1,jk-1)-vb(ji,jj-1,jk)) + &
                     &    avmv(ji,jj  ,jk) * (vn(ji,jj  ,jk-1)-vn(ji,jj  ,jk)) * (vb(ji,jj  ,jk-1)-vb(ji,jj  ,jk))
-                  zri = MAX( rn2(ji,jj,jk), 0. ) / (zcoef * (zdku + zdkv + 1.e-20) / avmt(ji,jj,jk)) ! local Richardson number
+                  zri = MAX( rn2b(ji,jj,jk), 0. ) * avm(ji,jj,jk)/ (zcoef * (zdku + zdkv + 1.e-20) ) ! local Richardson number
 # if defined key_cfg_1d
                   e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
@@ -770,7 +748,7 @@
       DO jk = 1, jpk
          avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
+         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
          avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
          avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
-         avmt(:,:,jk) = avmb(jk) * tmask(:,:,jk)
       END DO
       dissl(:,:,:) = 1.e-12
@@ -796,27 +774,55 @@
      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
      !
-     INTEGER ::   jit   ! dummy loop indices
+     INTEGER ::   jit, jk   ! dummy loop indices
+     INTEGER ::   id1, id2, id3, id4, id5, id6
      !!----------------------------------------------------------------------
      !
-     IF( TRIM(cdrw) == 'READ' ) THEN
-        IF( ln_rstart ) THEN
-           IF( iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 .AND. .NOT.(ln_rstke) ) THEN 
+     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+        !                                   ! ---------------
+        IF( ln_rstart ) THEN                   !* Read the restart file
+           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
+           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
+           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
+           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
+           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
+           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
+           !
+           IF( id1 > 0 ) THEN                       ! 'en' exists
               CALL iom_get( numror, jpdom_autoglo, 'en', en )
-           ELSE
-              IF( lwp .AND. iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 )   &
-                 &                       WRITE(numout,*) ' ===>>>> : previous run without tke scheme'
-              IF( lwp .AND. ln_rstke )   WRITE(numout,*) ' ===>>>> : We do not use en from the restart file'
-              IF( lwp                )   WRITE(numout,*) ' ===>>>> : en is set by iterative loop'
+              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
+                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
+                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
+                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
+                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
+                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
+              ELSE                                       ! one at least array is missing
+                 CALL tke_mlmc                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
+              ENDIF
+           ELSE                                     ! No TKE array found: initialisation
+              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
               en (:,:,:) = rn_emin * tmask(:,:,:)
-              DO jit = 2, nn_itke + 1
-                 CALL zdf_tke2( jit )
-              END DO
+              CALL tke_mlmc                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
+              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke2( jit )   ;   END DO
            ENDIF
-        ELSE
-           en(:,:,:) = rn_emin * tmask(:,:,:)      ! no restart: en set to its minimum value
+        ELSE                                   !* Start from rest
+           en(:,:,:) = rn_emin * tmask(:,:,:)
+           DO jk = 1, jpk                           ! set the Kz to the background value
+              avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
+              avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
+              avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
+              avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
+           END DO
         ENDIF
-     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
+        !
+     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+        !                                   ! -------------------
         IF(lwp) WRITE(numout,*) '---- tke2-rst ----'
-        CALL iom_rstput( kt, nitrst_tke2, numrow, 'en', en )
+        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
+        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt    )
+        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm    )
+        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu   )
+        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv   )
+        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
+        !
      ENDIF
      !
Index: trunk/NEMO/OPA_SRC/restart.F90
===================================================================
--- trunk/NEMO/OPA_SRC/restart.F90	(revision 1480)
+++ trunk/NEMO/OPA_SRC/restart.F90	(revision 1481)
@@ -23,8 +23,7 @@
    USE zpshde          ! partial step: hor. derivative (zps_hde routine)
    USE eosbn2          ! equation of state            (eos bn2 routine)
+   USE zdfddm          ! double diffusion mixing 
+   USE zdfmxl          ! mixed layer depth
    USE trdmld_oce      ! ocean active mixed layer tracers trends variables
-#if defined key_zdftke2
-   USE zdf_oce
-#endif
 
    IMPLICIT NONE
@@ -67,7 +66,4 @@
          lrst_oce = .FALSE.   
          nitrst = nitend
-#if defined key_zdftke2
-         nitrst_tke2 = nitrst + 1
-#endif
       ENDIF
       IF( MOD( kt - 1, nstock ) == 0 ) THEN   
@@ -76,7 +72,4 @@
          IF( nitrst > nitend )   nitrst = nitend   ! make sure we write a restart at the end of the run
       ENDIF
-#if defined key_zdftke2
-      IF ( nitrst_tke2 .NE. kt ) nitrst_tke2 = nitrst + 1
-#endif
       ! to get better performances with NetCDF format:
       ! we open and define the ocean restart file one time step before writing the data (-> at nitrst - 1)
@@ -119,10 +112,4 @@
       !!----------------------------------------------------------------------
 
-#if defined key_zdftke2
-      IF( kt == nitrst_tke2 ) THEN
-         CALL iom_close( numrow )     ! close the restart file (only at last time step)
-         IF( .NOT. lk_trdmld )   lrst_oce = .FALSE.
-      ELSE
-#endif
       !                                                                     ! the begining of the run [s]
       CALL iom_rstput( kt, nitrst, numrow, 'rdt'    , rdt               )   ! dynamics time step
@@ -138,37 +125,17 @@
       CALL iom_rstput( kt, nitrst, numrow, 'un'     , un      )
       CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vn      )
-      IF( lk_vvl ) CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn      )
       CALL iom_rstput( kt, nitrst, numrow, 'tn'     , tn      )
       CALL iom_rstput( kt, nitrst, numrow, 'sn'     , sn      )
       CALL iom_rstput( kt, nitrst, numrow, 'rotn'   , rotn    )
       CALL iom_rstput( kt, nitrst, numrow, 'hdivn'  , hdivn   )
-      CALL iom_rstput( kt, nitrst, numrow, 'rn2' , rn2  )
-
-#if defined key_zdftke2
-         CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop    )
-#endif
-      IF( nn_dynhpg_rst == 1 .OR. lk_vvl ) THEN
-         CALL iom_rstput( kt, nitrst, numrow, 'rhd' , rhd  )
-#if defined key_zdftke2
-         CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt  )
-      ENDIF
-#else
-         CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop )
-#endif
-         IF( ln_zps ) THEN
-            CALL iom_rstput( kt, nitrst, numrow, 'gtu' , gtu )
-            CALL iom_rstput( kt, nitrst, numrow, 'gsu' , gsu )
-            CALL iom_rstput( kt, nitrst, numrow, 'gru' , gru )
-            CALL iom_rstput( kt, nitrst, numrow, 'gtv' , gtv )
-            CALL iom_rstput( kt, nitrst, numrow, 'gsv' , gsv )
-            CALL iom_rstput( kt, nitrst, numrow, 'grv' , grv )
-         ENDIF
-      ENDIF
-#if ! defined key_zdftke2
+
+      CALL iom_rstput( kt, nitrst, numrow, 'rhop'  , rhop   )
+
+      IF( lk_zdfddm )   CALL iom_rstput( kt, nitrst, numrow, 'rrau' , rrau  )
+
       IF( kt == nitrst ) THEN
          CALL iom_close( numrow )     ! close the restart file (only at last time step)
          IF( .NOT. lk_trdmld )   lrst_oce = .FALSE.
       ENDIF
-#endif
       !
    END SUBROUTINE rst_write
@@ -232,5 +199,4 @@
       CALL iom_get( numror, jpdom_autoglo, 'un'   , un    )        ! now    i-component velocity
       CALL iom_get( numror, jpdom_autoglo, 'vn'   , vn    )        ! now    j-component velocity
-      IF( lk_vvl ) CALL iom_get( numror, jpdom_autoglo, 'wn'   , wn    )        ! now    k-component velocity
       CALL iom_get( numror, jpdom_autoglo, 'tn'   , tn    )        ! now    temperature
       CALL iom_get( numror, jpdom_autoglo, 'sn'   , sn    )        ! now    salinity
@@ -238,4 +204,14 @@
       CALL iom_get( numror, jpdom_autoglo, 'hdivn', hdivn )        ! now    horizontal divergence
 
+      CALL iom_get( numror, jpdom_autoglo, 'rhop', rhop )          ! now    mixed layer depth
+
+      IF( lk_zdfddm ) THEN
+         IF( iom_varid( numror, 'rrau', ldstop = .FALSE. ) > 0 ) THEN
+            CALL iom_get( numror, jpdom_autoglo, 'rrau' , rrau  )
+         ELSE
+            CALL eos_init                   ! read equation state type neos parameter
+            CALL eos( tb, sb, rhd, rhop )   ! compute rrau
+         ENDIF
+      ENDIF 
 
       IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
@@ -247,49 +223,4 @@
          hdivb(:,:,:) = hdivn(:,:,:)
       ENDIF
-      CALL eos_init            ! usefull to get the equation state type neos parameter
-      IF( iom_varid( numror, 'rn2', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'rn2' , rn2  )
-      ELSE
-         CALL bn2( tb, sb, rn2 )  ! before Brunt-Vaisala frequency   
-      ENDIF
-
-#if defined key_zdftke2
-      CALL eos( tb, sb, rhd, rhop )        ! before potential and in situ densities
-      IF( iom_varid( numror, 'rhd', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'rhd' , rhd  )
-      ENDIF
-      IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'rhop', rhop )
-      ENDIF
-#else
-      IF( iom_varid( numror, 'rhd', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'rhd' , rhd  )
-         CALL iom_get( numror, jpdom_autoglo, 'rhop', rhop )
-      ELSE
-         CALL eos( tb, sb, rhd, rhop )        ! before potential and in situ densities
-      ENDIF
-#endif
-#if defined key_zdftke2
-      IF( iom_varid( numror, 'avt', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'avt' , avt  )
-      ELSE
-         IF ( ln_dynhpg_imp ) avt (:,:,:) = 1.2e-5 * tmask(:,:,:)
-      ENDIF
-#endif
-
-      IF( ln_zps .AND. .NOT. lk_c1d ) THEN
-         IF( iom_varid( numror, 'gtu', ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numror, jpdom_autoglo, 'gtu' , gtu )
-            CALL iom_get( numror, jpdom_autoglo, 'gsu' , gsu )
-            CALL iom_get( numror, jpdom_autoglo, 'gru' , gru )
-            CALL iom_get( numror, jpdom_autoglo, 'gtv' , gtv )
-            CALL iom_get( numror, jpdom_autoglo, 'gsv' , gsv )
-            CALL iom_get( numror, jpdom_autoglo, 'grv' , grv )
-         ELSE
-            CALL zps_hde( nit000, tb , sb , rhd,   &  ! Partial steps: before Horizontal DErivative
-               &                  gtu, gsu, gru,   &  ! of t, s, rd at the bottom ocean level
-               &                  gtv, gsv, grv )
-         ENDIF
-      ENDIF
       !
    END SUBROUTINE rst_read
Index: trunk/NEMO/OPA_SRC/step.F90
===================================================================
--- trunk/NEMO/OPA_SRC/step.F90	(revision 1480)
+++ trunk/NEMO/OPA_SRC/step.F90	(revision 1481)
@@ -207,24 +207,8 @@
       ! N.B. ua, va, ta, sa arrays are used as workspace in this section
       !-----------------------------------------------------------------------
-      IF( neuler == 0 .AND. kstp == nit000 ) THEN
-                        CALL bn2( tn, sn, rn2 )              ! now    Brunt-Vaisala frequency
-                        rn2b(:,:,:) = rn2(:,:,:)
-      ELSE
-                        rn2b(:,:,:) = rn2(:,:,:)             ! before Brunt-Vaisala frequency  
-                        CALL bn2( tn, sn, rn2 )              ! now    Brunt-Vaisala frequency
-      ENDIF
-
-#if defined key_zdftke2
-      IF ( ln_dynhpg_imp ) THEN
-      !-----------------------------------------------------------------------
-      !  LATERAL PHYSICS
-      !-----------------------------------------------------------------------
-                               CALL zdf_mxl( kstp )                 ! mixed layer depth
-         IF( lk_ldfslp     )   CALL ldf_slp( kstp, rhd, rn2b )      ! before slope of the lateral mixing
-#  if defined key_traldf_c2d
-         IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )                 ! eddy induced velocity coefficient
-#  endif
-      ENDIF
-#endif
+
+                       CALL bn2( tb, sb, rn2b )              ! before Brunt-Vaisala frequency
+                       CALL bn2( tn, sn, rn2  )              ! now    Brunt-Vaisala frequency
+
       !-----------------------------------------------------------------------
       !  VERTICAL PHYSICS
@@ -260,18 +244,20 @@
                         CALL zdf_mxl( kstp )                 ! mixed layer depth
 
-#if defined key_zdftke2
-      IF( .NOT. ln_dynhpg_imp ) THEN
-                        CALL eos( tb, sb, rhd, rhop )        ! now (swap=before) in situ density for dynhpg module
-#endif
+      IF( lrst_oce .AND. lk_zdftke2 )   &                    ! write tke2 information in the restart file
+         &               CALL tke2_rst( kstp, 'WRITE' )
+
       !-----------------------------------------------------------------------
       !  LATERAL PHYSICS
       !-----------------------------------------------------------------------
-      IF( lk_ldfslp     )   CALL ldf_slp( kstp, rhd, rn2b )      ! before slope of the lateral mixing
+      IF( lk_ldfslp )   THEN 
+                           CALL eos( tb, sb, rhd, rhop )          ! before in situ density
+         IF( ln_zps )      CALL zps_hde( kstp, tb, sb, rhd,  &    ! Partial steps: before horizontal gradient
+            &                                gtu, gsu, gru,  &    ! of t, s, rd at the last ocean level
+            &                                gtv, gsv, grv )
+                               CALL ldf_slp( kstp, rhd, rn2b )      ! before slope of the lateral mixing
+      ENDIF
 #if defined key_traldf_c2d
       IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )                 ! eddy induced velocity coefficient
 #  endif
-#if defined key_zdftke2
-      ENDIF
-#endif
 
 #if defined key_top
@@ -307,20 +293,4 @@
                              CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields
 
-#if ! defined key_zdftke2
-      IF( ln_zdfnpc      )   CALL tra_npc    ( kstp )       ! update after fields by non-penetrative convection
-                             CALL tra_nxt    ( kstp )       ! tracer fields at next time step
-
-      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg
-                               CALL eos( ta, sa, rhd, rhop )          ! Time-filtered in situ density used in dynhpg module
-         IF( ln_zps    )       CALL zps_hde( kstp, ta, sa, rhd,   &   ! Partial steps: time filtered hor. gradient
-            &                                     gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
-            &                                     gtv, gsv, grv )
-      ELSE                                                  ! centered hpg (default case)
-                               CALL eos( tb, sb, rhd, rhop )          ! now (swap=before) in situ density for dynhpg module
-         IF( ln_zps    )       CALL zps_hde( kstp, tb, sb, rhd,   &   ! Partial steps: now horizontal gradient
-            &                                     gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
-            &                                     gtv, gsv, grv )
-      ENDIF
-#else
       IF( .NOT. ln_dynhpg_imp  ) THEN                       ! centered hpg (default case)
                                CALL eos( tn, sn, rhd, rhop )          ! now (swap=before) in situ density for dynhpg module
@@ -333,10 +303,9 @@
       IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg
                                CALL eos( ta, sa, rhd, rhop )          ! Time-filtered in situ density used in dynhpg module
-         IF( lk_ldfslp     )   CALL bn2( ta, sa, rn2 )                ! Time-filtered Brunt-Vaisala frequency
          IF( ln_zps    )       CALL zps_hde( kstp, ta, sa, rhd,   &   ! Partial steps: time filtered hor. gradient
             &                                     gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
             &                                     gtv, gsv, grv )
       ENDIF
-#endif
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Dynamics
@@ -400,24 +369,5 @@
 
       IF( lk_cpl )   CALL sbc_cpl_snd( kstp )                 ! coupled mode : field exchanges
-
-#if defined key_zdftke2
-      IF( ( kstp == nitend ).AND.( lrst_oce ) )  THEN
-
-         CALL day( kstp+1 )             ! Calendar
-         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-         ! Update data, open boundaries, surface boundary condition (including sea-ice)
-         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                            CALL sbc    ( kstp+1 )         ! Sea Boundary Condition (including sea-ice)
-         !-----------------------------------------------------------------------
-         !  VERTICAL PHYSICS
-         !-----------------------------------------------------------------------
-         ! N.B. ua, va, ta, sa arrays are used as workspace in this section
-         !-----------------------------------------------------------------------
-                            CALL bn2( tn, sn, rn2 )            ! now Brunt-Vaisala frequency
-         !                                                     ! Vertical eddy viscosity and diffusivity coefficients
-         IF( lk_zdftke2 )   CALL zdf_tke2 ( kstp+1 )           ! TKE2 closure scheme for Kz
-                            CALL rst_write( kstp+1 )           ! close the restart file
-      ENDIF
-#endif
+      !
       !
    END SUBROUTINE stp
