Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 4267)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 4278)
@@ -8,4 +8,5 @@
    !!            3.3  !  2010-09  (D. Storkey) add ice boundary conditions
    !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
+   !!            3.6  !  2012-01  (C. Rousset) add ice boundary conditions for lim3
    !!----------------------------------------------------------------------
 #if defined key_bdy 
@@ -62,4 +63,11 @@
       REAL(wp), POINTER, DIMENSION(:)     ::  hicif
       REAL(wp), POINTER, DIMENSION(:)     ::  hsnif
+#elif defined key_lim3
+      LOGICAL                         ::  ll_a_i
+      LOGICAL                         ::  ll_ht_i
+      LOGICAL                         ::  ll_ht_s
+      REAL, POINTER, DIMENSION(:,:)   ::  a_i   !: now ice leads fraction climatology
+      REAL, POINTER, DIMENSION(:,:)   ::  ht_i  !: Now ice  thickness climatology
+      REAL, POINTER, DIMENSION(:,:)   ::  ht_s  !: now snow thickness
 #endif
    END TYPE OBC_DATA
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 4267)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 4278)
@@ -11,4 +11,5 @@
    !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
    !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
+   !!            3.6  !  2012-01  (C. Rousset) add ice boundary conditions for lim3
    !!----------------------------------------------------------------------
 #if defined key_bdy
@@ -32,4 +33,8 @@
 #if defined key_lim2
    USE ice_2
+#elif defined key_lim3
+   USE par_ice
+   USE ice
+   USE limcat_1D          ! redistribute ice input into categories
 #endif
    USE sbcapr
@@ -50,4 +55,9 @@
 
    TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr   ! array of pointers to nbmap
+
+#if defined key_lim3
+   LOGICAL :: ll_bdylim3                  ! determine whether ice input is lim2 (F) or lim3 (T) type
+   INTEGER :: jfld_hti, jfld_hts, jfld_ai ! indices of ice thickness, snow thickness and concentration in bf structure
+#endif
 
 #  include "domzgr_substitute.h90"
@@ -78,5 +88,5 @@
                                                         ! etc.
       !!
-      INTEGER     ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd  ! local indices
+      INTEGER     ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd, jl  ! local indices
       INTEGER,          DIMENSION(jpbgrd) ::   ilen1 
       INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts
@@ -213,4 +223,38 @@
                   END DO 
                END IF
+            ENDIF
+#elif defined key_lim3
+            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 
+               ilen1(:) = nblen(:)
+               IF( dta%ll_a_i ) THEN
+                  igrd = 1   
+                  DO jl = 1, jpl
+                     DO ib = 1, ilen1(igrd)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        dta_bdy(ib_bdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1) 
+                     END DO
+                  END DO
+               ENDIF
+               IF( dta%ll_ht_i ) THEN
+                  igrd = 1   
+                  DO jl = 1, jpl
+                     DO ib = 1, ilen1(igrd)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        dta_bdy(ib_bdy)%ht_i (ib,jl) =  ht_i(ii,ij,jl) * tmask(ii,ij,1) 
+                     END DO
+                  END DO
+               ENDIF
+               IF( dta%ll_ht_s ) THEN
+                  igrd = 1   
+                  DO jl = 1, jpl
+                     DO ib = 1, ilen1(igrd)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        dta_bdy(ib_bdy)%ht_s (ib,jl) =  ht_s(ii,ij,jl) * tmask(ii,ij,1) 
+                     END DO
+                  END DO
+               ENDIF
             ENDIF
 #endif
@@ -346,4 +390,10 @@
 
                ENDIF
+#if defined key_lim3
+               IF( .NOT. ll_bdylim3 .AND. nn_ice_lim(ib_bdy) > 0 .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is lim2 type)
+     	          CALL lim_cat_1D ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
+                                  & dta_bdy(ib_bdy)%ht_i,     dta_bdy(ib_bdy)%ht_s,     dta_bdy(ib_bdy)%a_i     )
+               ENDIF
+#endif
             ENDIF
             jstart = jstart + dta%nread(1)
@@ -400,4 +450,9 @@
       INTEGER, POINTER, DIMENSION(:)         ::   nblen, nblenrim  ! short cuts
       TYPE(OBC_DATA), POINTER                ::   dta           ! short cut
+#if defined key_lim3
+      INTEGER, DIMENSION(3) ::   zdimsz   ! number of elements in each of the 4 dimensions (i.e. i,j,t,ice-cat) for an array
+      INTEGER               ::   zndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat)
+      INTEGER               ::   inum,id1 ! local integer
+#endif
       TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   blf_i         !  array of namelist information structures
       TYPE(FLD_N) ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! 
@@ -405,8 +460,12 @@
 #if defined key_lim2
       TYPE(FLD_N) ::   bn_frld, bn_hicif, bn_hsnif      !
+#elif defined key_lim3
+      TYPE(FLD_N) ::   bn_a_i, bn_ht_i, bn_ht_s      
 #endif
       NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 
 #if defined key_lim2
       NAMELIST/nambdy_dta/ bn_frld, bn_hicif, bn_hsnif
+#elif defined key_lim3
+      NAMELIST/nambdy_dta/ bn_a_i, bn_ht_i, bn_ht_s
 #endif
       NAMELIST/nambdy_dta/ ln_full_vel
@@ -425,6 +484,6 @@
                                ,nn_dyn3d_dta(ib_bdy)       &
                                ,nn_tra_dta(ib_bdy)         &
-#if defined key_lim2
-                               ,nn_ice_lim2_dta(ib_bdy)    &
+#if ( defined key_lim2 || defined key_lim3 )
+                              ,nn_ice_lim_dta(ib_bdy)    &
 #endif
                               )
@@ -446,6 +505,6 @@
             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2
          ENDIF
-#if defined key_lim2
-         IF( cn_ice_lim2(ib_bdy) /= 'none' .and. nn_ice_lim2_dta(ib_bdy) .eq. 1  ) THEN
+#if ( defined key_lim2 || defined key_lim3 )
+         IF( cn_ice_lim(ib_bdy) /= 'none' .and. nn_ice_lim_dta(ib_bdy) .eq. 1  ) THEN
             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3
          ENDIF
@@ -495,4 +554,8 @@
             bn_hicif   = FLD_N(  'bdy_hicif'   ,  24   , 'iicethic' , .false. , .false. ,   'yearly'  , ''       , ''        )
             bn_hsnif   = FLD_N(  'bdy_hsnif'   ,  24   , 'isnothic' , .false. , .false. ,   'yearly'  , ''       , ''        )
+#elif defined key_lim3
+            bn_a_i  = FLD_N(  'bdy_a_i'     ,  24   , 'ileadfra' , .false. , .false. ,   'yearly'  , ''       , ''        )
+            bn_ht_i = FLD_N(  'bdy_ht_i'    ,  24   , 'iicethic' , .false. , .false. ,   'yearly'  , ''       , ''        )
+            bn_ht_s = FLD_N(  'bdy_ht_s'    ,  24   , 'isnowthi' , .false. , .false. ,   'yearly'  , ''       , ''        )
 #endif
 
@@ -633,4 +696,49 @@
 
             ENDIF
+#elif defined key_lim3
+            ! sea ice
+            IF( nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN
+
+               ! Test for types of ice input (lim2 or lim3) 
+               CALL iom_open ( bn_a_i%clname, inum )
+               id1 = iom_varid ( inum, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. )
+               CALL iom_close ( inum )
+               !CALL fld_clopn ( bn_a_i, nyear, nmonth, nday, ldstop=.TRUE. )
+               !CALL iom_open ( bn_a_i %clname, inum )
+               !id1 = iom_varid ( bn_a_i%num, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. )
+      	       IF ( zndims == 4 ) THEN
+                 ll_bdylim3 = .TRUE.   ! lim3 input
+               ELSE
+                 ll_bdylim3 = .FALSE.  ! lim2 input      
+               ENDIF
+               ! End test
+
+               IF( dta%ll_a_i ) THEN
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_a_i
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 1
+                  ilen1(jfld) = nblen(igrid(jfld))
+                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF
+               ENDIF
+
+               IF( dta%ll_ht_i ) THEN
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_ht_i
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 1
+                  ilen1(jfld) = nblen(igrid(jfld))
+                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF
+               ENDIF
+
+               IF( dta%ll_ht_s ) THEN
+                  jfld = jfld + 1
+                   blf_i(jfld) = bn_ht_s
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 1
+                  ilen1(jfld) = nblen(igrid(jfld))
+                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF
+               ENDIF
+
 #endif
             ! Recalculate field counts
@@ -758,5 +866,5 @@
 
 #if defined key_lim2
-         IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN
+         IF (nn_ice_lim(ib_bdy) .gt. 0) THEN
             IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN
                ALLOCATE( dta_bdy(ib_bdy)%frld(nblen(1)) )
@@ -770,4 +878,33 @@
                jfld = jfld + 1
                dta_bdy(ib_bdy)%hsnif => bf(jfld)%fnow(:,1,1)
+            ENDIF
+         ENDIF
+#elif defined key_lim3
+         IF (nn_ice_lim(ib_bdy) .gt. 0) THEN
+            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN
+               ALLOCATE( dta_bdy(ib_bdy)%a_i (nblen(1),jpl) )
+               ALLOCATE( dta_bdy(ib_bdy)%ht_i(nblen(1),jpl) )
+               ALLOCATE( dta_bdy(ib_bdy)%ht_s(nblen(1),jpl) )
+            ELSE
+               IF ( ll_bdylim3 ) THEN ! case input is lim3 type
+                  jfld = jfld + 1
+                  dta_bdy(ib_bdy)%a_i  => bf(jfld)%fnow(:,1,:)
+                  jfld = jfld + 1
+                  dta_bdy(ib_bdy)%ht_i => bf(jfld)%fnow(:,1,:)
+                  jfld = jfld + 1
+                  dta_bdy(ib_bdy)%ht_s => bf(jfld)%fnow(:,1,:)
+               ELSE ! case input is lim2 type
+                  jfld_ai  = jfld + 1
+                  jfld_hti = jfld + 2
+                  jfld_hts = jfld + 3
+                  jfld     = jfld + 3
+                  ALLOCATE( dta_bdy(ib_bdy)%a_i (nblen(1),jpl) )
+                  ALLOCATE( dta_bdy(ib_bdy)%ht_i(nblen(1),jpl) )
+                  ALLOCATE( dta_bdy(ib_bdy)%ht_s(nblen(1),jpl) )
+                  dta_bdy(ib_bdy)%a_i (:,:) = 0.0
+                  dta_bdy(ib_bdy)%ht_i(:,:) = 0.0
+                  dta_bdy(ib_bdy)%ht_s(:,:) = 0.0
+               ENDIF
+
             ENDIF
          ENDIF
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 4278)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 4278)
@@ -0,0 +1,410 @@
+MODULE bdyice_lim
+   !!======================================================================
+   !!                       ***  MODULE  bdyice_lim  ***
+   !! Unstructured Open Boundary Cond. :  Open boundary conditions for sea-ice (LIM2 and LIM3)
+   !!======================================================================
+   !!  History :  3.3  !  2010-09 (D. Storkey)  Original code
+   !!             3.4  !  2011    (D. Storkey)  rewrite in preparation for OBC-BDY merge
+   !!              -   !  2012-01 (C. Rousset)  add lim3 and remove useless jk loop 
+   !!----------------------------------------------------------------------
+#if defined   key_bdy   &&  ( defined key_lim2 || defined key_lim3 )
+   !!----------------------------------------------------------------------
+   !!   'key_bdy'            and                 Unstructured Open Boundary Conditions
+   !!   'key_lim2'                                                 LIM-2 sea ice model
+   !!   'key_lim3'                                                 LIM-3 sea ice model
+   !!----------------------------------------------------------------------
+   !!   bdy_ice_lim        : Application of open boundaries to ice
+   !!   bdy_ice_frs        : Application of Flow Relaxation Scheme
+   !!----------------------------------------------------------------------
+   USE timing          ! Timing
+   USE phycst          ! physical constant
+   USE eosbn2          ! equation of state
+   USE oce             ! ocean dynamics and tracers variables
+#if defined key_lim2
+   USE par_ice_2
+   USE ice_2           ! LIM_2 ice variables
+#elif defined key_lim3
+   USE par_ice
+   USE ice             ! LIM_3 ice variables
+#endif 
+   USE par_oce         ! ocean parameters
+   USE dom_oce         ! ocean space and time domain variables 
+   USE dom_ice          ! sea-ice domain
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE bdy_oce         ! ocean open boundary conditions
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager  ! write to numout file
+   USE lib_mpp         ! distributed memory computing
+   USE lib_fortran     ! to use key_nosignedzero
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   bdy_ice_lim    ! routine called in sbcmod
+   PUBLIC   bdy_ice_lim_dyn ! routine called in limrhg
+
+   REAL(wp) ::   epsi20 = 1.e-20_wp  ! module constants
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: bdyice.F90 2715 2011-03-30 15:58:35Z rblod $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE bdy_ice_lim( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_ice_lim  ***
+      !!
+      !! ** Purpose : - Apply open boundary conditions for ice (LIM2 and LIM3)
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt     ! Main time step counter
+      !!
+      INTEGER               :: ib_bdy ! Loop index
+      DO ib_bdy=1, nb_bdy
+
+         SELECT CASE( cn_ice_lim(ib_bdy) )
+         CASE('none')
+            CYCLE
+         CASE('frs')
+            CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+         CASE DEFAULT
+            CALL ctl_stop( 'bdy_ice_lim : unrecognised option for open boundaries for ice fields' )
+         END SELECT
+      ENDDO
+
+   END SUBROUTINE bdy_ice_lim
+
+   SUBROUTINE bdy_ice_frs( idx, dta, kt, ib_bdy )
+      !!------------------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_ice_frs  ***
+      !!                    
+      !! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields in the case 
+      !!              of unstructured open boundaries.
+      !! 
+      !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three-
+      !!             dimensional baroclinic ocean model with realistic topography. Tellus, 365-382.
+      !!------------------------------------------------------------------------------
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      INTEGER,         INTENT(in) ::   kt   ! main time-step counter
+      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index      !!
+
+      INTEGER  ::   jb, jk, jgrd, jl   ! dummy loop indices
+      INTEGER  ::   ji, jj             ! local scalar
+      REAL(wp) ::   zwgt, zwgt1        ! local scalar
+      REAL(wp) ::   zinda, ztmelts
+      !!------------------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_ice_frs')
+      !
+      !IF( kt /= nit000 ) THEN
+      jgrd = 1      ! Everything is at T-points here
+      !
+#if defined key_lim2
+      DO jb = 1, idx%nblen(jgrd)
+         ji    = idx%nbi(jb,jgrd)
+         jj    = idx%nbj(jb,jgrd)
+         zwgt  = idx%nbw(jb,jgrd)
+         zwgt1 = 1.e0 - idx%nbw(jb,jgrd)
+         frld (ji,jj) = ( frld (ji,jj) * zwgt1 + dta%frld (jb) * zwgt ) * tmask(ji,jj,1)     ! Leads fraction 
+         hicif(ji,jj) = ( hicif(ji,jj) * zwgt1 + dta%hicif(jb) * zwgt ) * tmask(ji,jj,1)     ! Ice depth 
+         hsnif(ji,jj) = ( hsnif(ji,jj) * zwgt1 + dta%hsnif(jb) * zwgt ) * tmask(ji,jj,1)     ! Snow depth
+
+!         zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - frld(ji,jj) ) ) ! 0 if no ice
+!         !------------------------------
+!         ! Sea ice surface temperature
+!         !------------------------------
+!         sist(ji,jj)   = zinda * 270.0 + ( 1.0 - zinda ) * tfu(ji,jj)
+!         !-----------------------------------------------
+!         ! Ice/snow temperatures and energy stored in brines 
+!         !-----------------------------------------------
+!         !!! TO BE CONTIUNED (as LIM3 below) !!!
+!            zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH
+! 
+!               ! Recover in situ values.
+!               zindb         = MAX( rzero, SIGN( rone, zs0a(ji,jj) - epsi06 ) )
+!               zacrith       = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) )
+!               zs0a (ji,jj)  = zindb * MIN( zs0a(ji,jj), zacrith )
+!               hsnif(ji,jj)  = zindb * ( zs0sn(ji,jj) /MAX( zs0a(ji,jj), epsi16 ) )
+!               hicif(ji,jj)  = zindb * ( zs0ice(ji,jj)/MAX( zs0a(ji,jj), epsi16 ) )
+!               zindsn        = MAX( rzero, SIGN( rone, hsnif(ji,jj) - epsi06 ) )
+!               zindic        = MAX( rzero, SIGN( rone, hicif(ji,jj) - epsi03 ) )
+!               zindb         = MAX( zindsn, zindic )
+!               zs0a (ji,jj)  = zindb * zs0a(ji,jj)
+!               frld (ji,jj)  = 1.0 - zs0a(ji,jj)
+!               hsnif(ji,jj)  = zindsn * hsnif(ji,jj)
+!               hicif(ji,jj)  = zindic * hicif(ji,jj)
+!               zusvosn       = 1.0/MAX( hsnif(ji,jj) * zs0a(ji,jj), epsi16 )
+!               zusvoic       = 1.0/MAX( hicif(ji,jj) * zs0a(ji,jj), epsi16 )
+!               zignm         = MAX( rzero,  SIGN( rone, hsndif - hsnif(ji,jj) ) )
+!               zrtt          = 173.15 * rone 
+!               ztsn          =          zignm   * tbif(ji,jj,1)  &
+!                              + ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) ) 
+!               ztic1          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c1(ji,jj) ) , tfu(ji,jj) )
+!               ztic2          = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c2(ji,jj) ) , tfu(ji,jj) )
+! 
+!               tbif(ji,jj,1) = zindsn * ztsn  + ( 1.0 - zindsn ) * tfu(ji,jj)               
+!               tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj)
+!               tbif(ji,jj,3) = zindic * ztic2 + ( 1.0 - zindic ) * tfu(ji,jj)
+!               qstoif(ji,jj) = zindb  * xlic * zs0st(ji,jj) /  MAX( zs0a(ji,jj), epsi16 )
+
+      END DO 
+
+      CALL lbc_bdy_lnk( frld,  'T', 1., ib_bdy )                                         ! lateral boundary conditions
+      CALL lbc_bdy_lnk( hicif, 'T', 1., ib_bdy )
+      CALL lbc_bdy_lnk( hsnif, 'T', 1., ib_bdy )
+
+      vt_i(:,:) = hicif(:,:) * frld(:,:)
+      vt_s(:,:) = hsnif(:,:) * frld(:,:)
+      !
+#elif defined key_lim3
+
+      DO jl = 1, jpl
+         DO jb = 1, idx%nblen(jgrd)
+            ji    = idx%nbi(jb,jgrd)
+            jj    = idx%nbj(jb,jgrd)
+            zwgt  = idx%nbw(jb,jgrd)
+            zwgt1 = 1.e0 - idx%nbw(jb,jgrd)
+            a_i (ji,jj,jl) = ( a_i (ji,jj,jl) * zwgt1 + dta%a_i (jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Leads fraction 
+            ht_i(ji,jj,jl) = ( ht_i(ji,jj,jl) * zwgt1 + dta%ht_i(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice depth 
+            ht_s(ji,jj,jl) = ( ht_s(ji,jj,jl) * zwgt1 + dta%ht_s(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth
+         ENDDO
+         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy )
+
+         DO jb = 1, idx%nblen(jgrd)
+            ji    = idx%nbi(jb,jgrd)
+            jj    = idx%nbj(jb,jgrd)
+            ! clem : condition on ice thickness depends on the ice velocity
+            ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values 
+            IF ( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) THEN
+               ht_i(ji,jj,jl) = ht_i(ji+1,jj  ,jl)
+               a_i (ji,jj,jl) = a_i (ji+1,jj  ,jl)
+               ht_s(ji,jj,jl) = ht_s(ji+1,jj  ,jl)
+            ENDIF
+            IF ( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) THEN
+               ht_i(ji,jj,jl) = ht_i(ji-1,jj  ,jl)
+               a_i (ji,jj,jl) = a_i (ji-1,jj  ,jl)
+               ht_s(ji,jj,jl) = ht_s(ji-1,jj  ,jl)
+            ENDIF
+            IF ( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) THEN
+               ht_i(ji,jj,jl) = ht_i(ji  ,jj+1,jl)
+               a_i (ji,jj,jl) = a_i (ji  ,jj+1,jl)
+               ht_s(ji,jj,jl) = ht_s(ji  ,jj+1,jl)
+            ENDIF
+            IF ( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) THEN
+               ht_i(ji,jj,jl) = ht_i(ji  ,jj-1,jl)
+               a_i (ji,jj,jl) = a_i (ji  ,jj-1,jl)
+               ht_s(ji,jj,jl) = ht_s(ji  ,jj-1,jl)
+            ENDIF
+
+            zinda = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - a_i(ji,jj,jl) ) ) ! 0 if no ice
+            ! --------------------
+            ! Ice and snow volumes
+            ! --------------------
+            v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
+            v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl)
+
+            ! -------------
+            ! Ice salinity
+            !---------------
+            sm_i(ji,jj,jl)   = zinda * bulk_sal  + ( 1.0 - zinda ) * s_i_min
+            smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
+            
+            !----------
+            ! Ice age
+            !----------
+            o_i(ji,jj,jl)    = zinda * 1.0   + ( 1.0 - zinda )
+            oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl)
+            
+            !------------------------------
+            ! Sea ice surface temperature
+            !------------------------------
+            t_su(ji,jj,jl)   = zinda * 270.0 + ( 1.0 - zinda ) * t_bo(ji,jj)
+            
+            !------------------------------------
+            ! Snow temperature and heat content
+            !------------------------------------
+            DO jk = 1, nlay_s
+               t_s(ji,jj,jk,jl) = zinda * 270.00 + ( 1.0 - zinda ) * rtt
+               ! Snow energy of melting
+               e_s(ji,jj,jk,jl) = zinda * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
+               ! Change dimensions
+               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
+               ! Multiply by volume, so that heat content in 10^9 Joules
+               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
+                    v_s(ji,jj,jl)  / nlay_s
+            END DO !jk
+            
+            !-----------------------------------------------
+            ! Ice salinities, temperature and heat content 
+            !-----------------------------------------------
+            DO jk = 1, nlay_i
+               t_i(ji,jj,jk,jl) = zinda * 270.00 + ( 1.0 - zinda ) * rtt 
+               s_i(ji,jj,jk,jl) = zinda * bulk_sal + ( 1.0 - zinda ) * s_i_min
+               ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
+               
+               ! heat content per unit volume
+               e_i(ji,jj,jk,jl) = zinda * rhoic * &
+                    (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
+                    +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) &
+                    - rcp      * ( ztmelts - rtt ) )
+               
+               ! Correct dimensions to avoid big values
+               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
+               
+               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
+               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 
+                    area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i
+            END DO ! jk
+
+         END DO !jb
+ 
+
+         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )                                         ! lateral boundary conditions
+         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk(  v_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk(  v_s(:,:,jl), 'T', 1., ib_bdy )
+ 
+         CALL lbc_bdy_lnk( smv_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk(  sm_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk(  oa_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk(   o_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk(  t_su(:,:,jl), 'T', 1., ib_bdy )
+         DO jk = 1, nlay_s
+            CALL lbc_bdy_lnk(t_s(:,:,jk,jl), 'T', 1., ib_bdy )
+            CALL lbc_bdy_lnk(e_s(:,:,jk,jl), 'T', 1., ib_bdy )
+         END DO
+         DO jk = 1, nlay_i
+            CALL lbc_bdy_lnk(t_i(:,:,jk,jl), 'T', 1., ib_bdy )
+            CALL lbc_bdy_lnk(e_i(:,:,jk,jl), 'T', 1., ib_bdy )
+         END DO
+
+      END DO !jl
+    
+#endif
+      !      
+      !ENDIF !nit000/=0
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_frs')
+      !
+   END SUBROUTINE bdy_ice_frs
+
+
+   SUBROUTINE bdy_ice_lim_dyn( kn, pvar1, pvar2, pvar12 )
+      !!------------------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_ice_lim_dyn  ***
+      !!                    
+      !! ** Purpose : Apply dynamics boundary conditions for sea-ice in the cas of unstructured open boundaries.
+      !!      kn = 1: u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free
+      !!              if adjacent grid point is ice free, then u_ice and v_ice are equal to ocean velocities
+      !!      kn = 2: the stress tensor is set to 0 (i.e. pvar1, pvar2, pvar12)
+      !!
+      !! 2013-06 : C. Rousset
+      !!------------------------------------------------------------------------------
+      !!
+      INTEGER,  INTENT( in    )                           ::   kn     ! set up of ice vel (kn=1) or stress tensor (kn=2)
+      REAL(wp), INTENT( inout ), DIMENSION(:,:), OPTIONAL ::   pvar1, pvar2, pvar12  ! stress tensor components (zs1,zs2,zs12) 
+      INTEGER  ::   jb, jgrd   ! dummy loop indices
+      INTEGER  ::   ji, jj             ! local scalar
+      INTEGER  ::   ib_bdy ! Loop index
+      REAL(wp) ::   zmsk1, zmsk2, zflag 
+     !!------------------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_ice_lim_dyn')
+      !
+      DO ib_bdy=1, nb_bdy
+      !
+      SELECT CASE( nn_ice_lim(ib_bdy) )
+
+      CASE(jp_none)
+         CYCLE
+
+      CASE(jp_frs)
+
+         IF ( kn == 1 ) THEN ! set up of u_ice and v_ice
+           
+            jgrd = 2      ! u velocity
+            DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd)
+               ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
+               jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
+               zflag = idx_bdy(ib_bdy)%flagu(jb)
+
+               IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries
+                  ! one of the two zmsk is always 0 (because of zflag)
+                  zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice
+                  zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice
+                  
+                  ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
+                  u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + &
+                     &            u_ice(ji-1,jj) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + &
+                     &            u_oce(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
+               ELSE                             ! everywhere else
+                  u_ice(ji,jj) = u_oce(ji,jj)
+               ENDIF
+
+            ENDDO
+            
+            jgrd = 3      ! v velocity
+            DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd)
+               ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
+               jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
+               zflag = idx_bdy(ib_bdy)%flagv(jb)
+              
+               IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries
+                  ! one of the two zmsk is always 0 (because of zflag)
+                  zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice
+                  zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice
+                  
+                  ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
+                  v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + &
+                     &            v_ice(ji,jj-1) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + &
+                     &            v_oce(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
+               ELSE                             ! everywhere else
+                  v_ice(ji,jj) = v_oce(ji,jj)
+               ENDIF
+
+            ENDDO
+            
+         CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy )
+         CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy )
+ 
+         ELSEIF ( kn == 2 ) THEN ! set up of stress tensor (not sure it works)
+         
+            jgrd = 1      ! T grid
+            DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd)
+               ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
+               jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
+
+               pvar1 (ji,jj) = 0._wp
+               pvar2 (ji,jj) = 0._wp
+               pvar12(ji,jj) = 0._wp
+
+            ENDDO
+
+         ENDIF
+
+      CASE DEFAULT
+         CALL ctl_stop( 'bdy_ice_lim_dyn : unrecognised option for open boundaries for ice fields' )
+      END SELECT
+
+      ENDDO
+
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_lim_dyn')
+ 
+    END SUBROUTINE bdy_ice_lim_dyn
+
+#else
+   !!---------------------------------------------------------------------------------
+   !!   Default option                                                    Empty module
+   !!---------------------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE bdy_ice_lim( kt )      ! Empty routine
+      WRITE(*,*) 'bdy_ice_lim: You should not have seen this print! error?', kt
+   END SUBROUTINE bdy_ice_lim
+#endif
+
+   !!=================================================================================
+END MODULE bdyice_lim
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim2.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim2.F90	(revision 4267)
+++ 	(revision )
@@ -1,113 +1,0 @@
-MODULE bdyice_lim2
-   !!======================================================================
-   !!                       ***  MODULE  bdyice_lim2  ***
-   !! Unstructured Open Boundary Cond. :  Open boundary conditions for sea-ice (LIM2)
-   !!======================================================================
-   !!  History :  3.3  !  2010-09 (D. Storkey)  Original code
-   !!             3.4  !  2011    (D. Storkey) rewrite in preparation for OBC-BDY merge
-   !!             3.5  !  2012    (S. Mocavero, I. Epicoco) Optimization of BDY communications
-   !!----------------------------------------------------------------------
-#if defined   key_bdy   &&   defined key_lim2
-   !!----------------------------------------------------------------------
-   !!   'key_bdy'            and                 Unstructured Open Boundary Conditions
-   !!   'key_lim2'                                                 LIM-2 sea ice model
-   !!----------------------------------------------------------------------
-   !!   bdy_ice_lim_2      : Application of open boundaries to ice
-   !!   bdy_ice_frs        : Application of Flow Relaxation Scheme
-   !!----------------------------------------------------------------------
-   USE timing          ! Timing
-   USE oce             ! ocean dynamics and tracers variables
-   USE ice_2           ! LIM_2 ice variables
-   USE dom_oce         ! ocean space and time domain variables 
-   USE bdy_oce         ! ocean open boundary conditions
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE in_out_manager  ! write to numout file
-   USE lib_mpp         ! distributed memory computing
-   
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   bdy_ice_lim_2    ! routine called in sbcmod
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE bdy_ice_lim_2( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  SUBROUTINE bdy_ice_lim_2  ***
-      !!
-      !! ** Purpose : - Apply open boundary conditions for ice (LIM2)
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) :: kt     ! Main time step counter
-      !!
-      INTEGER               :: ib_bdy ! Loop index
-
-      DO ib_bdy=1, nb_bdy
-
-         SELECT CASE( cn_ice_lim2(ib_bdy) )
-         CASE('none')
-            CYCLE
-         CASE('frs')
-            CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
-         CASE DEFAULT
-            CALL ctl_stop( 'bdy_ice_lim_2 : unrecognised option for open boundaries for ice fields' )
-         END SELECT
-      ENDDO
-
-   END SUBROUTINE bdy_ice_lim_2
-
-   SUBROUTINE bdy_ice_frs( idx, dta, ib_bdy )
-      !!------------------------------------------------------------------------------
-      !!                 ***  SUBROUTINE bdy_ice_frs  ***
-      !!                    
-      !! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields in the case 
-      !!              of unstructured open boundaries. Currently only tested for LIM2.
-      !! 
-      !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three-
-      !!             dimensional baroclinic ocean model with realistic topography. Tellus, 365-382.
-      !!------------------------------------------------------------------------------
-      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
-      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
-      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
-      !!
-      INTEGER  ::   jb, jgrd   ! dummy loop indices
-      INTEGER  ::   ii, ij         ! local scalar
-      REAL(wp) ::   zwgt, zwgt1    ! local scalar
-      !!------------------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 ) CALL timing_start('bdy_ice_frs')
-      !
-      jgrd = 1      ! Everything is at T-points here
-      !
-      DO jb = 1, idx%nblen(jgrd)
-            ii    = idx%nbi(jb,jgrd)
-            ij    = idx%nbj(jb,jgrd)
-            zwgt  = idx%nbw(jb,jgrd)
-            zwgt1 = 1.e0 - idx%nbw(jb,jgrd)
-            frld (ii,ij) = ( frld (ii,ij) * zwgt1 + dta%frld (jb) * zwgt ) * tmask(ii,ij,1)     ! Leads fraction 
-            hicif(ii,ij) = ( hicif(ii,ij) * zwgt1 + dta%hicif(jb) * zwgt ) * tmask(ii,ij,1)     ! Ice depth 
-            hsnif(ii,ij) = ( hsnif(ii,ij) * zwgt1 + dta%hsnif(jb) * zwgt ) * tmask(ii,ij,1)     ! Snow depth
-      END DO 
-      CALL lbc_bdy_lnk( frld, 'T', 1., ib_bdy )                                         ! lateral boundary conditions
-      CALL lbc_bdy_lnk( hicif, 'T', 1., ib_bdy )   ;   CALL lbc_bdy_lnk( hsnif, 'T', 1., ib_bdy )
-      !      
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_frs')
-      !
-   END SUBROUTINE bdy_ice_frs
-#else
-   !!---------------------------------------------------------------------------------
-   !!   Default option                                                    Empty module
-   !!---------------------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE bdy_ice_lim_2( kt )      ! Empty routine
-      WRITE(*,*) 'bdy_ice_frs_lim_2: You should not have seen this print! error?', kt
-   END SUBROUTINE bdy_ice_lim_2
-#endif
-
-   !!=================================================================================
-END MODULE bdyice_lim2
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 4267)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 4278)
@@ -100,6 +100,6 @@
          &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &  
          &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
-#if defined key_lim2
-         &             cn_ice_lim2, nn_ice_lim2_dta,                           &
+#if ( defined key_lim2 || defined key_lim3 )
+         &             cn_ice_lim, nn_ice_lim_dta,                           &
 #endif
          &             ln_vol, nn_volctl, nn_rimwidth
@@ -141,7 +141,7 @@
       ln_dyn3d_dmp(:)   = .false.
       rn_time_dmp(:)    = 1.
-#if defined key_lim2
-      cn_ice_lim2(:)    = ''
-      nn_ice_lim2_dta(:)= -1  ! uninitialised flag
+#if ( defined key_lim2 || defined key_lim3 )
+      cn_ice_lim(:)    = ''
+      nn_ice_lim_dta(:)= -1  ! uninitialised flag
 #endif
       ln_vol            = .false.
@@ -333,5 +333,5 @@
 #if defined key_lim2
         IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
-        SELECT CASE( cn_ice_lim2(ib_bdy) )                  
+        SELECT CASE( cn_ice_lim(ib_bdy) )                  
           CASE('none')
              IF(lwp) WRITE(numout,*) '      no open boundary condition'        
@@ -344,11 +344,34 @@
              dta_bdy(ib_bdy)%ll_hicif = .true.
              dta_bdy(ib_bdy)%ll_hsnif = .true.
-          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice_lim2' )
+          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice_lim' )
         END SELECT
-        IF( cn_ice_lim2(ib_bdy) /= 'none' ) THEN 
-           SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   ! 
+        IF( cn_ice_lim(ib_bdy) /= 'none' ) THEN 
+           SELECT CASE( nn_ice_lim_dta(ib_bdy) )                   ! 
               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
               CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
-              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' )
+              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim_dta must be 0 or 1' )
+           END SELECT
+        ENDIF
+        IF(lwp) WRITE(numout,*)
+#elif defined key_lim3
+        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
+        SELECT CASE( cn_ice_lim(ib_bdy) )                  
+          CASE('none')
+             IF(lwp) WRITE(numout,*) '      no open boundary condition'        
+             dta_bdy(ib_bdy)%ll_a_i  = .false.
+             dta_bdy(ib_bdy)%ll_ht_i = .false.
+             dta_bdy(ib_bdy)%ll_ht_s = .false.
+          CASE('frs')
+             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+             dta_bdy(ib_bdy)%ll_a_i  = .true.
+             dta_bdy(ib_bdy)%ll_ht_i = .true.
+             dta_bdy(ib_bdy)%ll_ht_s = .true.
+          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice_lim' )
+        END SELECT
+        IF( cn_ice_lim(ib_bdy) /= 'none' ) THEN 
+           SELECT CASE( nn_ice_lim_dta(ib_bdy) )                   ! 
+              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
+              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
+              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim_dta must be 0 or 1' )
            END SELECT
         ENDIF
Index: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4267)
+++ branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4278)
@@ -89,5 +89,4 @@
       ! Update data, open boundaries, surface boundary condition (including sea-ice)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL sbc     ( kstp )        ! Sea Boundary Condition (including sea-ice)
       IF( lk_tide    )   CALL sbc_tide( kstp )
       IF( lk_obc     )   CALL obc_dta ( kstp )        ! update dynamic and tracer data at open boundaries
@@ -95,4 +94,6 @@
       IF( lk_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
 
+                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
+                                                      ! clem: moved here for bdy ice purpose
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       !  Ocean dynamics : hdiv, rot, ssh, e3, wn
