Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5329)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 5332)
@@ -472,25 +472,4 @@
          risfdep(:,:)=0.e0
          misfdep(:,:)=1
-         !
-         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code
-         IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN 
-            risfdep(:,:)=200.e0 
-            misfdep(:,:)=1 
-            ij0 = 1 ; ij1 = 40 
-            DO jj = mj0(ij0), mj1(ij1) 
-               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
-            END DO 
-            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
-         ! 
-         ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN
-         ! 
-            risfdep(:,:)=0.e0
-            misfdep(:,:)=1
-            ij0 = 1 ; ij1 = 40
-            DO jj = mj0(ij0), mj1(ij1)
-               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
-            END DO
-            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp
-         END IF
          !
          DEALLOCATE( idta, zdta )
@@ -969,5 +948,5 @@
       !!
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
-      INTEGER  ::   ik, it           ! temporary integers
+      INTEGER  ::   ik, it, ikb, ikt ! temporary integers
       LOGICAL  ::   ll_print         ! Allow  control print for debugging
       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
@@ -1152,17 +1131,18 @@
       IF ( ln_isfcav ) THEN
       ! (ISF) define e3uw (adapted for 2 cells in the water column)
-      ! Need to test if the modification of only mikt and mbkt levels is enough
-         DO jk = 2,jpk                          
-            DO jj = 1, jpjm1 
-               DO ji = 1, fs_jpim1   ! vector opt. 
-                  e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) &
-                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) )
-                  e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) &
-                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) )
-               END DO 
-            END DO 
+         DO jj = 2, jpjm1 
+            DO ji = 2, fs_jpim1   ! vector opt. 
+               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
+               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
+               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
+               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
+            END DO
          END DO
       END IF
-      
+
       CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
       CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
@@ -1538,5 +1518,5 @@
  
  ! remove single point "bay" on isf coast line in the ice shelf draft'
-         DO jk = 1, jpk
+         DO jk = 2, jpk
             WHERE (misfdep==0) misfdep=jpk
             zmask=0
Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 5329)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 5332)
@@ -110,9 +110,4 @@
          ELSEIF( cp_cfg == 'gyre' ) THEN         
             CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields
-        ELSEIF( cp_cfg == 'isomip' .OR. cp_cfg == 'isomip2') THEN
-            IF(lwp) WRITE(numout,*) 'Initialization of T+S for ISOMIP domain' 
-            tsn(:,:,:,jp_tem)=-1.9*tmask(:,:,:)          ! ISOMIP configuration : start from constant T+S fields 
-            tsn(:,:,:,jp_sal)=34.4*tmask(:,:,:)
-            tsb(:,:,:,:)=tsn(:,:,:,:)  
          ELSE                                    ! Initial T-S, U-V fields read in files
             IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000
Index: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 5329)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 5332)
@@ -171,10 +171,12 @@
             END DO
          END DO
+         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
+
          IF ( ln_isfcav ) THEN
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
                   ! (ISF) ========================================================================
-                  ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points 
-                  ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points)
+                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
+                  ikbv = mikv(ji,jj)         ! (1st wet ocean u- and v-points)
                   !
                   zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     &
@@ -183,6 +185,6 @@
                      &           + 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 )
+                  zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_tfeb2 )
+                  zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_tfeb2 )
               !
                   tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1))
@@ -202,7 +204,7 @@
                END DO
             END DO
-         END IF
-         !
-         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
+            CALL lbc_lnk( tfrua, 'U', 1. )   ;   CALL lbc_lnk( tfrva, 'V', 1. )      ! Lateral boundary condition
+         END IF
+         !
          !
          IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        &
@@ -277,9 +279,11 @@
             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
+         IF ( ln_isfcav ) THEN
+            IF(lwp) WRITE(numout,*) '      top    friction coef.   rn_bfri1  = ', rn_tfri1
+            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
+         END IF
          !
          IF(ln_bfr2d) THEN
@@ -295,4 +299,19 @@
          bfrua(:,:) = - bfrcoef2d(:,:)
          bfrva(:,:) = - bfrcoef2d(:,:)
+         !
+         IF ( ln_isfcav ) THEN
+            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(:,:)
+         END IF
          !
       CASE( 2 )
@@ -311,18 +330,20 @@
             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_isfcav ) THEN
+            IF(lwp) WRITE(numout,*) '      quadratic top    friction'
+            IF(lwp) WRITE(numout,*) '      friction coef.    rn_tfri2     = ', 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,*) '      top roughness     rn_tfrz0 [m] = ', rn_tfrz0
+            IF( rn_tfrz0<=0.e0 ) THEN
+               WRITE(ctmp1,*) '      top 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
+         END IF
          !
          IF(ln_bfr2d) THEN
@@ -336,4 +357,17 @@
             bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable
          ENDIF
+         
+         IF ( ln_isfcav ) THEN
+            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_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) )
+            ELSE
+               tfrcoef2d(:,:) = rn_tfri2  ! initialize tfrcoef2d to the namelist variable
+            ENDIF
+         END IF
          !
          IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all
@@ -346,4 +380,14 @@
                END DO
             END DO
+            IF ( ln_isfcav ) THEN
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     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
+            END IF
          ENDIF
          !
@@ -398,4 +442,28 @@
              zminbfr = MIN(  zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) )  )
              zmaxbfr = MAX(  zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) )  )
+! (ISF)
+             IF ( ln_isfcav ) THEN
+                ikbu = miku(ji,jj)       ! 1st wet 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,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbu
+                      WRITE(numout,*) 'TFR ', ABS( tfrcoef2d(ji,jj) ), zfru
+                   ENDIF
+                   ictu = ictu + 1
+                ENDIF
+                IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN
+                   IF( ln_ctl ) THEN
+                      WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbv
+                      WRITE(numout,*) 'TFR ', 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 IF
+! END ISF
          END DO
       END DO
@@ -405,12 +473,14 @@
          CALL mpp_min( zminbfr )
          CALL mpp_max( zmaxbfr )
+         IF ( ln_isfcav) CALL mpp_min( zmintfr )
+         IF ( ln_isfcav) CALL mpp_max( zmaxtfr )
       ENDIF
       IF( .NOT.ln_bfrimp) THEN
       IF( lwp .AND. ictu + ictv > 0 ) THEN
-         WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points '
-         WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '
+         WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictu, ' U-points '
+         WRITE(numout,*) ' Bottom/Top 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'
+         IF ( ln_isfcav ) WRITE(numout,*) ' Top friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr
+         WRITE(numout,*) ' Bottom/Top friction coefficient will be reduced where necessary'
       ENDIF
       ENDIF
