Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 4666)
@@ -40,4 +40,5 @@
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   bfrua, bfrva   !: Bottom friction coefficients set in zdfbfr
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   tfrua, tfrva   !: top friction coefficients set in zdfbfr
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avmu , avmv    !: vertical viscosity coef at uw- & vw-pts       [m2/s]
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm  , avt     !: vertical viscosity & diffusivity coef at w-pt [m2/s]
@@ -57,4 +58,5 @@
       ALLOCATE(avmb(jpk) , bfrua(jpi,jpj) ,                         &
          &     avtb(jpk) , bfrva(jpi,jpj) , avtb_2d(jpi,jpj) ,      &
+         &     tfrua(jpi, jpj), tfrva(jpi, jpj)              ,      &
          &     avmu(jpi,jpj,jpk), avm(jpi,jpj,jpk)           ,      &
          &     avmv(jpi,jpj,jpk), avt(jpi,jpj,jpk)           , STAT = zdf_oce_alloc )
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 4666)
@@ -41,8 +41,16 @@
    REAL(wp), PUBLIC ::   rn_bfrien    ! local factor to enhance coefficient bfri (PUBLIC for TAM)
    LOGICAL , PUBLIC ::   ln_bfr2d     ! logical switch for 2D enhancement (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfri1     ! top drag coefficient (linear case)  (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfri2     ! top drag coefficient (non linear case) (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfri2_max ! Maximum top drag coefficient (non linear case and ln_loglayer=T) (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfeb2     ! background top turbulent kinetic energy  [m2/s2] (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfrien    ! local factor to enhance coefficient tfri (PUBLIC for TAM)
+   LOGICAL , PUBLIC ::   ln_tfr2d     ! logical switch for 2D enhancement (PUBLIC for TAM)
+
    LOGICAL , PUBLIC ::   ln_loglayer  ! switch for log layer bfr coeff. (PUBLIC for TAM)
    REAL(wp), PUBLIC ::   rn_bfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM)
+   REAL(wp), PUBLIC ::   rn_tfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM)
    LOGICAL , PUBLIC ::   ln_bfrimp    ! logical switch for implicit bottom friction
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d            ! 2D bottom drag coefficient (PUBLIC for TAM)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d   ! 2D bottom/top drag coefficient (PUBLIC for TAM)
 
    !! * Substitutions
@@ -60,5 +68,5 @@
       !!                ***  FUNCTION zdf_bfr_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( bfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )
+      ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )
       !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_bfr_alloc )
@@ -88,5 +96,5 @@
       INTEGER  ::   ikbt, ikbu, ikbv             ! local integers
       REAL(wp) ::   zvu, zuv, zecu, zecv, ztmp   ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zbfrt
+      REAL(wp), POINTER, DIMENSION(:,:) ::  zbfrt, ztfrt
       !!----------------------------------------------------------------------
       !
@@ -101,5 +109,5 @@
       IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only
          !
-         CALL wrk_alloc( jpi, jpj, zbfrt )
+         CALL wrk_alloc( jpi, jpj, zbfrt, ztfrt )
 
          IF ( ln_loglayer.AND.lk_vvl ) THEN ! "log layer" bottom friction coefficient
@@ -120,4 +128,11 @@
                   zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp)
                   zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max)
+! (ISF)
+                  ikbt = mikt(ji,jj)
+! JC: possible WAD implementation should modify line below if layers vanish
+                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
+                  ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp)
+                  ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max)
+
                END DO
             END DO
@@ -125,4 +140,5 @@
          ELSE
             zbfrt(:,:) = bfrcoef2d(:,:)
+            ztfrt(:,:) = tfrcoef2d(:,:)
          ENDIF
 
@@ -150,4 +166,21 @@
                bfrua(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) * zecu
                bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) * zecv
+               !
+               ! (ISF) ========================================================================
+               ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points 
+               ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points)
+               !
+               zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     &
+                  &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  )
+               zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     &
+                  &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  )
+               !
+               zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 )
+               zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 )
+               !
+               tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu 
+               tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv
+               ! (ISF) END ====================================================================
+
             END DO
          END DO
@@ -158,5 +191,5 @@
          IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        &
             &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 )
-         CALL wrk_dealloc( jpi,jpj,zbfrt )
+         CALL wrk_dealloc( jpi,jpj,zbfrt, ztfrt )
       ENDIF
       !
@@ -183,8 +216,10 @@
       INTEGER   ::   ios
       REAL(wp)  ::   zminbfr, zmaxbfr   ! temporary scalars
+      REAL(wp)  ::   zmintfr, zmaxtfr   ! temporary scalars
       REAL(wp)  ::   ztmp, zfru, zfrv   !    -         -
       !!
       NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfri2_max, rn_bfeb2, rn_bfrz0, ln_bfr2d, &
-                    &  rn_bfrien, ln_bfrimp, ln_loglayer
+                    &          rn_tfri1, rn_tfri2, rn_tfri2_max, rn_tfeb2, rn_tfrz0, ln_tfr2d, &
+                    &  rn_bfrien, rn_tfrien, ln_bfrimp, ln_loglayer
       !!----------------------------------------------------------------------
       !
@@ -215,11 +250,18 @@
          bfrua(:,:) = 0.e0
          bfrva(:,:) = 0.e0
+         tfrua(:,:) = 0.e0
+         tfrva(:,:) = 0.e0
          !
       CASE( 1 )
          IF(lwp) WRITE(numout,*) '      linear botton friction'
-         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri1  = ', rn_bfri1
+         IF(lwp) WRITE(numout,*) '      bottom friction coef.   rn_bfri1  = ', rn_bfri1
          IF( ln_bfr2d ) THEN
             IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d
             IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien
+         ENDIF
+         IF(lwp) WRITE(numout,*) '      top    friction coef.   rn_bfri1  = ', rn_bfri1
+         IF( ln_tfr2d ) THEN
+            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d
+            IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien
          ENDIF
          !
@@ -236,4 +278,17 @@
          bfrua(:,:) = - bfrcoef2d(:,:)
          bfrva(:,:) = - bfrcoef2d(:,:)
+         !
+         IF(ln_tfr2d) THEN
+            ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
+            CALL iom_open('tfr_coef.nc',inum)
+            CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array
+            CALL iom_close(inum)
+            tfrcoef2d(:,:) = rn_tfri1 * ( 1 + rn_tfrien * tfrcoef2d(:,:) )
+         ELSE
+            tfrcoef2d(:,:) = rn_tfri1  ! initialize tfrcoef2d to the namelist variable
+         ENDIF
+         !
+         tfrua(:,:) = - tfrcoef2d(:,:)
+         tfrva(:,:) = - tfrcoef2d(:,:)
          !
       CASE( 2 )
@@ -252,4 +307,18 @@
             IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien
          ENDIF
+         IF(lwp) WRITE(numout,*) '      quadratic top    friction'
+         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri2  = ', rn_tfri2
+         IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_tfri2_max  = ', rn_tfri2_max
+         IF(lwp) WRITE(numout,*) '      background tke   rn_tfeb2  = ', rn_tfeb2
+         IF(lwp) WRITE(numout,*) '      log formulation   ln_tfr2d = ', ln_loglayer
+         IF(lwp) WRITE(numout,*) '      bottom roughness  rn_tfrz0 [m] = ', rn_tfrz0
+         IF( rn_tfrz0<=0.e0 ) THEN
+            WRITE(ctmp1,*) '      bottom roughness must be strictly positive'
+            CALL ctl_stop( ctmp1 )
+         ENDIF
+         IF( ln_tfr2d ) THEN
+            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d
+            IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien
+         ENDIF
          !
          IF(ln_bfr2d) THEN
@@ -263,4 +332,15 @@
             bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable
          ENDIF
+
+         IF(ln_tfr2d) THEN
+            ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
+            CALL iom_open('tfr_coef.nc',inum)
+            CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! bfrcoef2d is used as tmp array
+            CALL iom_close(inum)
+            !
+            tfrcoef2d(:,:) = rn_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) )
+         ELSE
+            tfrcoef2d(:,:) = rn_tfri2  ! initialize tfrcoef2d to the namelist variable
+         ENDIF
          !
          IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all
@@ -279,4 +359,8 @@
                   bfrcoef2d(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp)
                   bfrcoef2d(ji,jj) = MIN(bfrcoef2d(ji,jj), rn_bfri2_max)
+                  ikbt = mikt(ji,jj)
+                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_tfrz0 ))**2._wp
+                  tfrcoef2d(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp)
+                  tfrcoef2d(ji,jj) = MIN(tfrcoef2d(ji,jj), rn_tfri2_max)
                END DO
             END DO
@@ -308,4 +392,6 @@
       zminbfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient
       zmaxbfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient
+      zmintfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient
+      zmaxtfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient
       !
 #  if defined key_vectopt_loop
@@ -339,4 +425,26 @@
              zminbfr = MIN(  zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) )  )
              zmaxbfr = MAX(  zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) )  )
+! (ISF)
+             ikbu = miku(ji,jj)       ! deepest ocean level at u- and v-points
+             ikbv = mikv(ji,jj)
+             zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt
+             zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt
+             IF( ABS( tfrcoef2d(ji,jj) ) > zfru ) THEN
+                IF( ln_ctl ) THEN
+                   WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu
+                   WRITE(numout,*) 'BFR ', ABS( tfrcoef2d(ji,jj) ), zfru
+                ENDIF
+                ictu = ictu + 1
+             ENDIF
+             IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN
+                 IF( ln_ctl ) THEN
+                     WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbv
+                     WRITE(numout,*) 'BFR ', tfrcoef2d(ji,jj), zfrv
+                 ENDIF
+                 ictv = ictv + 1
+             ENDIF
+             zmintfr = MIN(  zmintfr, MIN( zfru, ABS( tfrcoef2d(ji,jj) ) )  )
+             zmaxtfr = MAX(  zmaxtfr, MIN( zfrv, ABS( tfrcoef2d(ji,jj) ) )  )
+
          END DO
       END DO
@@ -346,4 +454,6 @@
          CALL mpp_min( zminbfr )
          CALL mpp_max( zmaxbfr )
+         CALL mpp_min( zmintfr )
+         CALL mpp_max( zmaxtfr )
       ENDIF
       IF( .NOT.ln_bfrimp) THEN
@@ -352,4 +462,5 @@
          WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '
          WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr
+         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr
          WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary'
       ENDIF
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 4666)
@@ -107,12 +107,13 @@
 
       !                                                ! ===============
-      DO jk = 2, jpkm1                                 ! Horizontal slab
+      !                                                ! Horizontal slab
          !                                             ! ===============
+      DO jj = 1, jpj                                     ! indicators:
+         DO ji = 1, jpi
+            DO jk = mikt(ji,jj)+1, jpkm1                                 ! Horizontal slab
          ! Define the mask 
          ! ---------------
-         rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau
-
-         DO jj = 1, jpj                                     ! indicators:
-            DO ji = 1, jpi
+               rrau(ji,jj,jk) = MAX( 1.e-20, rrau(ji,jj,jk) )         ! only retains positive value of rrau
+
                ! stability indicator: msks=1 if rn2>0; 0 elsewhere
                IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
@@ -136,8 +137,6 @@
                ELSE                                                         ;   zmskd3(ji,jj) = 1._wp
                ENDIF
-            END DO
-         END DO
          ! mask zmsk in order to have avt and avs masked
-         zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
+               zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)
 
 
@@ -145,8 +144,4 @@
          ! ------------------
          ! Constant eddy coefficient: reset to the background value
-!CDIR NOVERRCHK
-         DO jj = 1, jpj
-!CDIR NOVERRCHK
-            DO ji = 1, jpi
                zinr = 1./rrau(ji,jj,jk)
                ! salt fingering
@@ -165,14 +160,16 @@
             END DO
          END DO
-
-
+      END DO
          ! Increase avmu, avmv if necessary
          ! --------------------------------
 !!gm to be changed following the definition of avm.
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector opt.
+            DO jk = miku(ji,jj)+1, jpkm1                                 ! Horizontal slab
                avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
                   &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
                   &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk)
+            END DO
+            DO jk = mikv(ji,jj)+1, jpkm1
                avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
                   &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 4666)
@@ -241,5 +241,9 @@
       DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+            IF (mikt(ji,jj) .GT. 1) THEN
+               en(ji,jj,mikt(ji,jj))=rn_emin
+            ELSE
+               en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+            END IF
          END DO
       END DO
@@ -342,7 +346,7 @@
       END DO
       !
-      DO jk = 2, jpkm1           !* Matrix and right hand side in en
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jk = mikt(ji,jj)+1, jpkm1           !* Matrix and right hand side in en
                zcof   = zfact1 * tmask(ji,jj,jk)
                zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
@@ -360,43 +364,26 @@
                !                                   ! right hand side in en
                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
-      END DO
-      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
-      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
+                  &                                 * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
+            END DO
+            !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
+            DO jk = mikt(ji,jj)+2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
                zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
             END DO
-         END DO
-      END DO
-      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-         DO ji = fs_2, fs_jpim1    ! vector opt.
-            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
-         END DO
-      END DO
-      DO jk = 3, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+            ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+            zd_lw(ji,jj,mikt(ji,jj)+1) = en(ji,jj,mikt(ji,jj)+1) - zd_lw(ji,jj,mikt(ji,jj)+1) * en(ji,jj,mikt(ji,jj))    ! Surface boudary conditions on tke
+            !
+            DO jk = mikt(ji,jj)+2, jpkm1
                zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
             END DO
-         END DO
-      END DO
-      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-         DO ji = fs_2, fs_jpim1    ! vector opt.
+            !
+            ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
             en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
-         END DO
-      END DO
-      DO jk = jpk-2, 2, -1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+            !
+            DO jk = jpk-2, mikt(ji,jj)+1, -1
                en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
             END DO
-         END DO
-      END DO
-      DO jk = 2, jpkm1                             ! set the minimum value of tke
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            !
+            DO jk = mikt(ji,jj), jpkm1                             ! set the minimum value of tke
                en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
             END DO
@@ -412,5 +399,5 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
-                     &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+                     &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -421,5 +408,5 @@
                jk = nmln(ji,jj)
                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
-                  &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+                  &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1)
             END DO
          END DO
@@ -433,9 +420,9 @@
                   ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
                   zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
-                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )    ! module of the mean stress 
+                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress 
                   zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
-                     &                                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)
+                     &                                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1)
                END DO
             END DO
@@ -506,20 +493,33 @@
       !
       IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
-         zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
-         zmxlm(:,:,1) = MAX(  rn_mxl0,  zraug * taum(:,:)  )
-      ELSE                          ! surface set to the minimum value
-         zmxlm(:,:,1) = rn_mxl0
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               IF (mikt(ji,jj) .GT. 1) THEN
+                  zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min
+               ELSE
+                  zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
+                  zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) )
+               END IF
+            END DO
+         END DO
+      ELSE 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1                         ! surface set to the minimum value
+               zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min)
+            END DO
+         END DO
       ENDIF
       zmxlm(:,:,jpk)  = rmxl_min     ! last level set to the interior minium value
       !
 !CDIR NOVERRCHK
-      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
-!CDIR NOVERRCHK
-         DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            !CDIR NOVERRCHK
+            DO jk = mikt(ji,jj)+1, jpkm1              ! interior value : l=sqrt(2*e/n^2)
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
                zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
             END DO
+            zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj))   ! surface set to the minimum value 
          END DO
       END DO
@@ -533,8 +533,8 @@
       !
       CASE ( 0 )           ! bounded by the distance to surface and bottom
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               DO jk = mikt(ji,jj)+1, jpkm1
+                  zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
                   &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) )
                   zmxlm(ji,jj,jk) = zemxl
@@ -545,7 +545,7 @@
          !
       CASE ( 1 )           ! bounded by the vertical scale factor
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               DO jk = mikt(ji,jj)+1, jpkm1
                   zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
                   zmxlm(ji,jj,jk) = zemxl
@@ -556,14 +556,10 @@
          !
       CASE ( 2 )           ! |dk[xml]| bounded by e3t :
-         DO jk = 2, jpkm1         ! from the surface to the bottom :
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom :
                   zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
                END DO
-            END DO
-         END DO
-         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+               DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface :
                   zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
                   zmxlm(ji,jj,jk) = zemxl
@@ -574,14 +570,10 @@
          !
       CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
-         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : lup
                   zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
                END DO
-            END DO
-         END DO
-         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+               DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : ldown
                   zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
                END DO
@@ -628,18 +620,20 @@
       CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
       !
-      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jk = miku(ji,jj)+1, jpkm1            !* vertical eddy viscosity at u- and v-points
+               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
+            END DO
+            DO jk = mikv(ji,jj)+1, jpkm1
+               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
+      !
+      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               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
-      !
-      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+               DO jk = mikt(ji,jj)+1, jpkm1
                   zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
                   !                                          ! shear
@@ -655,6 +649,6 @@
                   avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
 # if defined key_c1d
-                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
-                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
+                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number
+                  e_ric(ji,jj,jk) = zri   * tmask(ji,jj,jk)  ! c1d config. : save Ri
 # endif
               END DO
Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 4666)
@@ -126,5 +126,5 @@
       zkz(:,:) = 0.e0               !* Associated potential energy consummed over the whole water column
       DO jk = 2, jpkm1
-         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk)
+         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk) * tmask(:,:,jk-1)
       END DO
 
@@ -135,6 +135,10 @@
       END DO
 
-      DO jk = 2, jpkm1              !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
-         zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s
+      DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
+         DO ji = 1, jpi
+            DO jk = mikt(ji,jj)+1, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
+               zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s
+            END DO
+         END DO
       END DO
 
@@ -162,9 +166,16 @@
       !                          !   Update  mixing coefs  !                          
       !                          ! ----------------------- !
-      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
-         avt(:,:,jk) = avt(:,:,jk) + zav_tide(:,:,jk)
-         avm(:,:,jk) = avm(:,:,jk) + zav_tide(:,:,jk)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1  ! vector opt.
+      DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
+         DO ji = 1, jpi
+            DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
+               avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk)
+               avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1  ! vector opt.
+            DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
                avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
                avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
@@ -237,6 +248,6 @@
       zsum2(:,:) = 0.e0
       DO jk= 2, jpk
-         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk)
-         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk)                
+         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)
+         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)               
       END DO
       DO jj = 1, jpj
@@ -274,5 +285,5 @@
       zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column
       DO jk = 2, jpkm1
-         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk)
+         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)
       END DO
 
@@ -284,5 +295,5 @@
 
       DO jk = 2, jpkm1              ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s
-         zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. )   ! kz max = 120 cm2/s
+         zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) * tmask(:,:,jk) * tmask(:,:,jk-1)   ! kz max = 120 cm2/s
       END DO
 
@@ -377,5 +388,5 @@
       READ  ( numnam_cfg, namzdf_tmx, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namzdf_tmx )
+      WRITE ( numond, namzdf_tmx )
 
       IF(lwp) THEN                   ! Control print
@@ -414,5 +425,5 @@
       ! only the energy available for mixing is taken into account,
       ! (mixing efficiency tidal dissipation efficiency)
-      en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * tmask(:,:,1)
+      en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * lmask(:,:)
 
       ! Vertical structure (az_tmx)
@@ -443,7 +454,7 @@
          ztpc = 0.e0
          zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:)
-         DO jk= 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk= mikt(ji,jj)+1, jpkm1
                   ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
@@ -459,10 +470,10 @@
          zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )   
          zkz(:,:) = 0.e0
-         DO jk = 2, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk = mikt(ji,jj)+1, jpkm1
                zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* tmask(ji,jj,jk)
-            END DO
-         END DO
+               END DO
+            END DO
          END DO
          ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz
@@ -484,6 +495,10 @@
          WRITE(numout,*) '          Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) )
 
-         DO jk = 2, jpkm1
-            zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk = mikt(ji,jj)+1, jpkm1
+                  zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s
+               END DO
+            END DO
          END DO
          ztpc = 0.e0
