Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/closea.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/closea.F90	(revision 2113)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/closea.F90	(revision 2113)
@@ -0,0 +1,453 @@
+MODULE closea
+   !!======================================================================
+   !!                       ***  MODULE  closea  ***
+   !! Closed Seas  : specific treatments associated with closed seas
+   !!======================================================================
+   !! History :   8.2  !  00-05  (O. Marti)  Original code
+   !!             8.5  !  02-06  (E. Durand, G. Madec)  F90
+   !!             9.0  !  06-07  (G. Madec)  add clo_rnf, clo_ups, clo_bat
+   !!        NEMO 3.4  !  03-12  (P.G. Fogli) sbc_clo bug fix & mpp reproducibility
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dom_clo    : modification of the ocean domain for closed seas cases
+   !!   sbc_clo    : Special handling of closed seas
+   !!   clo_rnf    : set close sea outflows as river mouths (see sbcrnf)
+   !!   clo_ups    : set mixed centered/upstream scheme in closed sea (see traadv_cen2)
+   !!   clo_bat    : set to zero a field over closed sea (see domzrg)
+   !!----------------------------------------------------------------------
+   USE oce             ! dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE in_out_manager  ! I/O manager
+   USE sbc_oce         ! ocean surface boundary conditions
+   USE lib_fortran,    ONLY: glob_sum, DDPDD
+   USE lbclnk          ! lateral boundary condition - MPP exchanges
+   USE lib_mpp         ! MPP library
+   USE timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC dom_clo      ! routine called by domain module
+   PUBLIC sbc_clo      ! routine called by step module
+   PUBLIC clo_rnf      ! routine called by sbcrnf module
+   PUBLIC clo_ups      ! routine called in traadv_cen2(_jki) module
+   PUBLIC clo_bat      ! routine called in domzgr module
+
+   INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 4      !: number of closed sea
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncstt            !: Type of closed sea
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi1, ncsj1     !: south-west closed sea limits (i,j)
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi2, ncsj2     !: north-east closed sea limits (i,j)
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsnr            !: number of point where run-off pours
+   INTEGER, PUBLIC, DIMENSION(jpncs,4) ::   ncsir, ncsjr     !: Location of runoff
+
+   REAL(wp), DIMENSION (jpncs+1)       ::   surf             ! closed sea surface
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: closea.F90 3632 2012-11-22 15:28:42Z acc $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dom_clo
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_clo  ***
+      !!        
+      !! ** Purpose :   Closed sea domain initialization
+      !!
+      !! ** Method  :   if a closed sea is located only in a model grid point
+      !!                just the thermodynamic processes are applied.
+      !!
+      !! ** Action  :   ncsi1(), ncsj1() : south-west closed sea limits (i,j)
+      !!                ncsi2(), ncsj2() : north-east Closed sea limits (i,j)
+      !!                ncsir(), ncsjr() : Location of runoff
+      !!                ncsnr            : number of point where run-off pours
+      !!                ncstt            : Type of closed sea
+      !!                                   =0 spread over the world ocean
+      !!                                   =2 put at location runoff
+      !!----------------------------------------------------------------------
+      INTEGER ::   jc            ! dummy loop indices
+      !!----------------------------------------------------------------------
+      
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*)'dom_clo : closed seas '
+      IF(lwp) WRITE(numout,*)'~~~~~~~'
+
+      ! initial values
+      ncsnr(:) = 1  ;  ncsi1(:) = 1  ;  ncsi2(:) = 1  ;  ncsir(:,:) = 1
+      ncstt(:) = 0  ;  ncsj1(:) = 1  ;  ncsj2(:) = 1  ;  ncsjr(:,:) = 1
+
+      ! set the closed seas (in data domain indices)
+      ! -------------------
+
+      IF( cp_cfg == "orca" ) THEN
+         !
+         SELECT CASE ( jp_cfg )
+         !                                           ! =======================
+         CASE ( 1 )                                  ! ORCA_R1 configuration
+            !                                        ! =======================
+            ncsnr(1)   = 1    ; ncstt(1)   = 0           ! Caspian Sea
+            ncsir(1,1) = 1    ; ncsjr(1,1) = 1
+            ncsi1(1)   = 332  ; ncsj1(1)   = 203
+            ncsi2(1)   = 344  ; ncsj2(1)   = 235
+            !                                        
+            !                                        ! =======================
+         CASE ( 2 )                                  !  ORCA_R2 configuration
+            !                                        ! =======================
+            !                                            ! Caspian Sea
+            ncsnr(1)   =   1  ;  ncstt(1)   =   0           ! spread over the globe
+            ncsir(1,1) =   1  ;  ncsjr(1,1) =   1 
+            ncsi1(1)   =  11  ;  ncsj1(1)   = 103           !  
+            ncsi2(1)   =  17  ;  ncsj2(1)   = 112
+            !                                            ! Great North American Lakes
+            ncsnr(2)   =   1  ;  ncstt(2)   =   2           ! put at St Laurent mouth
+            ncsir(2,1) = 110  ;  ncsjr(2,1) = 111           
+            ncsi1(2)   =  97  ;  ncsj1(2)   = 107           !
+            ncsi2(2)   = 103  ;  ncsj2(2)   = 111
+            !                                            ! Black Sea (crossed by the cyclic boundary condition)
+            ncsnr(3:4) =   4  ;  ncstt(3:4) =   2           ! put in Med Sea (north of Aegean Sea)
+            ncsir(3:4,1) = 171;  ncsjr(3:4,1) = 106         !
+            ncsir(3:4,2) = 170;  ncsjr(3:4,2) = 106 
+            ncsir(3:4,3) = 171;  ncsjr(3:4,3) = 105 
+            ncsir(3:4,4) = 170;  ncsjr(3:4,4) = 105 
+            ncsi1(3)   = 174  ;  ncsj1(3)   = 107           ! 1 : west part of the Black Sea      
+            ncsi2(3)   = 181  ;  ncsj2(3)   = 112           !            (ie west of the cyclic b.c.)
+            ncsi1(4)   =   2  ;  ncsj1(4)   = 107           ! 2 : east part of the Black Sea 
+            ncsi2(4)   =   6  ;  ncsj2(4)   = 112           !           (ie east of the cyclic b.c.)
+             
+          
+
+            !                                        ! =======================
+         CASE ( 4 )                                  !  ORCA_R4 configuration
+            !                                        ! =======================
+            !                                            ! Caspian Sea
+            ncsnr(1)   =  1  ;  ncstt(1)   =  0  
+            ncsi1(1)   =  4  ;  ncsj1(1)   = 53 
+            ncsi2(1)   =  4  ;  ncsj2(1)   = 56
+            ncsir(1,1) =  1  ;  ncsjr(1,1) =  1
+            !                                            ! Great North American Lakes
+            ncsnr(2)   =  1  ;  ncstt(2)   =  2 
+            ncsi1(2)   = 49  ;  ncsj1(2)   = 55
+            ncsi2(2)   = 51  ;  ncsj2(2)   = 56
+            ncsir(2,1) = 57  ;  ncsjr(2,1) = 55
+            !                                            ! Black Sea
+            ncsnr(3)   =  4  ;  ncstt(3)   =  2  
+            ncsi1(3)   = 88  ;  ncsj1(3)   = 55 
+            ncsi2(3)   = 91  ;  ncsj2(3)   = 56
+            ncsir(3,1) = 86  ;  ncsjr(3,1) = 53
+            ncsir(3,2) = 87  ;  ncsjr(3,2) = 53 
+            ncsir(3,3) = 86  ;  ncsjr(3,3) = 52 
+            ncsir(3,4) = 87  ;  ncsjr(3,4) = 52
+            !                                            ! Baltic Sea
+            ncsnr(4)   =  1  ;  ncstt(4)   =  2
+            ncsi1(4)   = 75  ;  ncsj1(4)   = 59
+            ncsi2(4)   = 76  ;  ncsj2(4)   = 61
+            ncsir(4,1) = 84  ;  ncsjr(4,1) = 59 
+            !                                        ! =======================
+         CASE ( 025 )                                ! ORCA_R025 configuration
+            !                                        ! =======================
+            ncsnr(1)   = 1    ; ncstt(1)   = 0               ! Caspian + Aral sea
+            ncsi1(1)   = 1330 ; ncsj1(1)   = 645
+            ncsi2(1)   = 1400 ; ncsj2(1)   = 795
+            ncsir(1,1) = 1    ; ncsjr(1,1) = 1
+            !                                        
+            ncsnr(2)   = 1    ; ncstt(2)   = 0               ! Azov Sea 
+            ncsi1(2)   = 1284 ; ncsj1(2)   = 722
+            ncsi2(2)   = 1304 ; ncsj2(2)   = 747
+            ncsir(2,1) = 1    ; ncsjr(2,1) = 1
+            !
+         END SELECT
+         !
+      ENDIF
+
+      ! convert the position in local domain indices
+      ! --------------------------------------------
+      DO jc = 1, jpncs
+         ncsi1(jc)   = mi0( ncsi1(jc) )
+         ncsj1(jc)   = mj0( ncsj1(jc) )
+
+         ncsi2(jc)   = mi1( ncsi2(jc) )   
+         ncsj2(jc)   = mj1( ncsj2(jc) )  
+      END DO
+      !
+   END SUBROUTINE dom_clo
+
+
+   SUBROUTINE sbc_clo( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_clo  ***
+      !!                    
+      !! ** Purpose :   Special handling of closed seas
+      !!
+      !! ** Method  :   Water flux is forced to zero over closed sea
+      !!      Excess is shared between remaining ocean, or
+      !!      put as run-off in open ocean.
+      !!
+      !! ** Action  :   emp updated surface freshwater fluxes and associated heat content at kt
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean model time step
+      !
+      INTEGER             ::   ji, jj, jc, jn   ! dummy loop indices
+      REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon
+      REAL(wp)            ::   zze2, ztmp, zcorr     ! 
+      REAL(wp)            ::   zcoef, zcoef1         ! 
+      COMPLEX(wp)         ::   ctmp 
+      REAL(wp), DIMENSION(jpncs) ::   zfwf   ! 1D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('sbc_clo')
+      !                                                   !------------------!
+      IF( kt == nit000 ) THEN                             !  Initialisation  !
+         !                                                !------------------!
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*)'sbc_clo : closed seas '
+         IF(lwp) WRITE(numout,*)'~~~~~~~'
+
+         surf(:) = 0.e0_wp
+         !
+         surf(jpncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean
+         !
+         !                                        ! surface of closed seas 
+         IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation
+            DO jc = 1, jpncs
+               ctmp = CMPLX( 0.e0, 0.e0, wp )
+               DO jj = ncsj1(jc), ncsj2(jc)
+                  DO ji = ncsi1(jc), ncsi2(jc)
+                     ztmp = e1e2t(ji,jj) * tmask_i(ji,jj)
+                     CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+                  END DO 
+               END DO 
+               IF( lk_mpp )   CALL mpp_sum( ctmp )
+               surf(jc) = REAL(ctmp,wp)
+            END DO
+         ELSE                                          ! Standard calculation           
+            DO jc = 1, jpncs
+               DO jj = ncsj1(jc), ncsj2(jc)
+                  DO ji = ncsi1(jc), ncsi2(jc)
+                     surf(jc) = surf(jc) + e1e2t(ji,jj) * tmask_i(ji,jj)      ! surface of closed seas
+                  END DO 
+               END DO 
+            END DO 
+            IF( lk_mpp )   CALL mpp_sum ( surf, jpncs )       ! mpp: sum over all the global domain
+         ENDIF
+
+         IF(lwp) WRITE(numout,*)'     Closed sea surfaces'
+         DO jc = 1, jpncs
+            IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc)
+         END DO
+
+         ! jpncs+1 : surface of sea, closed seas excluded
+         DO jc = 1, jpncs
+            surf(jpncs+1) = surf(jpncs+1) - surf(jc)
+         END DO           
+         !
+      ENDIF
+      !
+      !                                                   !--------------------!
+      !                                                   !  update emp        !
+      zfwf = 0.e0_wp                                      !--------------------!
+      IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation
+         DO jc = 1, jpncs
+            ctmp = CMPLX( 0.e0, 0.e0, wp )
+            DO jj = ncsj1(jc), ncsj2(jc)
+               DO ji = ncsi1(jc), ncsi2(jc)
+                  ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)
+                  CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+               END DO  
+            END DO 
+            IF( lk_mpp )   CALL mpp_sum( ctmp )
+            zfwf(jc) = REAL(ctmp,wp)
+         END DO
+      ELSE                                          ! Standard calculation           
+         DO jc = 1, jpncs
+            DO jj = ncsj1(jc), ncsj2(jc)
+               DO ji = ncsi1(jc), ncsi2(jc)
+                  zfwf(jc) = zfwf(jc) + e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 
+               END DO  
+            END DO 
+         END DO
+         IF( lk_mpp )   CALL mpp_sum ( zfwf(:) , jpncs )       ! mpp: sum over all the global domain
+      ENDIF
+
+      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration
+         zze2    = ( zfwf(3) + zfwf(4) ) * 0.5_wp
+         zfwf(3) = zze2
+         zfwf(4) = zze2
+      ENDIF
+
+      zcorr = 0._wp
+
+      DO jc = 1, jpncs
+         !
+         ! The following if avoids the redistribution of the round off
+         IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN
+            !
+            IF( ncstt(jc) == 0 ) THEN           ! water/evap excess is shared by all open ocean
+               zcoef    = zfwf(jc) / surf(jpncs+1)
+               zcoef1   = rcp * zcoef
+               emp(:,:) = emp(:,:) + zcoef
+               qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
+               ! accumulate closed seas correction
+               zcorr    = zcorr    + zcoef
+               !
+            ELSEIF( ncstt(jc) == 1 ) THEN       ! Excess water in open sea, at outflow location, excess evap shared
+               IF ( zfwf(jc) <= 0.e0_wp ) THEN 
+                   DO jn = 1, ncsnr(jc)
+                     ji = mi0(ncsir(jc,jn))
+                     jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
+                     IF (      ji > 1 .AND. ji < jpi   &
+                         .AND. jj > 1 .AND. jj < jpj ) THEN 
+                         zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) )
+                         zcoef1     = rcp * zcoef
+                         emp(ji,jj) = emp(ji,jj) + zcoef
+                         qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
+                     ENDIF 
+                   END DO 
+               ELSE 
+                   zcoef    = zfwf(jc) / surf(jpncs+1)
+                   zcoef1   = rcp * zcoef
+                   emp(:,:) = emp(:,:) + zcoef
+                   qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
+                   ! accumulate closed seas correction
+                   zcorr    = zcorr    + zcoef
+               ENDIF
+            ELSEIF( ncstt(jc) == 2 ) THEN       ! Excess e-p-r (either sign) goes to open ocean, at outflow location
+               DO jn = 1, ncsnr(jc)
+                  ji = mi0(ncsir(jc,jn))
+                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
+                  IF(      ji > 1 .AND. ji < jpi    &
+                     .AND. jj > 1 .AND. jj < jpj ) THEN 
+                     zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) *  e1e2t(ji,jj) )
+                     zcoef1     = rcp * zcoef
+                     emp(ji,jj) = emp(ji,jj) + zcoef
+                     qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
+                  ENDIF 
+               END DO 
+            ENDIF 
+            !
+            DO jj = ncsj1(jc), ncsj2(jc)
+               DO ji = ncsi1(jc), ncsi2(jc)
+                  zcoef      = zfwf(jc) / surf(jc)
+                  zcoef1     = rcp * zcoef
+                  emp(ji,jj) = emp(ji,jj) - zcoef
+                  qns(ji,jj) = qns(ji,jj) + zcoef1 * sst_m(ji,jj)
+               END DO  
+            END DO 
+            !
+         END IF
+      END DO 
+
+      IF ( ABS(zcorr) > rsmall ) THEN      ! remove the global correction from the closed seas
+         DO jc = 1, jpncs                  ! only if it is large enough
+            DO jj = ncsj1(jc), ncsj2(jc)
+               DO ji = ncsi1(jc), ncsi2(jc)
+                  emp(ji,jj) = emp(ji,jj) - zcorr
+                  qns(ji,jj) = qns(ji,jj) + rcp * zcorr * sst_m(ji,jj)
+               END DO  
+             END DO 
+          END DO
+      ENDIF
+      !
+      emp (:,:) = emp (:,:) * tmask(:,:,1)
+      !
+      CALL lbc_lnk( emp , 'T', 1._wp )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('sbc_clo')
+      !
+   END SUBROUTINE sbc_clo
+
+
+   SUBROUTINE clo_rnf( p_rnfmsk )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_rnf  ***
+      !!                    
+      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
+      !!                to be the same as river mouth grid-points
+      !!
+      !! ** Method  :   set to 1 the runoff mask (mskrnf, see sbcrnf module)
+      !!                at the closed sea outflow grid-point.
+      !!
+      !! ** Action  :   update (p_)mskrnf (set 1 at closed sea outflow)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array)
+      !
+      INTEGER  ::   jc, jn      ! dummy loop indices
+      INTEGER  ::   ii, ij      ! temporary integer
+      !!----------------------------------------------------------------------
+      !
+      DO jc = 1, jpncs
+         IF( ncstt(jc) >= 1 ) THEN            ! runoff mask set to 1 at closed sea outflows
+             DO jn = 1, 4
+               ii = mi0( ncsir(jc,jn) )
+               ij = mj0( ncsjr(jc,jn) )
+               p_rnfmsk(ii,ij) = MAX( p_rnfmsk(ii,ij), 1.0_wp )
+            END DO 
+         ENDIF 
+      END DO 
+      !
+   END SUBROUTINE clo_rnf
+
+   
+   SUBROUTINE clo_ups( p_upsmsk )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_rnf  ***
+      !!                    
+      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
+      !!                to be the same as river mouth grid-points
+      !!
+      !! ** Method  :   set to 0.5 the upstream mask (upsmsk, see traadv_cen2 
+      !!                module) over the closed seas.
+      !!
+      !! ** Action  :   update (p_)upsmsk (set 0.5 over closed seas)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_upsmsk   ! upstream mask (upsmsk array)
+      !
+      INTEGER  ::   jc, ji, jj      ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      DO jc = 1, jpncs
+         DO jj = ncsj1(jc), ncsj2(jc)
+            DO ji = ncsi1(jc), ncsi2(jc)
+               p_upsmsk(ji,jj) = 0.5_wp         ! mixed upstream/centered scheme over closed seas
+            END DO 
+         END DO 
+       END DO 
+       !
+   END SUBROUTINE clo_ups
+   
+      
+   SUBROUTINE clo_bat( pbat, kbat )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE clo_bat  ***
+      !!                    
+      !! ** Purpose :   suppress closed sea from the domain
+      !!
+      !! ** Method  :   set to 0 the meter and level bathymetry (given in 
+      !!                arguments) over the closed seas.
+      !!
+      !! ** Action  :   set pbat=0 and kbat=0 over closed seas
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pbat   ! bathymetry in meters (bathy array)
+      INTEGER , DIMENSION(jpi,jpj), INTENT(inout) ::   kbat   ! bathymetry in levels (mbathy array)
+      !
+      INTEGER  ::   jc, ji, jj      ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      DO jc = 1, jpncs
+         DO jj = ncsj1(jc), ncsj2(jc)
+            DO ji = ncsi1(jc), ncsi2(jc)
+               pbat(ji,jj) = 0._wp   
+               kbat(ji,jj) = 0   
+            END DO 
+         END DO 
+       END DO 
+       !
+   END SUBROUTINE clo_bat
+
+   !!======================================================================
+END MODULE closea
+
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limistate.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limistate.F90	(revision 2113)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limistate.F90	(revision 2113)
@@ -0,0 +1,587 @@
+MODULE limistate
+   !!======================================================================
+   !!                     ***  MODULE  limistate  ***
+   !!              Initialisation of diagnostics ice variables
+   !!======================================================================
+   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
+   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3' :                                    LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   lim_istate      :  Initialisation of diagnostics ice variables
+   !!   lim_istate_init :  initialization of ice state and namelist read
+   !!----------------------------------------------------------------------
+   USE phycst           ! physical constant
+   USE oce              ! dynamics and tracers variables
+   USE dom_oce          ! ocean domain
+   USE sbc_oce          ! Surface boundary condition: ocean fields
+   USE sbc_ice          ! Surface boundary condition: ice fields
+   USE eosbn2           ! equation of state
+   USE ice              ! sea-ice variables
+   USE par_ice          ! ice parameters
+   USE dom_ice          ! sea-ice domain
+   USE in_out_manager   ! I/O manager
+   USE lbclnk           ! lateral boundary condition - MPP exchanges
+   USE lib_mpp          ! MPP library
+   USE wrk_nemo         ! work arrays
+   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   lim_istate      ! routine called by lim_init.F90
+
+   !                                  !!** init namelist (namiceini) **
+   REAL(wp) ::   ttest    = 2.0_wp     ! threshold water temperature for initial sea ice
+   REAL(wp) ::   hninn    = 0.5_wp     ! initial snow thickness in the north
+   REAL(wp) ::   hginn_u  = 2.5_wp     ! initial ice thickness in the north
+   REAL(wp) ::   aginn_u  = 0.7_wp     ! initial leads area in the north
+   REAL(wp) ::   hginn_d  = 5.0_wp     ! initial ice thickness in the north
+   REAL(wp) ::   aginn_d  = 0.25_wp    ! initial leads area in the north
+   REAL(wp) ::   hnins    = 0.1_wp     ! initial snow thickness in the south
+   REAL(wp) ::   hgins_u  = 1.0_wp     ! initial ice thickness in the south
+   REAL(wp) ::   agins_u  = 0.7_wp     ! initial leads area in the south
+   REAL(wp) ::   hgins_d  = 2.0_wp     ! initial ice thickness in the south
+   REAL(wp) ::   agins_d  = 0.2_wp     ! initial leads area in the south
+   REAL(wp) ::   sinn     = 6.301_wp   ! initial salinity 
+   REAL(wp) ::   sins     = 6.301_wp   !
+
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)
+   !! $Id: limistate.F90 3625 2012-11-21 13:19:18Z acc $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE lim_istate
+      !!-------------------------------------------------------------------
+      !!                    ***  ROUTINE lim_istate  ***
+      !!
+      !! ** Purpose :   defined the sea-ice initial state
+      !!
+      !! ** Method  :   restart from a state defined in a binary file
+      !!                or from arbitrary sea-ice conditions
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
+      REAL(wp) ::   zeps6, zeps, ztmelts, epsi06   ! local scalars
+      REAL(wp) ::   zvol, zare, zh, zh1, zh2, zh3, zan, zbn, zas, zbs 
+      REAL(wp), POINTER, DIMENSION(:)   ::   zgfactorn, zhin 
+      REAL(wp), POINTER, DIMENSION(:)   ::   zgfactors, zhis
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zidto      ! ice indicator
+      !--------------------------------------------------------------------
+
+      CALL wrk_alloc( jpm, zgfactorn, zgfactors, zhin, zhis )
+      CALL wrk_alloc( jpi, jpj, zidto )
+
+      !--------------------------------------------------------------------
+      ! 1) Preliminary things 
+      !--------------------------------------------------------------------
+      epsi06 = 1.e-6_wp
+
+      CALL lim_istate_init     !  reading the initials parameters of the ice
+
+!!gm  in lim2  the initialisation if only done if required in the namelist :
+!!gm      IF( .NOT. ln_limini ) THEN
+!!gm  this should be added in lim3 namelist...
+
+      !--------------------------------------------------------------------
+      ! 2) Ice initialization (hi,hs,frld,t_su,sm_i,t_i,t_s)              | 
+      !--------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
+
+      t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius]
+
+      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
+         DO ji = 1, jpi
+            IF( tsn(ji,jj,1,jp_tem)  - t_bo(ji,jj) >= ttest ) THEN   ;   zidto(ji,jj) = 0.e0      ! no ice
+            ELSE                                                     ;   zidto(ji,jj) = 1.e0      !    ice
+            ENDIF
+         END DO
+      END DO
+
+      t_bo(:,:) = t_bo(:,:) + rt0                          ! t_bo converted from Celsius to Kelvin (rt0 over land)
+
+      ! constants for heat contents
+      zeps   = 1.e-20_wp
+      zeps6  = 1.e-06_wp
+
+      ! zgfactor for initial ice distribution
+      zgfactorn(:) = 0._wp
+      zgfactors(:) = 0._wp
+
+      ! first ice type
+      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
+         zhin (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
+         zgfactorn(1) = zgfactorn(1) + exp(-(zhin(1)-hginn_u)*(zhin(1)-hginn_u) * 0.5_wp )
+         zhis (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
+         zgfactors(1) = zgfactors(1) + exp(-(zhis(1)-hgins_u)*(zhis(1)-hgins_u) * 0.5_wp )
+      END DO ! jl
+      zgfactorn(1) = aginn_u / zgfactorn(1)
+      zgfactors(1) = agins_u / zgfactors(1)
+
+      ! -------------
+      ! new distribution, polynom of second order, conserving area and volume
+      zh1 = 0._wp
+      zh2 = 0._wp
+      zh3 = 0._wp
+      DO jl = 1, jpl
+         zh = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
+         zh1 = zh1 + zh
+         zh2 = zh2 + zh * zh
+         zh3 = zh3 + zh * zh * zh
+      END DO
+      IF(lwp) WRITE(numout,*) ' zh1 : ', zh1
+      IF(lwp) WRITE(numout,*) ' zh2 : ', zh2
+      IF(lwp) WRITE(numout,*) ' zh3 : ', zh3
+
+      zvol = aginn_u * hginn_u
+      zare = aginn_u
+      IF( jpl >= 2 ) THEN
+         zbn = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3)
+         zan = ( zare - zbn*zh1 ) / zh2
+      ENDIF
+
+      IF(lwp) WRITE(numout,*) ' zvol: ', zvol
+      IF(lwp) WRITE(numout,*) ' zare: ', zare
+      IF(lwp) WRITE(numout,*) ' zbn : ', zbn 
+      IF(lwp) WRITE(numout,*) ' zan : ', zan 
+
+      zvol = agins_u * hgins_u
+      zare = agins_u
+      IF( jpl >= 2 ) THEN
+         zbs = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3)
+         zas = ( zare - zbs*zh1 ) / zh2
+      ENDIF
+
+      IF(lwp) WRITE(numout,*) ' zvol: ', zvol
+      IF(lwp) WRITE(numout,*) ' zare: ', zare
+      IF(lwp) WRITE(numout,*) ' zbn : ', zbn 
+      IF(lwp) WRITE(numout,*) ' zan : ', zan 
+
+      !end of new lines
+      ! -------------
+!!!
+      ! retour a LIMA_MEC
+      !     ! second ice type
+      !     zdummy  = hi_max(ice_cat_bounds(2,1)-1)
+      !     hi_max(ice_cat_bounds(2,1)-1) = 0.0
+
+      !     ! here to change !!!!
+      !     jm = 2
+      !     DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
+      !        zhin (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
+      !        zhin (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + &
+      !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0
+      !        zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0)
+      !        zhis (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
+      !        zhis (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + &
+      !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0
+      !        zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0)
+      !     END DO ! jl
+      !     zgfactorn(2) = aginn_d / zgfactorn(2)
+      !     zgfactors(2) = agins_d / zgfactors(2)
+      !     hi_max(ice_cat_bounds(2,1)-1) = zdummy
+      ! END retour a LIMA_MEC
+!!!
+
+!!gm  optimisation :  loop over the ice categories inside the ji, jj loop !!!
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            !--- Northern hemisphere
+            !----------------------------------------------------------------
+            IF( fcor(ji,jj) >= 0._wp ) THEN    
+
+               !-----------------------
+               ! Ice area / thickness
+               !-----------------------
+
+               IF ( jpl .EQ. 1) THEN ! one category
+
+                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
+                     a_i(ji,jj,jl)    = zidto(ji,jj) * aginn_u
+                     ht_i(ji,jj,jl)   = zidto(ji,jj) * hginn_u
+                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
+                  END DO
+
+               ELSE ! several categories
+
+                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
+                     zhin(1)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
+                     a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(1) * exp(-(zhin(1)-hginn_u)* & 
+                        (zhin(1)-hginn_u)/2.0) , epsi06)
+                     ! new line
+                     a_i(ji,jj,jl)    = zidto(ji,jj) * ( zan * zhin(1) * zhin(1) + zbn * zhin(1) )
+                     ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(1) 
+                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
+                  END DO
+
+               ENDIF
+
+
+!!!
+               ! retour a LIMA_MEC
+               !              !ridged ice
+               !              zdummy  = hi_max(ice_cat_bounds(2,1)-1)
+               !              hi_max(ice_cat_bounds(2,1)-1) = 0.0
+               !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories
+               !                 zhin(2)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
+               !                 a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* &
+               !                                    (zhin(2)-hginn_d)/2.0) , epsi06)
+               !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(2) 
+               !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
+               !              END DO
+               !              hi_max(ice_cat_bounds(2,1)-1) = zdummy
+
+               !              !rafted ice
+               !              jl = 6
+               !              a_i(ji,jj,jl)       = 0.0
+               !              ht_i(ji,jj,jl)      = 0.0
+               !              v_i(ji,jj,jl)       = 0.0
+               ! END retour a LIMA_MEC
+!!!
+
+               DO jl = 1, jpl
+
+                  !-------------
+                  ! Snow depth
+                  !-------------
+                  ht_s(ji,jj,jl)   = zidto(ji,jj) * hninn
+                  v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl)
+
+                  !---------------
+                  ! Ice salinity
+                  !---------------
+                  sm_i(ji,jj,jl)   = zidto(ji,jj) * sinn  + ( 1.0 - zidto(ji,jj) ) * 0.1
+                  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)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) )
+                  oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl)
+
+                  !------------------------------
+                  ! Sea ice surface temperature
+                  !------------------------------
+
+                  t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj)
+
+                  !------------------------------------
+                  ! Snow temperature and heat content
+                  !------------------------------------
+
+                  DO jk = 1, nlay_s
+                     t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt
+                     ! Snow energy of melting
+                     e_s(ji,jj,jk,jl) = zidto(ji,jj) * 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) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
+                     s_i(ji,jj,jk,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1
+                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
+
+                     ! heat content per unit volume
+                     e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * &
+                        (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
+                        +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &
+                        - 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 ! jl 
+
+            ELSE ! on fcor 
+
+               !--- Southern hemisphere
+               !----------------------------------------------------------------
+
+               !-----------------------
+               ! Ice area / thickness
+               !-----------------------
+
+               IF ( jpl .EQ. 1) THEN ! one category
+
+                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
+                     a_i(ji,jj,jl)    = zidto(ji,jj) * agins_u
+                     ht_i(ji,jj,jl)   = zidto(ji,jj) * hgins_u
+                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
+                  END DO
+
+               ELSE ! several categories
+
+                  !level ice
+                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !over thickness categories
+
+                     zhis(1)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
+                     a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactors(1) * exp(-(zhis(1)-hgins_u) * & 
+                        (zhis(1)-hgins_u)/2.0) , epsi06 )
+                     ! new line square distribution volume conserving
+                     a_i(ji,jj,jl)    = zidto(ji,jj) * ( zas * zhis(1) * zhis(1) + zbs * zhis(1) )
+                     ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(1) 
+                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
+
+                  END DO ! jl
+
+               ENDIF
+
+!!!
+               ! retour a LIMA_MEC
+               !              !ridged ice
+               !              zdummy  = hi_max(ice_cat_bounds(2,1)-1)
+               !              hi_max(ice_cat_bounds(2,1)-1) = 0.0
+               !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories
+               !                 zhis(2)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
+               !                 a_i(ji,jj,jl) = zidto(ji,jj)*MAX( zgfactors(2)   &
+               !                    &          * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 )
+               !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(2) 
+               !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
+               !              END DO
+               !              hi_max(ice_cat_bounds(2,1)-1) = zdummy
+
+               !              !rafted ice
+               !              jl = 6
+               !              a_i(ji,jj,jl)       = 0.0
+               !              ht_i(ji,jj,jl)      = 0.0
+               !              v_i(ji,jj,jl)       = 0.0
+               ! END retour a LIMA_MEC
+!!!
+
+               DO jl = 1, jpl !over thickness categories
+
+                  !---------------
+                  ! Snow depth
+                  !---------------
+
+                  ht_s(ji,jj,jl)   = zidto(ji,jj) * hnins
+                  v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl)
+
+                  !---------------
+                  ! Ice salinity
+                  !---------------
+
+                  sm_i(ji,jj,jl)   = zidto(ji,jj) * sins  + ( 1.0 - zidto(ji,jj) ) * 0.1
+                  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)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) )
+                  oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl)
+
+                  !------------------------------
+                  ! Sea ice surface temperature
+                  !------------------------------
+
+                  t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj)
+
+                  !----------------------------------
+                  ! Snow temperature / heat content
+                  !----------------------------------
+
+                  DO jk = 1, nlay_s
+                     t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt
+                     ! Snow energy of melting
+                     e_s(ji,jj,jk,jl) = zidto(ji,jj) * 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
+
+                  !---------------------------------------------
+                  ! Ice temperature, salinity and heat content
+                  !---------------------------------------------
+
+                  DO jk = 1, nlay_i
+                     t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
+                     s_i(ji,jj,jk,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1
+                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
+
+                     ! heat content per unit volume
+                     e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * &
+                        (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
+                        +   lfus  * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &
+                        - 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 ! jl 
+
+            ENDIF ! on fcor
+
+         END DO
+      END DO
+
+      !--------------------------------------------------------------------
+      ! 3) Global ice variables for output diagnostics                    | 
+      !--------------------------------------------------------------------
+
+      fsbbq (:,:)     = 0.e0
+      u_ice (:,:)     = 0.e0
+      v_ice (:,:)     = 0.e0
+      stress1_i(:,:)  = 0.0
+      stress2_i(:,:)  = 0.0
+      stress12_i(:,:) = 0.0
+
+      !--------------------------------------------------------------------
+      ! 4) Moments for advection
+      !--------------------------------------------------------------------
+
+      sxopw (:,:) = 0.e0 
+      syopw (:,:) = 0.e0 
+      sxxopw(:,:) = 0.e0 
+      syyopw(:,:) = 0.e0 
+      sxyopw(:,:) = 0.e0
+
+      sxice (:,:,:)  = 0.e0   ;   sxsn (:,:,:)  = 0.e0   ;   sxa  (:,:,:)  = 0.e0
+      syice (:,:,:)  = 0.e0   ;   sysn (:,:,:)  = 0.e0   ;   sya  (:,:,:)  = 0.e0
+      sxxice(:,:,:)  = 0.e0   ;   sxxsn(:,:,:)  = 0.e0   ;   sxxa (:,:,:)  = 0.e0
+      syyice(:,:,:)  = 0.e0   ;   syysn(:,:,:)  = 0.e0   ;   syya (:,:,:)  = 0.e0
+      sxyice(:,:,:)  = 0.e0   ;   sxysn(:,:,:)  = 0.e0   ;   sxya (:,:,:)  = 0.e0
+
+      sxc0  (:,:,:)  = 0.e0   ;   sxe  (:,:,:,:)= 0.e0   
+      syc0  (:,:,:)  = 0.e0   ;   sye  (:,:,:,:)= 0.e0   
+      sxxc0 (:,:,:)  = 0.e0   ;   sxxe (:,:,:,:)= 0.e0   
+      syyc0 (:,:,:)  = 0.e0   ;   syye (:,:,:,:)= 0.e0   
+      sxyc0 (:,:,:)  = 0.e0   ;   sxye (:,:,:,:)= 0.e0   
+
+      sxsal  (:,:,:)  = 0.e0
+      sysal  (:,:,:)  = 0.e0
+      sxxsal (:,:,:)  = 0.e0
+      syysal (:,:,:)  = 0.e0
+      sxysal (:,:,:)  = 0.e0
+
+      !--------------------------------------------------------------------
+      ! 5) Lateral boundary conditions                                    | 
+      !--------------------------------------------------------------------
+
+      DO jl = 1, jpl
+         CALL lbc_lnk( a_i(:,:,jl)  , 'T', 1. )
+         CALL lbc_lnk( v_i(:,:,jl)  , 'T', 1. )
+         CALL lbc_lnk( v_s(:,:,jl)  , 'T', 1. )
+         CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. )
+         CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. )
+         !
+         CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. )
+         CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. )
+         CALL lbc_lnk( sm_i(:,:,jl) , 'T', 1. )
+         CALL lbc_lnk( o_i(:,:,jl)  , 'T', 1. )
+         CALL lbc_lnk( t_su(:,:,jl) , 'T', 1. )
+         DO jk = 1, nlay_s
+            CALL lbc_lnk(t_s(:,:,jk,jl), 'T', 1. )
+            CALL lbc_lnk(e_s(:,:,jk,jl), 'T', 1. )
+         END DO
+         DO jk = 1, nlay_i
+            CALL lbc_lnk(t_i(:,:,jk,jl), 'T', 1. )
+            CALL lbc_lnk(e_i(:,:,jk,jl), 'T', 1. )
+         END DO
+         !
+         a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl)
+      END DO
+      
+      at_i (:,:) = 0.0_wp
+      DO jl = 1, jpl
+         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
+      END DO
+
+      CALL lbc_lnk( at_i , 'T', 1. )
+      at_i(:,:) = tms(:,:) * at_i(:,:)                       ! put 0 over land
+      !
+      CALL lbc_lnk( fsbbq  , 'T', 1. )
+      !
+      CALL wrk_dealloc( jpm, zgfactorn, zgfactors, zhin, zhis )
+      CALL wrk_dealloc( jpi, jpj, zidto )
+      !
+      !--------------------------------------------------------------------
+      ! 6) ????                                                           | 
+      !--------------------------------------------------------------------
+      tn_ice (:,:,:) = t_su (:,:,:)
+      !
+   END SUBROUTINE lim_istate
+
+
+   SUBROUTINE lim_istate_init
+      !!-------------------------------------------------------------------
+      !!                   ***  ROUTINE lim_istate_init  ***
+      !!        
+      !! ** Purpose : Definition of initial state of the ice 
+      !!
+      !! ** Method :   Read the namiceini namelist and check the parameter 
+      !!             values called at the first timestep (nit000)
+      !!
+      !! ** input  :   namelist namiceini
+      !!-----------------------------------------------------------------------------
+      NAMELIST/namiceini/ ttest, hninn, hginn_u, aginn_u, hginn_d, aginn_d, hnins,   &
+         &                hgins_u, agins_u, hgins_d, agins_d, sinn, sins
+      !!-----------------------------------------------------------------------------
+      !
+      REWIND ( numnam_ice )               ! Read Namelist namiceini 
+      READ   ( numnam_ice , namiceini )
+      !
+      IF(lwp) THEN                        ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   threshold water temp. for initial sea-ice    ttest      = ', ttest
+         WRITE(numout,*) '   initial snow thickness in the north          hninn      = ', hninn
+         WRITE(numout,*) '   initial undef ice thickness in the north     hginn_u    = ', hginn_u
+         WRITE(numout,*) '   initial undef ice concentr. in the north     aginn_u    = ', aginn_u          
+         WRITE(numout,*) '   initial  def  ice thickness in the north     hginn_d    = ', hginn_d
+         WRITE(numout,*) '   initial  def  ice concentr. in the north     aginn_d    = ', aginn_d          
+         WRITE(numout,*) '   initial snow thickness in the south          hnins      = ', hnins 
+         WRITE(numout,*) '   initial undef ice thickness in the north     hgins_u    = ', hgins_u
+         WRITE(numout,*) '   initial undef ice concentr. in the north     agins_u    = ', agins_u          
+         WRITE(numout,*) '   initial  def  ice thickness in the north     hgins_d    = ', hgins_d
+         WRITE(numout,*) '   initial  def  ice concentr. in the north     agins_d    = ', agins_d          
+         WRITE(numout,*) '   initial  ice salinity       in the north     sinn       = ', sinn
+         WRITE(numout,*) '   initial  ice salinity       in the south     sins       = ', sins
+      ENDIF
+      !
+   END SUBROUTINE lim_istate_init
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :         Empty module          NO LIM sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE lim_istate          ! Empty routine
+   END SUBROUTINE lim_istate
+#endif
+
+   !!======================================================================
+END MODULE limistate
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limrst.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limrst.F90	(revision 2113)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limrst.F90	(revision 2113)
@@ -0,0 +1,585 @@
+MODULE limrst
+   !!======================================================================
+   !!                     ***  MODULE  limrst  ***
+   !! Ice restart :  write the ice restart file
+   !!======================================================================
+   !! History:   -   ! 2005-04 (M. Vancoppenolle) Original code
+   !!           3.0  ! 2008-03 (C. Ethe) restart files in using IOM interface
+   !!           4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3' :                                   LIM sea-ice model
+   !!----------------------------------------------------------------------
+   !!   lim_rst_opn   : open ice restart file
+   !!   lim_rst_write : write of the restart file 
+   !!   lim_rst_read  : read  the restart file 
+   !!----------------------------------------------------------------------
+   USE ice            ! sea-ice variables
+   USE par_ice        ! sea-ice parameters
+   USE dom_oce        ! ocean domain
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE sbc_ice        ! Surface boundary condition: ice fields
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! work arrays
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   lim_rst_opn    ! routine called by icestep.F90
+   PUBLIC   lim_rst_write  ! routine called by icestep.F90
+   PUBLIC   lim_rst_read   ! routine called by iceini.F90
+
+   LOGICAL, PUBLIC ::   lrst_ice         !: logical to control the ice restart write 
+   INTEGER, PUBLIC ::   numrir, numriw   !: logical unit for ice restart (read and write)
+
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)
+   !! $Id: limrst.F90 3625 2012-11-21 13:19:18Z acc $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE lim_rst_opn( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  lim_rst_opn  ***
+      !!
+      !! ** purpose  :   output of sea-ice variable in a netcdf file
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt       ! number of iteration
+      !
+      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
+      CHARACTER(LEN=50)   ::   clname   ! ice output restart file name
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 )   lrst_ice = .FALSE.   ! default definition
+
+      ! in order to get better performances with NetCDF format, we open and define the ice restart file 
+      ! one ice time step before writing the data (-> at nitrst - 2*nn_fsbc + 1), except if we write ice 
+      ! restart files every ice time step or if an ice restart file was writen at nitend - 2*nn_fsbc + 1
+      IF( kt == nitrst - 2*nn_fsbc + 1 .OR. nstock == nn_fsbc    &
+         &                             .OR. ( kt == nitend - nn_fsbc + 1 .AND. .NOT. lrst_ice ) ) THEN
+         ! beware of the format used to write kt (default is i8.8, that should be large enough...)
+         IF( nitrst > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
+         ELSE                           ;   WRITE(clkt, '(i8.8)') nitrst
+         ENDIF
+         ! create the file
+         clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_icerst_out)
+         IF(lwp) THEN
+            WRITE(numout,*)
+            SELECT CASE ( jprstlib )
+            CASE ( jprstdimg )   ;   WRITE(numout,*) '             open ice restart binary file: '//clname
+            CASE DEFAULT         ;   WRITE(numout,*) '             open ice restart NetCDF file: '//clname
+            END SELECT
+            IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN   
+               WRITE(numout,*)         '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp
+            ELSE   ;   WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp
+            ENDIF
+         ENDIF
+         !
+         CALL iom_open( clname, numriw, ldwrt = .TRUE., kiolib = jprstlib )
+         lrst_ice = .TRUE.
+      ENDIF
+      !
+   END SUBROUTINE lim_rst_opn
+
+
+   SUBROUTINE lim_rst_write( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  lim_rst_write  ***
+      !!
+      !! ** purpose  :   output of sea-ice variable in a netcdf file
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt     ! number of iteration
+      !!
+      INTEGER ::   ji, jj, jk ,jl   ! dummy loop indices
+      INTEGER ::   iter
+      CHARACTER(len=15) ::   znam
+      CHARACTER(len=1)  ::   zchar, zchar1
+      REAL(wp), POINTER, DIMENSION(:,:) :: z2d
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi, jpj, z2d )
+
+      iter = kt + nn_fsbc - 1   ! ice restarts are written at kt == nitrst - nn_fsbc + 1
+
+      IF( iter == nitrst ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'lim_rst_write : write ice restart file  kt =', kt
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'         
+      ENDIF
+
+      ! Write in numriw (if iter == nitrst)
+      ! ------------------ 
+      !                                                                        ! calendar control
+      CALL iom_rstput( iter, nitrst, numriw, 'nn_fsbc', REAL( nn_fsbc, wp ) )      ! time-step 
+      CALL iom_rstput( iter, nitrst, numriw, 'kt_ice' , REAL( iter   , wp ) )      ! date
+
+      ! Prognostic variables 
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I1)') jl
+         znam = 'v_i'//'_htc'//zchar
+         z2d(:,:) = v_i(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'v_s'//'_htc'//zchar
+         z2d(:,:) = v_s(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'smv_i'//'_htc'//zchar
+         z2d(:,:) = smv_i(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'oa_i'//'_htc'//zchar
+         z2d(:,:) = oa_i(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'a_i'//'_htc'//zchar
+         z2d(:,:) = a_i(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 't_su'//'_htc'//zchar
+         z2d(:,:) = t_su(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+      END DO
+
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I1)') jl
+         znam = 'tempt_sl1'//'_htc'//zchar
+         z2d(:,:) = e_s(:,:,1,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+      END DO
+
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I1)') jl
+         DO jk = 1, nlay_i 
+            WRITE(zchar1,'(I1)') jk
+            znam = 'tempt'//'_il'//zchar1//'_htc'//zchar
+            z2d(:,:) = e_i(:,:,jk,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         END DO
+      END DO
+
+      CALL iom_rstput( iter, nitrst, numriw, 'u_ice'     , u_ice      )
+      CALL iom_rstput( iter, nitrst, numriw, 'v_ice'     , v_ice      )
+      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq'     , fsbbq      )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
+
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I1)') jl
+         znam = 'sxice'//'_htc'//zchar
+         z2d(:,:) = sxice(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'syice'//'_htc'//zchar
+         z2d(:,:) = syice(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxxice'//'_htc'//zchar
+         z2d(:,:) = sxxice(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'syyice'//'_htc'//zchar
+         z2d(:,:) = syyice(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxyice'//'_htc'//zchar
+         z2d(:,:) = sxyice(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxsn'//'_htc'//zchar
+         z2d(:,:) = sxsn(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sysn'//'_htc'//zchar
+         z2d(:,:) = sysn(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxxsn'//'_htc'//zchar
+         z2d(:,:) = sxxsn(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'syysn'//'_htc'//zchar
+         z2d(:,:) = syysn(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxysn'//'_htc'//zchar
+         z2d(:,:) = sxysn(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxa'//'_htc'//zchar
+         z2d(:,:) = sxa(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sya'//'_htc'//zchar
+         z2d(:,:) = sya(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxxa'//'_htc'//zchar
+         z2d(:,:) = sxxa(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'syya'//'_htc'//zchar
+         z2d(:,:) = syya(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxya'//'_htc'//zchar
+         z2d(:,:) = sxya(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxc0'//'_htc'//zchar
+         z2d(:,:) = sxc0(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'syc0'//'_htc'//zchar
+         z2d(:,:) = syc0(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxxc0'//'_htc'//zchar
+         z2d(:,:) = sxxc0(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'syyc0'//'_htc'//zchar
+         z2d(:,:) = syyc0(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxyc0'//'_htc'//zchar
+         z2d(:,:) = sxyc0(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxsal'//'_htc'//zchar
+         z2d(:,:) = sxsal(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sysal'//'_htc'//zchar
+         z2d(:,:) = sysal(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxxsal'//'_htc'//zchar
+         z2d(:,:) = sxxsal(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'syysal'//'_htc'//zchar
+         z2d(:,:) = syysal(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxysal'//'_htc'//zchar
+         z2d(:,:) = sxysal(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxage'//'_htc'//zchar
+         z2d(:,:) = sxage(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'syage'//'_htc'//zchar
+         z2d(:,:) = syage(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxxage'//'_htc'//zchar
+         z2d(:,:) = sxxage(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'syyage'//'_htc'//zchar
+         z2d(:,:) = syyage(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         znam = 'sxyage'//'_htc'//zchar
+         z2d(:,:) = sxyage(:,:,jl)
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+      END DO
+
+      CALL iom_rstput( iter, nitrst, numriw, 'sxopw ' ,  sxopw  )
+      CALL iom_rstput( iter, nitrst, numriw, 'syopw ' ,  syopw  )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxxopw' ,  sxxopw )
+      CALL iom_rstput( iter, nitrst, numriw, 'syyopw' ,  syyopw )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxyopw' ,  sxyopw )
+
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I1)') jl
+         DO jk = 1, nlay_i 
+            WRITE(zchar1,'(I1)') jk
+            znam = 'sxe'//'_il'//zchar1//'_htc'//zchar
+            z2d(:,:) = sxe(:,:,jk,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sye'//'_il'//zchar1//'_htc'//zchar
+            z2d(:,:) = sye(:,:,jk,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar
+            z2d(:,:) = sxxe(:,:,jk,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'syye'//'_il'//zchar1//'_htc'//zchar
+            z2d(:,:) = syye(:,:,jk,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxye'//'_il'//zchar1//'_htc'//zchar
+            z2d(:,:) = sxye(:,:,jk,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         END DO
+      END DO
+
+      IF( iter == nitrst ) THEN
+         CALL iom_close( numriw )                         ! close the restart file
+         lrst_ice = .FALSE.
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, z2d )
+      !
+   END SUBROUTINE lim_rst_write
+
+
+   SUBROUTINE lim_rst_read
+      !!----------------------------------------------------------------------
+      !!                    ***  lim_rst_read  ***
+      !!
+      !! ** purpose  :   read of sea-ice variable restart in a netcdf file
+      !!----------------------------------------------------------------------
+      INTEGER :: ji, jj, jk, jl, indx
+      REAL(wp) ::   zfice, ziter
+      REAL(wp) ::   zs_inf, z_slope_s, zsmax, zsmin, zalpha, zindb   ! local scalars used for the salinity profile
+      REAL(wp), POINTER, DIMENSION(:)  ::   zs_zero 
+      REAL(wp), POINTER, DIMENSION(:,:) ::   z2d
+      CHARACTER(len=15) ::   znam
+      CHARACTER(len=1)  ::   zchar, zchar1
+      INTEGER           ::   jlibalt = jprstlib
+      LOGICAL           ::   llok
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( nlay_i, zs_zero )
+      CALL wrk_alloc( jpi, jpj, z2d )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'lim_rst_read : read ice NetCDF restart file'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+      ENDIF
+
+      IF ( jprstlib == jprstdimg ) THEN
+        ! eventually read netcdf file (monobloc)  for restarting on different number of processors
+        ! if {cn_icerst_in}.nc exists, then set jlibalt to jpnf90
+        INQUIRE( FILE = TRIM(cn_icerst_in)//'.nc', EXIST = llok )
+        IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF
+      ENDIF
+
+      CALL iom_open ( cn_icerst_in, numrir, kiolib = jprstlib )
+
+      CALL iom_get( numrir, 'nn_fsbc', zfice )
+      CALL iom_get( numrir, 'kt_ice' , ziter )    
+      IF(lwp) WRITE(numout,*) '   read ice restart file at time step    : ', ziter
+      IF(lwp) WRITE(numout,*) '   in any case we force it to nit000 - 1 : ', nit000 - 1
+
+      !Control of date
+
+      IF( ( nit000 - INT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 )   &
+         &     CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nit000 in ice restart',  &
+         &                   '   verify the file or rerun with the value 0 for the',        &
+         &                   '   control of time parameter  nrstdt' )
+      IF( INT(zfice) /= nn_fsbc          .AND. ABS( nrstdt ) == 1 )   &
+         &     CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nn_fsbc in ice restart',  &
+         &                   '   verify the file or rerun with the value 0 for the',         &
+         &                   '   control of time parameter  nrstdt' )
+
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I1)') jl
+         znam = 'v_i'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         v_i(:,:,jl) = z2d(:,:)
+         znam = 'v_s'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         v_s(:,:,jl) = z2d(:,:) 
+         znam = 'smv_i'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         smv_i(:,:,jl) = z2d(:,:)
+         znam = 'oa_i'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         oa_i(:,:,jl) = z2d(:,:)
+         znam = 'a_i'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         a_i(:,:,jl) = z2d(:,:)
+         znam = 't_su'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         t_su(:,:,jl) = z2d(:,:)
+         tn_ice (:,:,:) = t_su (:,:,:)
+      END DO
+
+      DO jl = 1, jpl 
+         CALL lbc_lnk( smv_i(:,:,jl) , 'T' ,  1. )
+         CALL lbc_lnk( v_i  (:,:,jl) , 'T' ,  1. )
+         CALL lbc_lnk( a_i  (:,:,jl) , 'T' ,  1. )
+      END DO
+
+      ! we first with bulk ice salinity
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zindb          = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - 1.0e-4 ) ) 
+               sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX(v_i(ji,jj,jl),1.0e-6) * zindb
+               ht_i(ji,jj,jl) = v_i(ji,jj,jl)   / MAX(a_i(ji,jj,jl),1.0e-6) * zindb
+            END DO
+         END DO
+      END DO
+
+      DO jk = 1, nlay_i
+         s_i(:,:,jk,:) = sm_i(:,:,:)
+      END DO
+
+      IF( num_sal == 2 ) THEN      ! Salinity profile
+         DO jl = 1, jpl
+            DO jk = 1, nlay_i
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     zs_inf        = sm_i(ji,jj,jl)
+                     z_slope_s     = 2._wp * sm_i(ji,jj,jl) / MAX( 0.01_wp , ht_i(ji,jj,jl) )
+                     !- slope of the salinity profile
+                     zs_zero(jk)   = z_slope_s * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) / REAL(nlay_i,wp)
+                     zsmax = 4.5_wp
+                     zsmin = 3.5_wp
+                     IF(     sm_i(ji,jj,jl) < zsmin ) THEN
+                        zalpha = 1._wp
+                     ELSEIF( sm_i(ji,jj,jl) < zsmax ) THEN
+                        zalpha = sm_i(ji,jj,jl) / ( zsmin - zsmax ) + zsmax / ( zsmax - zsmin )
+                     ELSE
+                        zalpha = 0._wp
+                     ENDIF
+                     s_i(ji,jj,jk,jl) = zalpha * zs_zero(jk) + ( 1._wp - zalpha ) * zs_inf
+                  END DO
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I1)') jl
+         znam = 'tempt_sl1'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         e_s(:,:,1,jl) = z2d(:,:)
+      END DO
+
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I1)') jl
+         DO jk = 1, nlay_i 
+            WRITE(zchar1,'(I1)') jk
+            znam = 'tempt'//'_il'//zchar1//'_htc'//zchar
+            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+            e_i(:,:,jk,jl) = z2d(:,:)
+         END DO
+      END DO
+
+      CALL iom_get( numrir, jpdom_autoglo, 'u_ice'     , u_ice      )
+      CALL iom_get( numrir, jpdom_autoglo, 'v_ice'     , v_ice      )
+      CALL iom_get( numrir, jpdom_autoglo, 'fsbbq'     , fsbbq      )
+      CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
+      CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
+      CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
+
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I1)') jl
+         znam = 'sxice'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxice(:,:,jl) = z2d(:,:)
+         znam = 'syice'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         syice(:,:,jl) = z2d(:,:)
+         znam = 'sxxice'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxxice(:,:,jl) = z2d(:,:)
+         znam = 'syyice'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         syyice(:,:,jl) = z2d(:,:)
+         znam = 'sxyice'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxyice(:,:,jl) = z2d(:,:)
+         znam = 'sxsn'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxsn(:,:,jl) = z2d(:,:)
+         znam = 'sysn'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sysn(:,:,jl) = z2d(:,:)
+         znam = 'sxxsn'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxxsn(:,:,jl) = z2d(:,:)
+         znam = 'syysn'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         syysn(:,:,jl) = z2d(:,:)
+         znam = 'sxysn'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxysn(:,:,jl) = z2d(:,:)
+         znam = 'sxa'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxa(:,:,jl) = z2d(:,:)
+         znam = 'sya'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sya(:,:,jl) = z2d(:,:)
+         znam = 'sxxa'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxxa(:,:,jl) = z2d(:,:)
+         znam = 'syya'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         syya(:,:,jl) = z2d(:,:)
+         znam = 'sxya'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxya(:,:,jl) = z2d(:,:)
+         znam = 'sxc0'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxc0(:,:,jl) = z2d(:,:)
+         znam = 'syc0'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         syc0(:,:,jl) = z2d(:,:)
+         znam = 'sxxc0'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxxc0(:,:,jl) = z2d(:,:)
+         znam = 'syyc0'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         syyc0(:,:,jl) = z2d(:,:)
+         znam = 'sxyc0'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxyc0(:,:,jl) = z2d(:,:)
+         znam = 'sxsal'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxsal(:,:,jl) = z2d(:,:)
+         znam = 'sysal'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sysal(:,:,jl) = z2d(:,:)
+         znam = 'sxxsal'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxxsal(:,:,jl) = z2d(:,:)
+         znam = 'syysal'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         syysal(:,:,jl) = z2d(:,:)
+         znam = 'sxysal'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxysal(:,:,jl) = z2d(:,:)
+         znam = 'sxage'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxage(:,:,jl) = z2d(:,:)
+         znam = 'syage'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         syage(:,:,jl) = z2d(:,:)
+         znam = 'sxxage'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxxage(:,:,jl) = z2d(:,:)
+         znam = 'syyage'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         syyage(:,:,jl) = z2d(:,:)
+         znam = 'sxyage'//'_htc'//zchar
+         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+         sxyage(:,:,jl)= z2d(:,:)
+      END DO
+
+      CALL iom_get( numrir, jpdom_autoglo, 'sxopw ' ,  sxopw  )
+      CALL iom_get( numrir, jpdom_autoglo, 'syopw ' ,  syopw  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxxopw' ,  sxxopw )
+      CALL iom_get( numrir, jpdom_autoglo, 'syyopw' ,  syyopw )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxyopw' ,  sxyopw )
+
+      DO jl = 1, jpl 
+         WRITE(zchar,'(I1)') jl
+         DO jk = 1, nlay_i 
+            WRITE(zchar1,'(I1)') jk
+            znam = 'sxe'//'_il'//zchar1//'_htc'//zchar
+            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+            sxe(:,:,jk,jl) = z2d(:,:)
+            znam = 'sye'//'_il'//zchar1//'_htc'//zchar
+            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+            sye(:,:,jk,jl) = z2d(:,:)
+            znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar
+            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+            sxxe(:,:,jk,jl) = z2d(:,:)
+            znam = 'syye'//'_il'//zchar1//'_htc'//zchar
+            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+            syye(:,:,jk,jl) = z2d(:,:)
+            znam = 'sxye'//'_il'//zchar1//'_htc'//zchar
+            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+            sxye(:,:,jk,jl) = z2d(:,:)
+         END DO
+      END DO
+      !
+      CALL iom_close( numrir )
+      !
+      CALL wrk_dealloc( nlay_i, zs_zero )
+      CALL wrk_dealloc( jpi, jpj, z2d )
+      !
+   END SUBROUTINE lim_rst_read
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :       Empty module            NO LIM sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE lim_rst_read             ! Empty routine
+   END SUBROUTINE lim_rst_read
+   SUBROUTINE lim_rst_write            ! Empty routine
+   END SUBROUTINE lim_rst_write
+#endif
+
+   !!======================================================================
+END MODULE limrst
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limsbc.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limsbc.F90	(revision 2113)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/limsbc.F90	(revision 2113)
@@ -0,0 +1,464 @@
+MODULE limsbc
+   !!======================================================================
+   !!                       ***  MODULE limsbc   ***
+   !!           computation of the flux at the sea ice/ocean interface
+   !!======================================================================
+   !! History :   -   ! 2006-07 (M. Vancoppelle)  LIM3 original code
+   !!            3.0  ! 2008-03 (C. Tallandier)  surface module
+   !!             -   ! 2008-04 (C. Tallandier)  split in 2 + new ice-ocean coupling
+   !!            3.3  ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea
+   !!                 !                  + simplification of the ice-ocean stress calculation
+   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
+   !!            3.5  ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                    LIM 3.0 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   lim_sbc_alloc : allocate the limsbc arrays
+   !!   lim_sbc_init  : initialisation
+   !!   lim_sbc_flx   : updates mass, heat and salt fluxes at the ocean surface
+   !!   lim_sbc_tau   : update i- and j-stresses, and its modulus at the ocean surface
+   !!----------------------------------------------------------------------
+   USE par_oce          ! ocean parameters
+   USE par_ice          ! ice parameters
+   USE dom_oce          ! ocean domain
+   USE sbc_ice          ! Surface boundary condition: sea-ice fields
+   USE sbc_oce          ! Surface boundary condition: ocean fields
+   USE phycst           ! physical constants
+   USE albedo           ! albedo parameters
+   USE ice              ! LIM sea-ice variables
+   USE lbclnk           ! ocean lateral boundary condition
+   USE in_out_manager   ! I/O manager
+   USE lib_mpp          ! MPP library
+   USE wrk_nemo         ! work arrays
+   USE prtctl           ! Print control
+   USE cpl_oasis3, ONLY : lk_cpl
+   USE oce,        ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n
+   USE dom_ice,    ONLY : tms
+   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   lim_sbc_init   ! called by ice_init
+   PUBLIC   lim_sbc_flx    ! called by sbc_ice_lim
+   PUBLIC   lim_sbc_tau    ! called by sbc_ice_lim
+
+   REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values
+   REAL(wp)  ::   rzero  = 0._wp    
+   REAL(wp)  ::   rone   = 1._wp
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2]
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean velocity   [m/s]
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   soce_0  , sice_0     ! cst SSS and ice salinity (levitating sea-ice) 
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)
+   !! $Id: limsbc.F90 3905 2013-05-24 13:40:39Z cetlod $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION lim_sbc_alloc()
+      !!-------------------------------------------------------------------
+      !!             ***  ROUTINE lim_sbc_alloc ***
+      !!-------------------------------------------------------------------
+      ALLOCATE( soce_0(jpi,jpj) , utau_oce(jpi,jpj) ,                       &
+         &      sice_0(jpi,jpj) , vtau_oce(jpi,jpj) , tmod_io(jpi,jpj), STAT=lim_sbc_alloc)
+         !
+      IF( lk_mpp             )   CALL mpp_sum( lim_sbc_alloc )
+      IF( lim_sbc_alloc /= 0 )   CALL ctl_warn('lim_sbc_alloc: failed to allocate arrays')
+   END FUNCTION lim_sbc_alloc
+
+
+   SUBROUTINE lim_sbc_flx( kt )
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE lim_sbc_flx ***
+      !!  
+      !! ** Purpose :   Update the surface ocean boundary condition for heat 
+      !!              salt and mass over areas where sea-ice is non-zero
+      !!         
+      !! ** Action  : - computes the heat and freshwater/salt fluxes
+      !!              at the ice-ocean interface.
+      !!              - Update the ocean sbc
+      !!     
+      !! ** Outputs : - qsr     : sea heat flux:     solar 
+      !!              - qns     : sea heat flux: non solar
+      !!              - emp     : freshwater budget: volume flux 
+      !!              - sfx     : salt flux 
+      !!              - fr_i    : ice fraction
+      !!              - tn_ice  : sea-ice surface temperature
+      !!              - alb_ice : sea-ice alberdo (lk_cpl=T)
+      !!
+      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
+      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt    ! number of iteration
+      !
+      INTEGER  ::   ji, jj, jl           ! dummy loop indices
+      INTEGER  ::   ierr, ifvt, i1mfr, idfr           ! local integer
+      INTEGER  ::   iflt, ial , iadv , ifral, ifrdv   !   -      -
+      REAL(wp) ::   zinda, zemp, zemp_snow, zfmm      ! local scalars
+      REAL(wp) ::   zemp_snw                          !   -      -
+      REAL(wp) ::   zfcm1 , zfcm2                     !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace
+      !!---------------------------------------------------------------------
+      
+      IF( lk_cpl )   CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp )
+
+      !------------------------------------------!
+      !      heat flux at the ocean surface      !
+      !------------------------------------------!
+      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
+      ! changed to old_frld and old ht_i
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zinda   = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
+            ifvt    = zinda  *  MAX( rzero , SIGN( rone, -phicif  (ji,jj) ) )  !subscripts are bad here
+            i1mfr   = 1.0 - MAX( rzero , SIGN( rone ,  - ( at_i(ji,jj)       ) ) )
+            idfr    = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(ji,jj) ) - pfrld(ji,jj) ) )
+            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )
+            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
+            iadv    = ( 1  - i1mfr ) * zinda
+            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
+            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
+
+            ! switch --- 1.0 ---------------- 0.0 --------------------
+            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            ! zinda   | if pfrld = 1       | if pfrld < 1            |
+            !  -> ifvt| if pfrld old_ht_i
+            ! i1mfr   | if frld = 1        | if frld  < 1            |
+            ! idfr    | if frld <= pfrld    | if frld > pfrld        |
+            ! iflt    | 
+            ! ial     |
+            ! iadv    |
+            ! ifral
+            ! ifrdv
+
+            !   computation the solar flux at ocean surface
+            IF (lk_cpl) THEN 
+               zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj)
+               ! MV 12 Oct 2012
+               zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj)
+               DO jl = 1, jpl
+                  zfcm1 = zfcm1 - qsr_ice(ji,jj,jl) * a_i(ji,jj,jl)
+               END DO
+               ! END MV 12 Oct 2012
+            ELSE
+               zfcm1   = pfrld(ji,jj) * qsr(ji,jj)  + ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj)
+            END IF
+            ! fstric     Solar flux transmitted trough the ice
+            ! qsr        Net short wave heat flux on free ocean
+            ! new line
+            fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj)
+
+            !  computation the non solar heat flux at ocean surface
+            zfcm2 = - zfcm1                                                                     & ! ???
+               &    + iflt    * fscmbq(ji,jj)                                                   & ! total ablation: heat given to the ocean
+               &    + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   &
+               &    + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice   &
+               &    + fhmec(ji,jj)                                                              & ! snow melt when ridging
+               &    + fheat_mec(ji,jj)                                                          & ! ridge formation
+               &    + fheat_res(ji,jj)                                                            ! residual heat flux
+            ! qcmif   Energy needed to bring the ocean surface layer until its freezing (ok)
+            ! qldif   heat balance of the lead (or of the open ocean)
+            ! qfvbq   latent heat uptake/release after accretion/ablation
+            ! qdtcn   Energy from the turbulent oceanic heat flux heat flux coming in the lead
+
+            IF( num_sal == 2 )   zfcm2 = zfcm2 + fhbri(ji,jj)    ! add contribution due to brine drainage 
+
+            ! bottom radiative component is sent to the computation of the oceanic heat flux
+            fsbbq(ji,jj) = ( 1._wp - ( ifvt + iflt ) ) * fscmbq(ji,jj)     
+
+            ! used to compute the oceanic heat flux at the next time step
+            qsr(ji,jj) = zfcm1                                       ! solar heat flux 
+            qns(ji,jj) = zfcm2 - fdtcn(ji,jj)                        ! non solar heat flux
+            !                           ! fdtcn : turbulent oceanic heat flux
+
+!!gm   this IF prevents the vertorisation of the whole loop
+            IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN
+               WRITE(numout,*) ' lim_sbc : heat fluxes '
+               WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx)
+               WRITE(numout,*) ' pfrld     : ', pfrld(jiindx,jjindx)
+               WRITE(numout,*) ' fstric    : ', fstric (jiindx,jjindx)
+               WRITE(numout,*)
+               WRITE(numout,*) ' qns       : ', qns(jiindx,jjindx)
+               WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx)
+               WRITE(numout,*) ' ifral     : ', ifral
+               WRITE(numout,*) ' ial       : ', ial  
+               WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx)
+               WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx)
+               WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice
+               WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice
+               WRITE(numout,*) ' ifrdv     : ', ifrdv
+               WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx)
+               WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx)
+               WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice
+               WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice
+               WRITE(numout,*) ' '
+               WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx)
+               WRITE(numout,*) ' fhmec     : ', fhmec(jiindx,jjindx)
+               WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx)
+               WRITE(numout,*) ' fhbri     : ', fhbri(jiindx,jjindx)
+               WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)
+            ENDIF
+!!gm   end
+         END DO
+      END DO
+
+      !------------------------------------------!
+      !      mass flux at the ocean surface      !
+      !------------------------------------------!
+
+!!gm   optimisation: this loop have to be merged with the previous one
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
+            !  ------------------------------------------------------------------------------------- 
+            !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
+            !  Thus  FW  flux  =  External ( E-P+snow melt)
+            !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
+            !                     Associated to Ice formation AND Ice melting
+            !                     Even if i see Ice melting as a FW and SALT flux
+            !        
+
+            !  computing freshwater exchanges at the ice/ocean interface
+            IF (lk_cpl) THEN 
+               zemp = - emp_tot(ji,jj) + emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !
+                  &   - rdm_snw(ji,jj) / rdt_ice
+            ELSE
+               zemp =   emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction
+                  &   - tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean
+                  &   + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice
+                  &   - fmmec(ji,jj)                                         ! snow falling when ridging
+            ENDIF
+
+            ! mass flux at the ocean/ice interface (sea ice fraction)
+            zemp_snw = rdm_snw(ji,jj) * r1_rdtice                         ! snow melting = pure water that enters the ocean
+            zfmm     = rdm_ice(ji,jj) * r1_rdtice                         ! Freezing minus mesting  
+
+            fmmflx(ji,jj) = zfmm                                     ! F/M mass flux save at least for biogeochemical model
+
+            emp(ji,jj) = zemp + zemp_snw + zfmm  ! mass flux + F/M mass flux (always ice/ocean mass exchange)
+            
+            !  correcting brine salt fluxes   (zinda = 1  if pfrld=1 , =0 otherwise)
+            zinda        = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
+            sfx_bri(ji,jj) = zinda * sfx_bri(ji,jj)
+         END DO
+      END DO
+
+      !------------------------------------------!
+      !      salt flux at the ocean surface      !
+      !------------------------------------------!
+
+      IF( num_sal == 2 ) THEN      ! variable ice salinity: brine drainage included in the salt flux
+         sfx(:,:) = sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) + sfx_bri(:,:)
+      ELSE                         ! constant ice salinity:
+         sfx(:,:) = sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:)
+      ENDIF
+      !-----------------------------------------------!
+      !   mass of snow and ice per unit area          !
+      !-----------------------------------------------!
+      IF( nn_ice_embd /= 0 ) THEN                               ! embedded sea-ice (mass required)
+         snwice_mass_b(:,:) = snwice_mass(:,:)                  ! save mass from the previous ice time step
+         !                                                      ! new mass per unit area
+         snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
+         !                                                      ! time evolution of snow+ice mass
+         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice
+      ENDIF
+
+      !-----------------------------------------------!
+      !   Storing the transmitted variables           !
+      !-----------------------------------------------!
+      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction            
+      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                      
+
+      !------------------------------------------------!
+      !    Computation of snow/ice and ocean albedo    !
+      !------------------------------------------------!
+      IF( lk_cpl ) THEN          ! coupled case
+         CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb )                  ! snow/ice albedo
+         !
+         alb_ice(:,:,:) =  0.5_wp * zalbp(:,:,:) + 0.5_wp * zalb (:,:,:)   ! Ice albedo (mean clear and overcast skys)
+      ENDIF
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl( tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns , clinfo2=' qns     : ' )
+         CALL prt_ctl( tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=sfx , clinfo2=' sfx     : ' )
+         CALL prt_ctl( tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ' )
+         CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl )
+      ENDIF
+      !
+      IF( lk_cpl )   CALL wrk_dealloc( jpi, jpj, jpl, zalb, zalbp )
+      ! 
+   END SUBROUTINE lim_sbc_flx
+
+
+   SUBROUTINE lim_sbc_tau( kt , pu_oce, pv_oce )
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE lim_sbc_tau ***
+      !!  
+      !! ** Purpose : Update the ocean surface stresses due to the ice
+      !!         
+      !! ** Action  : * at each ice time step (every nn_fsbc time step):
+      !!                - compute the modulus of ice-ocean relative velocity 
+      !!                  (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
+      !!                      tmod_io = rhoco * | U_ice-U_oce |
+      !!                - update the modulus of stress at ocean surface
+      !!                      taum = frld * taum + (1-frld) * tmod_io * | U_ice-U_oce |
+      !!              * at each ocean time step (every kt): 
+      !!                  compute linearized ice-ocean stresses as
+      !!                      Utau = tmod_io * | U_ice - pU_oce |
+      !!                using instantaneous current ocean velocity (usually before)
+      !!
+      !!    NB: - ice-ocean rotation angle no more allowed
+      !!        - here we make an approximation: taum is only computed every ice time step
+      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids 
+      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
+      !!
+      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
+      !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
+      !!---------------------------------------------------------------------
+      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
+      !!
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
+      REAL(wp) ::   zat_v, zvtau_ice, zv_t          !   -      -
+      !!---------------------------------------------------------------------
+      !
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
+!CDIR NOVERRCHK
+         DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
+!CDIR NOVERRCHK
+            DO ji = fs_2, fs_jpim1
+               !                                               ! 2*(U_ice-U_oce) at T-point
+               zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   
+               zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) 
+               !                                              ! |U_ice-U_oce|^2
+               zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )
+               !                                               ! update the ocean stress modulus
+               taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * rhoco * zmodt
+               tmod_io(ji,jj) = rhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point
+            END DO
+         END DO
+         CALL lbc_lnk( taum, 'T', 1. )   ;   CALL lbc_lnk( tmod_io, 'T', 1. )
+         !
+         utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
+         vtau_oce(:,:) = vtau(:,:)
+         !
+      ENDIF
+      !
+      !                                      !==  every ocean time-step  ==!
+      !
+      DO jj = 2, jpjm1                                !* update the stress WITHOUT a ice-ocean rotation angle
+         DO ji = fs_2, fs_jpim1   ! Vect. Opt.
+            zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp   ! ice area at u and V-points
+            zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
+            !                                                   ! linearized quadratic drag formulation
+            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
+            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
+            !                                                   ! stresses at the ocean surface
+            utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
+            vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
+         END DO
+      END DO
+      CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
+         &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask )
+      !  
+   END SUBROUTINE lim_sbc_tau
+
+
+   SUBROUTINE lim_sbc_init
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE lim_sbc_init  ***
+      !!             
+      !! ** Purpose : Preparation of the file ice_evolu for the output of
+      !!      the temporal evolution of key variables
+      !!
+      !! ** input   : Namelist namicedia
+      !!-------------------------------------------------------------------
+      !
+      INTEGER  ::   ji, jj                          ! dummy loop indices
+      REAL(wp) ::   zcoefu, zcoefv, zcoeff          ! local scalar
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'lim_sbc_init : LIM-3 sea-ice - surface boundary condition'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   '
+
+      !                                      ! allocate lim_sbc array
+      IF( lim_sbc_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' )
+      !
+      soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case
+      sice_0(:,:) = sice
+      !
+      IF( cp_cfg == "orca" ) THEN            ! decrease ocean & ice reference salinities in the Baltic sea 
+         WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &
+            &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
+            soce_0(:,:) = 4._wp
+            sice_0(:,:) = 2._wp
+         END WHERE
+      ENDIF
+      !                                      ! embedded sea ice
+      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
+         snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )
+         snwice_mass_b(:,:) = snwice_mass(:,:)
+      ELSE
+         snwice_mass  (:,:) = 0.0_wp         ! no mass exchanges
+         snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges
+      ENDIF
+      IF( nn_ice_embd == 2  .AND.         &  ! full embedment (case 2) & no restart
+         &  .NOT. ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area
+         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
+         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
+         !
+         ! Note: Changed the initial values of sshb and sshn=>  need to recompute ssh[u,v,f]_[b,n] 
+         !       which were previously set in domvvl
+         IF ( lk_vvl ) THEN            ! Is this necessary? embd 2 should be restricted to vvl only???
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
+                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
+                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
+                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
+                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
+                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
+                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
+                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
+                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
+                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
+                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
+                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
+               END DO
+            END DO
+            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
+            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1      ! NO Vector Opt.
+                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
+                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
+                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
+                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
+               END DO
+            END DO
+            CALL lbc_lnk( sshf_n, 'F', 1. )
+          ENDIF
+      ENDIF
+      !
+   END SUBROUTINE lim_sbc_init
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :        Dummy module       NO LIM 3.0 sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE lim_sbc           ! Dummy routine
+   END SUBROUTINE lim_sbc
+#endif 
+
+   !!======================================================================
+END MODULE limsbc
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbc_ice.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbc_ice.F90	(revision 2113)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbc_ice.F90	(revision 2113)
@@ -0,0 +1,151 @@
+MODULE sbc_ice
+   !!======================================================================
+   !!                 ***  MODULE  sbc_ice  ***
+   !! Surface module - LIM-3: parameters & variables defined in memory
+   !!======================================================================
+   !! History :  3.0  ! 2006-08  (G. Madec)  Surface module
+   !!            3.2  ! 2009-06  (S. Masson) merge with ice_oce
+   !!            3.3.1! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation
+   !!            3.4  ! 2011-11  (C. Harris) CICE added as an option
+   !!----------------------------------------------------------------------
+#if defined key_lim3 || defined key_lim2 || defined key_cice
+   !!----------------------------------------------------------------------
+   !!   'key_lim2' or 'key_lim3' :             LIM-2 or LIM-3 sea-ice model
+   !!----------------------------------------------------------------------
+   USE par_oce          ! ocean parameters
+# if defined key_lim3
+   USE par_ice          ! LIM-3 parameters
+# endif
+# if defined key_lim2
+   USE par_ice_2        ! LIM-2 parameters
+# endif
+# if defined key_cice 
+   USE ice_domain_size, only: ncat 
+#endif
+   USE lib_mpp          ! MPP library
+   USE in_out_manager   ! I/O manager
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC sbc_ice_alloc ! called in iceini(_2).F90
+
+   CHARACTER (len=8), PUBLIC :: cn_iceflx = 'none'                !: Flux handling over ice categories
+   LOGICAL, PUBLIC :: ln_iceflx_ave    = .FALSE. ! Average heat fluxes over all ice categories
+   LOGICAL, PUBLIC :: ln_iceflx_linear = .FALSE. ! Redistribute mean heat fluxes over all ice categories, using ice temperature and albedo
+
+# if defined  key_lim2
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim2    = .TRUE.   !: LIM-2 ice model
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3    = .FALSE.  !: no LIM-3
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE
+#  if defined key_lim2_vp
+   CHARACTER(len=1), PUBLIC, PARAMETER ::   cp_ice_msh = 'I'      !: VP : 'I'-grid ice-velocity (B-grid lower left corner)
+#  else
+   CHARACTER(len=1), PUBLIC, PARAMETER ::   cp_ice_msh = 'C'      !: EVP: 'C'-grid ice-velocity
+#  endif
+# endif
+# if defined  key_lim3
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim2    = .FALSE.  !: no LIM-2
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3    = .TRUE.   !: LIM-3 ice model
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE 
+   CHARACTER(len=1), PUBLIC, PARAMETER ::   cp_ice_msh = 'C'      !: 'C'-grid ice-velocity
+# endif
+# if defined  key_cice
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim2    = .FALSE.  !: no LIM-2
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3    = .FALSE.  !: no LIM-3
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .TRUE.   !: CICE ice model
+   CHARACTER(len=1), PUBLIC            ::   cp_ice_msh = 'F'      !: 'F'-grid ice-velocity
+# endif
+
+#if defined key_lim3 || defined key_lim2 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qns_ice   !: non solar heat flux over ice                  [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qsr_ice   !: solar heat flux over ice                      [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qla_ice   !: latent flux over ice                          [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dqla_ice  !: latent sensibility over ice                 [W/m2/K]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dqns_ice  !: non solar heat flux over ice (LW+SEN+LA)    [W/m2/K]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_ice    !: ice surface temperature                          [K]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   alb_ice   !: albedo of ice
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   utau_ice  !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vtau_ice  !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr1_i0    !: 1st Qsr fraction penetrating inside ice cover    [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr2_i0    !: 2nd Qsr fraction penetrating inside ice cover    [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice   !: sublimation-snow budget over ice             [kg/m2]
+
+
+# if defined key_lim3
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tatm_ice  !: air temperature
+# endif
+
+#elif defined key_cice
+   !
+   ! for consistency with LIM, these are declared with three dimensions
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qlw_ice            !: incoming long-wave
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qla_ice            !: latent flux over ice           [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qsr_ice            !: solar heat flux over ice       [W/m2]
+   !
+   ! other forcing arrays are two dimensional
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ss_iou             !: x ice-ocean surface stress at NEMO U point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ss_iov             !: y ice-ocean surface stress at NEMO V point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice            !: sublimation-snow budget over ice    [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tatm_ice           !: air temperature
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qatm_ice           !: specific humidity
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wndi_ice           !: i wind at T point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wndj_ice           !: j wind at T point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nfrzmlt            !: NEMO frzmlt
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iu              !: ice fraction at NEMO U point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iv              !: ice fraction at NEMO V point
+   !
+   ! finally, arrays corresponding to different ice categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i                !: category ice fraction
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   topmelt            !: category topmelt
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   botmelt            !: category botmelt
+#endif  
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id: sbc_ice.F90 3625 2012-11-21 13:19:18Z acc $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION sbc_ice_alloc()
+      !!----------------------------------------------------------------------
+      !!                     ***  FUNCTION sbc_ice_alloc  ***
+      !!----------------------------------------------------------------------
+#if defined key_lim3 || defined key_lim2
+      ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice (jpi,jpj,jpl) ,     &
+         &      qla_ice (jpi,jpj,jpl) , dqla_ice(jpi,jpj,jpl) ,     &
+         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl) ,     &
+         &      alb_ice (jpi,jpj,jpl) ,                             &
+         &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     ,     &
+         &      fr1_i0  (jpi,jpj)     , fr2_i0  (jpi,jpj)     ,     &
+#if defined key_lim3
+         &      emp_ice(jpi,jpj)      , tatm_ice(jpi,jpj)     , STAT= sbc_ice_alloc )
+#else
+         &      emp_ice(jpi,jpj)                              , STAT= sbc_ice_alloc )
+#endif
+#elif defined key_cice
+      ALLOCATE( qla_ice(jpi,jpj,1)    , qlw_ice(jpi,jpj,1)    , qsr_ice(jpi,jpj,1)    , &
+                wndi_ice(jpi,jpj)     , tatm_ice(jpi,jpj)     , qatm_ice(jpi,jpj)     , &
+                wndj_ice(jpi,jpj)     , nfrzmlt(jpi,jpj)      , ss_iou(jpi,jpj)       , &
+                ss_iov(jpi,jpj)       , fr_iu(jpi,jpj)        , fr_iv(jpi,jpj)        , &
+                a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat), STAT= sbc_ice_alloc )
+#endif
+         !
+      IF( lk_mpp            )   CALL mpp_sum ( sbc_ice_alloc )
+      IF( sbc_ice_alloc > 0 )   CALL ctl_warn('sbc_ice_alloc: allocation of arrays failed')
+   END FUNCTION sbc_ice_alloc
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option                      NO LIM 2.0 or 3.0 or CICE sea-ice model
+   !!----------------------------------------------------------------------
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim2    = .FALSE.  !: no LIM-2 ice model
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3    = .FALSE.  !: no LIM-3 ice model
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE  ice model
+   CHARACTER(len=1), PUBLIC, PARAMETER ::   cp_ice_msh = '-'      !: no grid ice-velocity
+#endif
+
+   !!======================================================================
+END MODULE sbc_ice
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbc_oce.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbc_oce.F90	(revision 2113)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbc_oce.F90	(revision 2113)
@@ -0,0 +1,169 @@
+MODULE sbc_oce
+   !!======================================================================
+   !!                       ***  MODULE  sbc_oce  ***
+   !! Surface module :   variables defined in core memory 
+   !!======================================================================
+   !! History :  3.0  ! 2006-06  (G. Madec)  Original code
+   !!             -   ! 2008-08  (G. Madec)  namsbc moved from sbcmod
+   !!            3.3  ! 2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
+   !!             -   ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
+   !!            3.3  ! 2010-10  (J. Chanut, C. Bricaud)  add the surface pressure forcing
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_oce_alloc : allocation of sbc arrays
+   !!   sbc_tau2wnd   : wind speed estimated from wind stress
+   !!----------------------------------------------------------------------
+   USE par_oce        ! ocean parameters
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc_oce_alloc   ! routine called in sbcmod.F90
+   PUBLIC   sbc_tau2wnd     ! routine called in several sbc modules
+   
+   !!----------------------------------------------------------------------
+   !!           Namelist for the Ocean Surface Boundary Condition
+   !!----------------------------------------------------------------------
+   !                                            !!* namsbc namelist *
+   LOGICAL , PUBLIC ::   ln_ana      = .FALSE.   !: analytical boundary condition flag
+   LOGICAL , PUBLIC ::   ln_flx      = .FALSE.   !: flux      formulation
+   LOGICAL , PUBLIC ::   ln_blk_clio = .FALSE.   !: CLIO bulk formulation
+   LOGICAL , PUBLIC ::   ln_blk_core = .FALSE.   !: CORE bulk formulation
+   LOGICAL , PUBLIC ::   ln_blk_mfs  = .FALSE.   !: MFS  bulk formulation
+   LOGICAL , PUBLIC ::   ln_cpl      = .FALSE.   !: coupled   formulation (overwritten by key_sbc_coupled )
+   LOGICAL , PUBLIC ::   ln_dm2dc    = .FALSE.   !: Daily mean to Diurnal Cycle short wave (qsr)
+   LOGICAL , PUBLIC ::   ln_rnf      = .FALSE.   !: runoffs / runoff mouths
+   LOGICAL , PUBLIC ::   ln_ssr      = .FALSE.   !: Sea Surface restoring on SST and/or SSS      
+   LOGICAL , PUBLIC ::   ln_apr_dyn  = .FALSE.   !: Atmospheric pressure forcing used on dynamics (ocean & ice)
+   LOGICAL , PUBLIC ::   ln_icebergs = .FALSE.   !: Icebergs
+   INTEGER , PUBLIC ::   nn_ice      = 0         !: flag for ice in the surface boundary condition (=0/1/2/3)
+   INTEGER , PUBLIC ::   nn_ice_embd = 0         !: flag for levitating/embedding sea-ice in the ocean
+   !                                             !: =0 levitating ice (no mass exchange, concentration/dilution effect)
+   !                                             !: =1 levitating ice with mass and salt exchange but no presure effect
+   !                                             !: =2 embedded sea-ice (full salt and mass exchanges and pressure)
+   INTEGER , PUBLIC ::   nn_fwb      = 0         !: FreshWater Budget: 
+   !                                             !:  = 0 unchecked 
+   !                                             !:  = 1 global mean of e-p-r set to zero at each nn_fsbc time step
+   !                                             !:  = 2 annual global mean of e-p-r set to zero
+   LOGICAL , PUBLIC ::   ln_wave     = .FALSE.   !: true if some coupling with wave model
+   LOGICAL , PUBLIC ::   ln_cdgw     = .FALSE.   !: true if neutral drag coefficient from wave model
+   LOGICAL , PUBLIC ::   ln_sdw      = .FALSE.   !: true if 3d stokes drift from wave model
+
+   !!----------------------------------------------------------------------
+   !!              Ocean Surface Boundary Condition fields
+   !!----------------------------------------------------------------------
+   LOGICAL , PUBLIC ::   lhftau = .FALSE.        !: HF tau used in TKE: mean(stress module) - module(mean stress)
+   !!                                   !!   now    ! before   !!
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum              !: module of sea surface stress (at T-point)    [N/m2] 
+   !! wndm is used onmpute surface gases exchanges in ice-free ocean or leads
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsr               !: sea heat flux:     solar                     [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qns    , qns_b    !: sea heat flux: non solar                     [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsr_tot           !: total     solar heat flux (over sea and ice) [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qns_tot           !: total non solar heat flux (over sea and ice) [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp    , emp_b    !: freshwater budget: volume flux               [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx    , sfx_b    !: salt flux                                    [PSU/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp_tot           !: total E-P over ocean and ice                 [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmflx            !: freshwater budget: freezing/melting          [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff   [Kg/m2/s]  
+   !!
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  sbc_tsc, sbc_tsc_b  !: sbc content trend                      [K.m/s] jpi,jpj,jpts
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  qsr_hc , qsr_hc_b   !: heat content trend due to qsr flux     [K.m/s] jpi,jpj,jpk
+   !!
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tprecip           !: total precipitation                          [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sprecip           !: solid precipitation                          [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fr_i              !: ice fraction = 1 - lead fraction      (between 0 to 1)
+#if defined key_cpl_carbon_cycle
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm]
+#endif
+
+   !!----------------------------------------------------------------------
+   !!                     Sea Surface Mean fields
+   !!----------------------------------------------------------------------
+   INTEGER , PUBLIC                     ::   nn_fsbc   !: frequency of sbc computation (as well as sea-ice model)
+   INTEGER , PUBLIC                     ::   nn_ice_dyn = 1   !: frequency of ice rheology/dynamics computation
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssu_m     !: mean (nn_fsbc time-step) surface sea i-current (U-point) [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssv_m     !: mean (nn_fsbc time-step) surface sea j-current (V-point) [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sst_m     !: mean (nn_fsbc time-step) surface sea temperature     [Celsius]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sss_m     !: mean (nn_fsbc time-step) surface sea salinity            [psu]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssh_m     !: mean (nn_fsbc time-step) sea surface height                [m]
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id: sbc_oce.F90 3905 2013-05-24 13:40:39Z cetlod $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION sbc_oce_alloc()
+      !!---------------------------------------------------------------------
+      !!                  ***  FUNCTION sbc_oce_alloc  ***
+      !!---------------------------------------------------------------------
+      INTEGER :: ierr(4)
+      !!---------------------------------------------------------------------
+      ierr(:) = 0
+      !
+      ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , taum(jpi,jpj) ,     &
+         &      vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , STAT=ierr(1) ) 
+         !
+      ALLOCATE( qns_tot(jpi,jpj) , qns  (jpi,jpj) , qns_b(jpi,jpj),        &
+         &      qsr_tot(jpi,jpj) , qsr  (jpi,jpj) ,                        &
+         &      emp    (jpi,jpj) , emp_b(jpi,jpj) ,                        &
+         &      sfx    (jpi,jpj) , sfx_b(jpi,jpj) , emp_tot(jpi,jpj), fmmflx(jpi,jpj), STAT=ierr(2) )
+         !
+      ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     &
+         &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) )
+         !
+      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     &
+#if defined key_cpl_carbon_cycle
+         &      atm_co2(jpi,jpj) ,                                        &
+#endif
+         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) ,                       &
+         &      ssv_m  (jpi,jpj) , sss_m  (jpi,jpj), ssh_m(jpi,jpj) , STAT=ierr(4) )
+         !
+      sbc_oce_alloc = MAXVAL( ierr )
+      IF( lk_mpp            )   CALL mpp_sum ( sbc_oce_alloc )
+      IF( sbc_oce_alloc > 0 )   CALL ctl_warn('sbc_oce_alloc: allocation of arrays failed')
+      !
+   END FUNCTION sbc_oce_alloc
+
+
+   SUBROUTINE sbc_tau2wnd
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc_tau2wnd  ***
+      !!                   
+      !! ** Purpose : Estimation of wind speed as a function of wind stress   
+      !!
+      !! ** Method  : |tau|=rhoa*Cd*|U|^2
+      !!---------------------------------------------------------------------
+      USE dom_oce         ! ocean space and time domain
+      USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
+      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
+      REAL(wp) ::   ztx, zty, ztau, zcoef ! temporary variables
+      INTEGER  ::   ji, jj                ! dummy indices
+      !!---------------------------------------------------------------------
+      zcoef = 0.5 / ( zrhoa * zcdrag ) 
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vect. opt.
+            ztx = utau(ji-1,jj  ) + utau(ji,jj) 
+            zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
+            ztau = SQRT( ztx * ztx + zty * zty )
+            wndm(ji,jj) = SQRT ( ztau * zcoef ) * tmask(ji,jj,1)
+         END DO
+      END DO
+      CALL lbc_lnk( wndm(:,:) , 'T', 1. )
+      !
+   END SUBROUTINE sbc_tau2wnd
+
+   !!======================================================================
+END MODULE sbc_oce
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcice_if.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcice_if.F90	(revision 2113)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcice_if.F90	(revision 2113)
@@ -0,0 +1,148 @@
+MODULE sbcice_if
+   !!======================================================================
+   !!                       ***  MODULE  sbcice  ***
+   !! Surface module :  update surface ocean boundary condition over ice
+   !!                   covered area using ice-if model
+   !!======================================================================
+   !! History :  3.0  !  2006-06  (G. Madec)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_ice_if    : update sbc in ice-covered area
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE eosbn2         ! equation of state
+   USE sbc_oce        ! surface boundary condition: ocean fields
+   USE sbccpl
+   USE fldread        ! read input field
+   USE iom            ! I/O manager library
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc_ice_if      ! routine called in sbcmod
+
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ice   ! structure of input ice-cover (file informations, fields read)
+   
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: sbcice_if.F90 3625 2012-11-21 13:19:18Z acc $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sbc_ice_if( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_ice_if  ***
+      !!
+      !! ** Purpose :   handle surface boundary condition over ice cover area
+      !!      when sea-ice model are not used
+      !!
+      !! ** Method  : - read sea-ice cover climatology
+      !!              - blah blah blah, ...
+      !!
+      !! ** Action  :   utau, vtau : remain unchanged
+      !!                taum, wndm : remain unchanged
+      !!                qns, qsr   : update heat flux below sea-ice
+      !!                emp, sfx   : update freshwater flux below sea-ice
+      !!                fr_i       : update the ice fraction
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      INTEGER  ::   ierror     ! return error code
+      REAL(wp) ::   ztrp, zsice, zt_fzp, zfr_obs
+      REAL(wp) ::   zqri, zqrj, zqrp, zqi
+      !!
+      CHARACTER(len=100) ::   cn_dir              ! Root directory for location of ice-if files
+      TYPE(FLD_N)        ::   sn_ice              ! informations about the fields to be read
+      NAMELIST/namsbc_iif/ cn_dir, sn_ice
+      !!---------------------------------------------------------------------
+      !                                         ! ====================== !
+      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
+         !                                      ! ====================== !
+         ! set file information
+         cn_dir = './'        ! directory in which the model is executed
+         ! ... default values (NB: frequency positive => hours, negative => months)
+         !             !   file    ! frequency !  variable  ! time intep !  clim  ! 'yearly' or ! weights  ! rotation   !
+         !             !   name    !  (hours)  !   name     !   (T/F)    !  (T/F) !  'monthly'  ! filename ! pairs      ! 
+         sn_ice = FLD_N('ice_cover',    -1    ,  'ice_cov' ,  .true.    , .true. ,   'yearly'  , ''       , ''         )
+
+         REWIND ( numnam )               ! ... read in namlist namiif
+         READ   ( numnam, namsbc_iif )
+
+         ALLOCATE( sf_ice(1), STAT=ierror )
+         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_if: unable to allocate sf_ice structure' )
+         ALLOCATE( sf_ice(1)%fnow(jpi,jpj,1) )
+         IF( sn_ice%ln_tint )   ALLOCATE( sf_ice(1)%fdta(jpi,jpj,1,2) )
+
+         ! fill sf_ice with sn_ice and control print
+         CALL fld_fill( sf_ice, (/ sn_ice /), cn_dir, 'sbc_ice_if', 'ice-if sea-ice model', 'namsbc_iif' )
+         !
+      ENDIF
+
+      CALL fld_read( kt, nn_fsbc, sf_ice )           ! Read input fields and provides the
+      !                                              ! input fields at the current time-step
+      
+      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
+         !
+         ztrp = -40.             ! restoring terme for temperature (w/m2/k)
+         zsice = - 0.04 / 0.8    ! ratio of isohaline compressibility over isotherme compressibility
+                                 ! ( d rho / dt ) / ( d rho / ds )      ( s = 34, t = -1.8 )
+         
+         fr_i(:,:) = tfreez( sss_m ) * tmask(:,:,1)      ! sea surface freezing temperature [Celcius]
+
+! OM : probleme. a_i pas defini dans les cas lim3 et cice
+#if defined key_coupled && defined key_lim2
+         a_i(:,:,1) = fr_i(:,:)         
+#endif
+
+         ! Flux and ice fraction computation
+!CDIR COLLAPSE
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               zt_fzp  = fr_i(ji,jj)                        ! freezing point temperature
+               zfr_obs = sf_ice(1)%fnow(ji,jj,1)            ! observed ice cover
+               !                                            ! ocean ice fraction (0/1) from the freezing point temperature
+               IF( sst_m(ji,jj) <= zt_fzp ) THEN   ;   fr_i(ji,jj) = 1.e0
+               ELSE                                ;   fr_i(ji,jj) = 0.e0
+               ENDIF
+
+               tsn(ji,jj,1,jp_tem) = MAX( tsn(ji,jj,1,jp_tem), zt_fzp )     ! avoid over-freezing point temperature
+
+               qsr(ji,jj) = ( 1. - zfr_obs ) * qsr(ji,jj)   ! solar heat flux : zero below observed ice cover
+
+               !                                            ! non solar heat flux : add a damping term 
+               !      # ztrp*(t-(tgel-1.))  if observed ice and no opa ice   (zfr_obs=1 fr_i=0)
+               !      # ztrp*min(0,t-tgel)  if observed ice and opa ice      (zfr_obs=1 fr_i=1)
+               zqri = ztrp * ( tsb(ji,jj,1,jp_tem) - ( zt_fzp - 1.) )
+               zqrj = ztrp * MIN( 0., tsb(ji,jj,1,jp_tem) - zt_fzp )
+               zqrp = ( zfr_obs * ( (1. - fr_i(ji,jj) ) * zqri    &
+                 &                 +      fr_i(ji,jj)   * zqrj ) ) * tmask(ji,jj,1)
+
+               !                                            ! non-solar heat flux 
+               !      # qns unchanged              if no climatological ice              (zfr_obs=0)
+               !      # qns = zqrp                 if climatological ice and no opa ice  (zfr_obs=1, fr_i=0)
+               !      # qns = zqrp -2(-4) watt/m2  if climatological ice and opa ice     (zfr_obs=1, fr_i=1)
+               !                                   (-2=arctic, -4=antarctic)   
+               zqi = -3. + SIGN( 1.e0, ff(ji,jj) )
+               qns(ji,jj) = ( ( 1.- zfr_obs ) * qns(ji,jj)                             &
+                  &          +      zfr_obs   * fr_i(ji,jj) * zqi ) * tmask(ji,jj,1)   &
+                  &       + zqrp
+            END DO
+         END DO
+         !
+      ENDIF
+      !
+   END SUBROUTINE sbc_ice_if
+
+   !!======================================================================
+END MODULE sbcice_if
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcice_lim.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcice_lim.F90	(revision 2113)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcice_lim.F90	(revision 2113)
@@ -0,0 +1,799 @@
+MODULE sbcice_lim
+   !!======================================================================
+   !!                       ***  MODULE  sbcice_lim  ***
+   !! Surface module :  update the ocean surface boundary condition over ice
+   !!       &           covered area using LIM sea-ice model
+   !! Sea-Ice model  :  LIM-3 Sea ice model time-stepping
+   !!=====================================================================
+   !! History :  2.0  ! 2006-12  (M. Vancoppenolle) Original code
+   !!            3.0  ! 2008-02  (C. Talandier)  Surface module from icestp.F90
+   !!             -   ! 2008-04  (G. Madec)  sltyle and lim_ctl routine
+   !!            3.3  ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
+   !!            3.4  ! 2011-01  (A Porter)  dynamical allocation
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3' :                                  LIM 3.0 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   sbc_ice_lim  : sea-ice model time-stepping and update ocean sbc over ice-covered area
+   !!   lim_ctl       : alerts in case of ice model crash
+   !!   lim_prt_state : ice control print at a given grid point
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE par_ice         ! sea-ice parameters
+   USE ice             ! LIM-3: ice variables
+   USE iceini          ! LIM-3: ice initialisation
+   USE dom_ice         ! LIM-3: ice domain
+
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE sbc_ice         ! Surface boundary condition: ice   fields
+   USE sbcblk_core     ! Surface boundary condition: CORE bulk
+   USE sbcblk_clio     ! Surface boundary condition: CLIO bulk
+   USE sbccpl          ! Surface boundary condition: coupled interface
+   USE albedo          ! ocean & ice albedo
+
+   USE phycst          ! Define parameters for the routines
+   USE eosbn2          ! equation of state
+   USE limdyn          ! Ice dynamics
+   USE limtrp          ! Ice transport
+   USE limthd          ! Ice thermodynamics
+   USE limitd_th       ! Thermodynamics on ice thickness distribution 
+   USE limitd_me       ! Mechanics on ice thickness distribution
+   USE limsbc          ! sea surface boundary condition
+   USE limdia          ! Ice diagnostics
+   USE limwri          ! Ice outputs
+   USE limrst          ! Ice restarts
+   USE limupdate       ! update of global variables
+   USE limvar          ! Ice variables switch
+
+   USE c1d             ! 1D vertical configuration
+   USE lbclnk          ! lateral boundary condition - MPP link
+   USE lib_mpp         ! MPP library
+   USE wrk_nemo        ! work arrays
+   USE iom             ! I/O manager library
+   USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   REAL (kind=wp), SAVE :: rdt_ice_dyn, rdt_ice_thermo
+
+   PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90
+   
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , UCL NEMO Consortium (2011)
+   !! $Id: sbcice_lim.F90 3625 2012-11-21 13:19:18Z acc $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   FUNCTION fice_cell_ave ( ptab)
+      !!--------------------------------------------------------------------------
+      !! * Compute average over categories, for grid cell (ice covered and free ocean)
+      !!--------------------------------------------------------------------------
+      REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave
+      REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab
+      INTEGER :: jl ! Dummy loop index
+      
+      fice_cell_ave (:,:) = 0.0_wp
+      
+      DO jl = 1, jpl
+         fice_cell_ave (:,:) = fice_cell_ave (:,:) &
+            &                  + a_i (:,:,jl) * ptab (:,:,jl)
+      END DO
+      
+   END FUNCTION fice_cell_ave
+   
+   FUNCTION fice_ice_ave ( ptab)
+      !!--------------------------------------------------------------------------
+      !! * Compute average over categories, for ice covered part of grid cell
+      !!--------------------------------------------------------------------------
+      REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave
+      REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab
+
+      fice_ice_ave (:,:) = 0.0_wp
+      WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)
+
+   END FUNCTION fice_ice_ave
+
+   !!======================================================================
+
+   SUBROUTINE sbc_ice_lim( kt, kblk )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_ice_lim  ***
+      !!                   
+      !! ** Purpose :   update the ocean surface boundary condition via the 
+      !!                Louvain la Neuve Sea Ice Model time stepping 
+      !!
+      !! ** Method  :   ice model time stepping
+      !!              - call the ice dynamics routine 
+      !!              - call the ice advection/diffusion routine 
+      !!              - call the ice thermodynamics routine 
+      !!              - call the routine that computes mass and 
+      !!                heat fluxes at the ice/ocean interface
+      !!              - save the outputs 
+      !!              - save the outputs for restart when necessary
+      !!
+      !! ** Action  : - time evolution of the LIM sea-ice model
+      !!              - update all sbc variables below sea-ice:
+      !!                utau, vtau, taum, wndm, qns , qsr, emp , sfx 
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt      ! ocean time step
+      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE)
+      !!
+      INTEGER  ::   jl      ! dummy loop index
+      REAL(wp) ::   zcoef   ! local scalar
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice      ! mean albedo of ice (for coupled)
+
+      REAL(wp), POINTER, DIMENSION(:,:) :: zalb_ice_all    ! Mean albedo over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: ztem_ice_all    ! Mean temperature over all categories
+      
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_ice_all   ! Mean solar heat flux over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_ice_all   ! Mean non solar heat flux over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_ice_all   ! Mean latent heat flux over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all  ! Mean d(qns)/dT over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqla_ice_all  ! Mean d(qla)/dT over all categories
+      !!----------------------------------------------------------------------
+
+      !- O.M. : why do we allocate all these arrays even when MOD( kt-1, nn_fsbc ) /= 0 ?????
+
+      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
+
+      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+         CALL wrk_alloc( jpi,jpj,jpl, zalb_ice)
+      END IF
+      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+         CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
+      ENDIF
+
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition' 
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping'
+         !
+         CALL ice_init
+         !
+         rdt_ice_dyn    = rdt_ice * REAL ( nn_ice_dyn, wp)
+         rdt_ice_thermo = rdt_ice
+         !
+         IF(lwp) WRITE(numout,*) 'nn_ice_dyn     : ', nn_ice_dyn
+         IF(lwp) WRITE(numout,*) 'rdt_ice        : ', rdt_ice
+         IF(lwp) WRITE(numout,*) 'rdt_ice_thermo : ', rdt_ice_thermo
+         IF(lwp) WRITE(numout,*) 'rdt_ice_dyn    : ', rdt_ice_dyn
+         !
+         IF( ln_nicep ) THEN      ! control print at a given point
+            jiindx = 44   ;   jjindx = 140
+            WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
+         ENDIF
+      ENDIF
+
+      !                                        !----------------------!
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  !
+         !                                     !----------------------!
+         !                                           !  Bulk Formulea !
+         !                                           !----------------!
+         !
+         u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point
+         v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean)
+         !
+         t_bo(:,:) = tfreez( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin]
+         !                                           ! (set to rt0 over land)
+         CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo
+
+         DO jl = 1, jpl
+            t_su(:,:,jl) = t_su(:,:,jl) +  rt0 * ( 1. - tmask(:,:,1) )
+         END DO
+
+         IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) )
+         
+         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+            !
+            ! Compute mean albedo and temperature
+            zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) ) 
+            ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) ) 
+            !
+         ENDIF
+                                                     ! Bulk formulea - provides the following fields:
+         ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2]
+         ! qsr_ice , qns_ice  : solar & non solar heat flux over ice   (T-point)         [W/m2]
+         ! qla_ice            : latent heat flux over ice              (T-point)         [W/m2]
+         ! dqns_ice, dqla_ice : non solar & latent heat sensistivity   (T-point)         [W/m2]
+         ! tprecip , sprecip  : total & solid precipitation            (T-point)         [Kg/m2/s]
+         ! fr1_i0  , fr2_i0   : 1sr & 2nd fraction of qsr penetration in ice             [%]
+         !
+         SELECT CASE( kblk )
+         CASE( 3 )                                       ! CLIO bulk formulation
+            CALL blk_ice_clio( t_su , zalb_ice_cs, zalb_ice_os,                           &
+               &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   &
+               &                      qla_ice    , dqns_ice   , dqla_ice  ,               &
+               &                      tprecip    , sprecip    ,                           &
+               &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  )
+            !         
+         CASE( 4 )                                       ! CORE bulk formulation
+            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice_cs,               &
+               &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   &
+               &                      qla_ice   , dqns_ice  , dqla_ice   ,               &
+               &                      tprecip   , sprecip   ,                            &
+               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  )
+            !
+         CASE ( 5 )
+            zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) )
+            
+            CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )
+
+            CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=tn_ice    )
+
+            ! Latent heat flux is forced to 0 in coupled :
+            !  it is included in qns (non-solar heat flux)
+            qla_ice  (:,:,:) = 0.0e0_wp
+            dqla_ice (:,:,:) = 0.0e0_wp
+            !
+         END SELECT
+
+         ! Average over all categories
+         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+
+            z_qns_ice_all  (:,:) = fice_ice_ave ( qns_ice  (:,:,:) )
+            z_qsr_ice_all  (:,:) = fice_ice_ave ( qsr_ice  (:,:,:) )
+            z_dqns_ice_all (:,:) = fice_ice_ave ( dqns_ice (:,:,:) )
+            z_qla_ice_all  (:,:) = fice_ice_ave ( qla_ice  (:,:,:) )
+            z_dqla_ice_all (:,:) = fice_ice_ave ( dqla_ice (:,:,:) )
+
+            DO jl = 1, jpl
+               dqns_ice (:,:,jl) = z_dqns_ice_all (:,:)
+               dqla_ice (:,:,jl) = z_dqla_ice_all (:,:)
+            END DO
+            !
+            IF ( ln_iceflx_ave ) THEN
+               DO jl = 1, jpl
+                  qns_ice  (:,:,jl) = z_qns_ice_all  (:,:)
+                  qsr_ice  (:,:,jl) = z_qsr_ice_all  (:,:)
+                  qla_ice  (:,:,jl) = z_qla_ice_all  (:,:)
+               END DO
+            END IF
+            !
+            IF ( ln_iceflx_linear ) THEN
+               DO jl = 1, jpl
+                  qns_ice  (:,:,jl) = z_qns_ice_all(:,:) + z_dqns_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:))
+                  qla_ice  (:,:,jl) = z_qla_ice_all(:,:) + z_dqla_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:))
+                  qsr_ice  (:,:,jl) = (1.0e0_wp-zalb_ice(:,:,jl)) / (1.0e0_wp-zalb_ice_all(:,:)) * z_qsr_ice_all(:,:)
+               END DO
+            END IF
+         END IF
+
+         !                                           !----------------------!
+         !                                           ! LIM-3  time-stepping !
+         !                                           !----------------------!
+         ! 
+         numit = numit + nn_fsbc                     ! Ice model time step
+         !
+         !                                           ! Store previous ice values
+!!gm : remark   old_...   should becomes ...b  as tn versus tb  
+         old_a_i  (:,:,:)   = a_i  (:,:,:)     ! ice area
+         old_e_i  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy
+         old_v_i  (:,:,:)   = v_i  (:,:,:)     ! ice volume
+         old_v_s  (:,:,:)   = v_s  (:,:,:)     ! snow volume 
+         old_e_s  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy
+         old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content
+         old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content
+
+         !                                           ! intialisation to zero    !!gm is it truly necessary ???
+         d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp
+         d_v_i_thd  (:,:,:)   = 0._wp   ;   d_v_i_trp  (:,:,:)   = 0._wp
+         d_e_i_thd  (:,:,:,:) = 0._wp   ;   d_e_i_trp  (:,:,:,:) = 0._wp
+         d_v_s_thd  (:,:,:)   = 0._wp   ;   d_v_s_trp  (:,:,:)   = 0._wp
+         d_e_s_thd  (:,:,:,:) = 0._wp   ;   d_e_s_trp  (:,:,:,:) = 0._wp
+         d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp
+         d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp
+         !
+         sfx    (:,:) = 0._wp
+         sfx_bri(:,:) = 0._wp   ;   sfx_mec  (:,:) = 0._wp   ;   sfx_res  (:,:) = 0._wp
+         fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp
+         fhmec  (:,:) = 0._wp   ;   
+         fmmec  (:,:) = 0._wp     
+         focea2D(:,:) = 0._wp
+         fsup2D (:,:) = 0._wp
+         ! 
+         diag_sni_gr(:,:) = 0._wp   ;   diag_lat_gr(:,:) = 0._wp
+         diag_bot_gr(:,:) = 0._wp   ;   diag_dyn_gr(:,:) = 0._wp
+         diag_bot_me(:,:) = 0._wp   ;   diag_sur_me(:,:) = 0._wp
+         ! dynamical invariants
+         delta_i(:,:) = 0._wp       ;   divu_i(:,:) = 0._wp       ;   shear_i(:,:) = 0._wp
+
+                          CALL lim_rst_opn( kt )     ! Open Ice restart file
+         !
+         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print
+         !
+         IF( .NOT. lk_c1d ) THEN                     ! Ice dynamics & transport (except in 1D case)
+            ! - Test OM, call ice dynamics every nn_ice_dyn ice time steps
+            IF( MOD( kt-1, nn_fsbc*nn_ice_dyn ) == 0 ) THEN
+                          rdt_ice = rdt_ice_dyn
+                          CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics )
+                          CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion )
+                          CALL lim_var_agg(1)             ! aggregate categories, requested
+                          CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting
+         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print
+                          CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting)
+            ENDIF
+         END IF
+                          rdt_ice = rdt_ice_thermo
+         !                                           ! Ice thermodynamics 
+                          CALL lim_var_glo2eqv            ! equivalent variables
+                          CALL lim_var_agg(1)             ! aggregate ice categories
+                          CALL lim_var_bv                 ! bulk brine volume (diag)
+                          CALL lim_thd( kt )              ! Ice thermodynamics 
+                          zcoef = rdt_ice /rday           !  Ice natural aging
+                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef
+                          CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin)
+         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print
+                          CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  !
+         !
+         !                                           ! Global variables update
+                          CALL lim_var_agg( 1 )           ! requested by limupdate
+                          CALL lim_update                 ! Global variables update
+                          CALL lim_var_glo2eqv            ! equivalent variables (outputs)
+                          CALL lim_var_agg(2)             ! aggregate ice thickness categories
+         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 2, ' - Final state - ' )   ! control print
+         !
+                          CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes
+         !
+         IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print
+         !
+         !                                           ! Diagnostics and outputs 
+         IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   &
+            &             CALL lim_dia 
+                          CALL lim_wri( 1  )              ! Ice outputs 
+         IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file 
+                          CALL lim_var_glo2eqv            ! ???
+         !
+         IF( ln_nicep )   CALL lim_ctl               ! alerts in case of model crash
+         !
+      ENDIF                                    ! End sea-ice time step only
+
+      !                                        !--------------------------!
+      !                                        !  at all ocean time step  !
+      !                                        !--------------------------!
+      !                                               
+      !                                              ! Update surface ocean stresses (only in ice-dynamic case)
+      !                                                   ! otherwise the atm.-ocean stresses are used everywhere
+      IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents
+      
+!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
+      !
+      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
+      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+         CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice)
+      END IF
+      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+         CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
+      ENDIF
+      !
+   END SUBROUTINE sbc_ice_lim
+
+
+   SUBROUTINE lim_ctl
+      !!-----------------------------------------------------------------------
+      !!                   ***  ROUTINE lim_ctl *** 
+      !!                 
+      !! ** Purpose :   Alerts in case of model crash
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk,  jl   ! dummy loop indices
+      INTEGER  ::   inb_altests       ! number of alert tests (max 20)
+      INTEGER  ::   ialert_id         ! number of the current alert
+      REAL(wp) ::   ztmelts           ! ice layer melting point
+      CHARACTER (len=30), DIMENSION(20)      ::   cl_alname   ! name of alert
+      INTEGER           , DIMENSION(20)      ::   inb_alp     ! number of alerts positive
+      !!-------------------------------------------------------------------
+
+      inb_altests = 10
+      inb_alp(:)  =  0
+
+      ! Alert if incompatible volume and concentration
+      ialert_id = 2 ! reference number of this alert
+      cl_alname(ialert_id) = ' Incompat vol and con         '    ! name of the alert
+
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN
+                  WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration '
+                  WRITE(numout,*) ' at_i     ', at_i(ji,jj)
+                  WRITE(numout,*) ' Point - category', ji, jj, jl
+                  WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl)
+                  WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl)
+                  WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
+                  WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
+                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+      ! Alerte if very thick ice
+      ialert_id = 3 ! reference number of this alert
+      cl_alname(ialert_id) = ' Very thick ice               ' ! name of the alert
+      jl = jpl 
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF(   ht_i(ji,jj,jl)  >  50._wp   ) THEN
+               CALL lim_prt_state( ji, jj, 2, ' ALERTE 3 :   Very thick ice ' )
+               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+            ENDIF
+         END DO
+      END DO
+
+      ! Alert if very fast ice
+      ialert_id = 4 ! reference number of this alert
+      cl_alname(ialert_id) = ' Very fast ice               ' ! name of the alert
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 0.5  .AND.  &
+               &  at_i(ji,jj) > 0._wp   ) THEN
+               CALL lim_prt_state( ji, jj, 1, ' ALERTE 4 :   Very fast ice ' )
+               WRITE(numout,*) ' ice strength             : ', strength(ji,jj)
+               WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj) 
+               WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj)
+               WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj) 
+               WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj)
+               WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj)
+               WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj)
+               WRITE(numout,*) ' sst                      : ', sst_m(ji,jj)
+               WRITE(numout,*) ' sss                      : ', sss_m(ji,jj)
+               WRITE(numout,*) 
+               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+            ENDIF
+         END DO
+      END DO
+
+      ! Alert if there is ice on continents
+      ialert_id = 6 ! reference number of this alert
+      cl_alname(ialert_id) = ' Ice on continents           ' ! name of the alert
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF(   tms(ji,jj) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN 
+               CALL lim_prt_state( ji, jj, 1, ' ALERTE 6 :   Ice on continents ' )
+               WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 
+               WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
+               WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
+               WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj)
+               WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj)
+               WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1)
+               WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj)
+               WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj)
+               !
+               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+            ENDIF
+         END DO
+      END DO
+
+!
+!     ! Alert if very fresh ice
+      ialert_id = 7 ! reference number of this alert
+      cl_alname(ialert_id) = ' Very fresh ice               ' ! name of the alert
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+!!gm  test twice sm_i ...  ????  bug?
+               IF( ( ( ABS( sm_i(ji,jj,jl) ) < 0.5 )   .OR.  &
+                     ( ABS( sm_i(ji,jj,jl) ) < 0.5 ) ) .AND. &
+                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN
+!                 CALL lim_prt_state(ji,jj,1, ' ALERTE 7 :   Very fresh ice ' )
+!                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
+!                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
+!                 WRITE(numout,*) ' s_i_newice           : ', s_i_newice(ji,jj,1:jpl)
+!                 WRITE(numout,*) 
+                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+               ENDIF
+            END DO
+         END DO
+      END DO
+!
+
+!     ! Alert if too old ice
+      ialert_id = 9 ! reference number of this alert
+      cl_alname(ialert_id) = ' Very old   ice               ' ! name of the alert
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( ( ( ABS( o_i(ji,jj,jl) ) > rdt_ice ) .OR. &
+                      ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. &
+                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN
+                  CALL lim_prt_state( ji, jj, 1, ' ALERTE 9 :   Wrong ice age ')
+                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+               ENDIF
+            END DO
+         END DO
+      END DO
+ 
+      ! Alert on salt flux
+      ialert_id = 5 ! reference number of this alert
+      cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN
+               CALL lim_prt_state( ji, jj, 3, ' ALERTE 5 :   High salt flux ' )
+               DO jl = 1, jpl
+                  WRITE(numout,*) ' Category no: ', jl
+                  WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)   
+                  WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl) 
+                  WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)   
+                  WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl) 
+                  WRITE(numout,*) ' '
+               END DO
+               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+            ENDIF
+         END DO
+      END DO
+
+      ! Alert if qns very big
+      ialert_id = 8 ! reference number of this alert
+      cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF(   ABS( qns(ji,jj) ) > 1500._wp  .AND.  at_i(ji,jj) > 0._wp  )  THEN
+               !
+               WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux'
+               WRITE(numout,*) ' ji, jj    : ', ji, jj
+               WRITE(numout,*) ' qns       : ', qns(ji,jj)
+               WRITE(numout,*) ' sst       : ', sst_m(ji,jj)
+               WRITE(numout,*) ' sss       : ', sss_m(ji,jj)
+               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj)
+               WRITE(numout,*) ' qldif     : ', qldif(ji,jj)
+               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice
+               WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice
+               WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj)
+               WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj)
+               WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice
+               WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice
+               WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj) 
+               WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj) 
+               WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj) 
+               WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 
+               WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj) 
+               !
+               CALL lim_prt_state( ji, jj, 2, '   ')
+               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+               !
+            ENDIF
+         END DO
+      END DO
+      !+++++
+ 
+      ! Alert if very warm ice
+      ialert_id = 10 ! reference number of this alert
+      cl_alname(ialert_id) = ' Very warm ice                ' ! name of the alert
+      inb_alp(ialert_id) = 0
+      DO jl = 1, jpl
+         DO jk = 1, nlay_i
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rtt
+                  IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-6   &
+                     &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN
+                     WRITE(numout,*) ' ALERTE 10 :   Very warm ice'
+                     WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl
+                     WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl)
+                     WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
+                     WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl)
+                     WRITE(numout,*) ' ztmelts : ', ztmelts
+                     inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+                  ENDIF
+               END DO
+            END DO
+         END DO
+      END DO
+
+      ialert_id = 1                                 ! reference number of this alert
+      cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert
+      WRITE(numout,*)
+      WRITE(numout,*) ' All alerts at the end of ice model '
+      DO ialert_id = 1, inb_altests
+         WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! '
+      END DO
+      !
+   END SUBROUTINE lim_ctl
+ 
+   
+   SUBROUTINE lim_prt_state( ki, kj, kn, cd1 )
+      !!-----------------------------------------------------------------------
+      !!                   ***  ROUTINE lim_prt_state *** 
+      !!                 
+      !! ** Purpose :   Writes global ice state on the (i,j) point 
+      !!                in ocean.ouput 
+      !!                3 possibilities exist 
+      !!                n = 1/-1 -> simple ice state (plus Mechanical Check if -1)
+      !!                n = 2    -> exhaustive state
+      !!                n = 3    -> ice/ocean salt fluxes
+      !!
+      !! ** input   :   point coordinates (i,j) 
+      !!                n : number of the option
+      !!-------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices
+      CHARACTER(len=*), INTENT(in) ::   cd1           !
+      !!
+      INTEGER :: jl
+      !!-------------------------------------------------------------------
+
+      WRITE(numout,*) cd1             ! print title
+
+      !----------------
+      !  Simple state
+      !----------------
+
+      IF ( kn == 1 .OR. kn == -1 ) THEN
+         WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
+         WRITE(numout,*) ' Simple state '
+         WRITE(numout,*) ' masks s,u,v   : ', tms(ki,kj), tmu(ki,kj), tmv(ki,kj)
+         WRITE(numout,*) ' lat - long    : ', gphit(ki,kj), glamt(ki,kj)
+         WRITE(numout,*) ' Time step     : ', numit
+         WRITE(numout,*) ' - Ice drift   '
+         WRITE(numout,*) '   ~~~~~~~~~~~ '
+         WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ki-1,kj)
+         WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ki,kj)
+         WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ki,kj-1)
+         WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ki,kj)
+         WRITE(numout,*) ' strength      : ', strength(ki,kj)
+         WRITE(numout,*)
+         WRITE(numout,*) ' - Cell values '
+         WRITE(numout,*) '   ~~~~~~~~~~~ '
+         WRITE(numout,*) ' cell area     : ', area(ki,kj)
+         WRITE(numout,*) ' at_i          : ', at_i(ki,kj)       
+         WRITE(numout,*) ' vt_i          : ', vt_i(ki,kj)       
+         WRITE(numout,*) ' vt_s          : ', vt_s(ki,kj)       
+         DO jl = 1, jpl
+            WRITE(numout,*) ' - Category (', jl,')'
+            WRITE(numout,*) ' a_i           : ', a_i(ki,kj,jl)
+            WRITE(numout,*) ' ht_i          : ', ht_i(ki,kj,jl)
+            WRITE(numout,*) ' ht_s          : ', ht_s(ki,kj,jl)
+            WRITE(numout,*) ' v_i           : ', v_i(ki,kj,jl)
+            WRITE(numout,*) ' v_s           : ', v_s(ki,kj,jl)
+            WRITE(numout,*) ' e_s           : ', e_s(ki,kj,1,jl)/1.0e9
+            WRITE(numout,*) ' e_i           : ', e_i(ki,kj,1:nlay_i,jl)/1.0e9
+            WRITE(numout,*) ' t_su          : ', t_su(ki,kj,jl)
+            WRITE(numout,*) ' t_snow        : ', t_s(ki,kj,1,jl)
+            WRITE(numout,*) ' t_i           : ', t_i(ki,kj,1:nlay_i,jl)
+            WRITE(numout,*) ' sm_i          : ', sm_i(ki,kj,jl)
+            WRITE(numout,*) ' smv_i         : ', smv_i(ki,kj,jl)
+            WRITE(numout,*)
+            WRITE(numout,*) ' Pathological case : ', patho_case(ki,kj,jl)
+         END DO
+      ENDIF
+      IF( kn == -1 ) THEN
+         WRITE(numout,*) ' Mechanical Check ************** '
+         WRITE(numout,*) ' Check what means ice divergence '
+         WRITE(numout,*) ' Total ice concentration ', at_i (ki,kj)
+         WRITE(numout,*) ' Total lead fraction     ', ato_i(ki,kj)
+         WRITE(numout,*) ' Sum of both             ', ato_i(ki,kj) + at_i(ki,kj)
+         WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ki,kj) + at_i(ki,kj) - 1.00
+      ENDIF
+
+
+      !--------------------
+      !  Exhaustive state
+      !--------------------
+
+      IF ( kn .EQ. 2 ) THEN
+         WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
+         WRITE(numout,*) ' Exhaustive state '
+         WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj)
+         WRITE(numout,*) ' Time step ', numit
+         WRITE(numout,*) 
+         WRITE(numout,*) ' - Cell values '
+         WRITE(numout,*) '   ~~~~~~~~~~~ '
+         WRITE(numout,*) ' cell area     : ', area(ki,kj)
+         WRITE(numout,*) ' at_i          : ', at_i(ki,kj)       
+         WRITE(numout,*) ' vt_i          : ', vt_i(ki,kj)       
+         WRITE(numout,*) ' vt_s          : ', vt_s(ki,kj)       
+         WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ki-1,kj)
+         WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ki,kj)
+         WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ki,kj-1)
+         WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ki,kj)
+         WRITE(numout,*) ' strength      : ', strength(ki,kj)
+         WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ki,kj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ki,kj)
+         WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ki,kj)  , ' old_v_ice     : ', old_v_ice(ki,kj)  
+         WRITE(numout,*)
+
+         DO jl = 1, jpl
+              WRITE(numout,*) ' - Category (',jl,')'
+              WRITE(numout,*) '   ~~~~~~~~         ' 
+              WRITE(numout,*) ' ht_i       : ', ht_i(ki,kj,jl)             , ' ht_s       : ', ht_s(ki,kj,jl)
+              WRITE(numout,*) ' t_i        : ', t_i(ki,kj,1:nlay_i,jl)
+              WRITE(numout,*) ' t_su       : ', t_su(ki,kj,jl)             , ' t_s        : ', t_s(ki,kj,1,jl)
+              WRITE(numout,*) ' sm_i       : ', sm_i(ki,kj,jl)             , ' o_i        : ', o_i(ki,kj,jl)
+              WRITE(numout,*) ' a_i        : ', a_i(ki,kj,jl)              , ' old_a_i    : ', old_a_i(ki,kj,jl)   
+              WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ki,kj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ki,kj,jl) 
+              WRITE(numout,*) ' v_i        : ', v_i(ki,kj,jl)              , ' old_v_i    : ', old_v_i(ki,kj,jl)   
+              WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ki,kj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ki,kj,jl) 
+              WRITE(numout,*) ' v_s        : ', v_s(ki,kj,jl)              , ' old_v_s    : ', old_v_s(ki,kj,jl)  
+              WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ki,kj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ki,kj,jl)
+              WRITE(numout,*) ' e_i1       : ', e_i(ki,kj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ki,kj,1,jl)/1.0e9 
+              WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ki,kj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ki,kj,1,jl)/1.0e9
+              WRITE(numout,*) ' e_i2       : ', e_i(ki,kj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ki,kj,2,jl)/1.0e9  
+              WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ki,kj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ki,kj,2,jl)/1.0e9
+              WRITE(numout,*) ' e_snow     : ', e_s(ki,kj,1,jl)            , ' old_e_snow : ', old_e_s(ki,kj,1,jl) 
+              WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ki,kj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ki,kj,1,jl)
+              WRITE(numout,*) ' smv_i      : ', smv_i(ki,kj,jl)            , ' old_smv_i  : ', old_smv_i(ki,kj,jl)   
+              WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ki,kj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ki,kj,jl) 
+              WRITE(numout,*) ' oa_i       : ', oa_i(ki,kj,jl)             , ' old_oa_i   : ', old_oa_i(ki,kj,jl)
+              WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ki,kj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ki,kj,jl)
+              WRITE(numout,*) ' Path. case : ', patho_case(ki,kj,jl)
+        END DO !jl
+
+        WRITE(numout,*)
+        WRITE(numout,*) ' - Heat / FW fluxes '
+        WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
+!       WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ki,kj)
+!       WRITE(numout,*) ' sfx        : ', sfx      (ki,kj)
+!       WRITE(numout,*) ' fsalt_res  : ', fsalt_res(ki,kj)
+        WRITE(numout,*) ' fmmec      : ', fmmec    (ki,kj)
+        WRITE(numout,*) ' fhmec      : ', fhmec    (ki,kj)
+        WRITE(numout,*) ' fhbri      : ', fhbri    (ki,kj)
+        WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ki,kj)
+        WRITE(numout,*) 
+        WRITE(numout,*) ' sst        : ', sst_m(ki,kj)  
+        WRITE(numout,*) ' sss        : ', sss_m(ki,kj)  
+        WRITE(numout,*) 
+        WRITE(numout,*) ' - Stresses '
+        WRITE(numout,*) '   ~~~~~~~~ '
+        WRITE(numout,*) ' utau_ice   : ', utau_ice(ki,kj) 
+        WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ki,kj)
+        WRITE(numout,*) ' utau       : ', utau    (ki,kj) 
+        WRITE(numout,*) ' vtau       : ', vtau    (ki,kj)
+        WRITE(numout,*) ' oc. vel. u : ', u_oce   (ki,kj)
+        WRITE(numout,*) ' oc. vel. v : ', v_oce   (ki,kj)
+     ENDIF
+
+     !---------------------
+     ! Salt / heat fluxes
+     !---------------------
+
+     IF ( kn .EQ. 3 ) THEN
+        WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj
+        WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
+        WRITE(numout,*) ' - Salt / Heat Fluxes '
+        WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
+        WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj)
+        WRITE(numout,*) ' Time step ', numit
+        WRITE(numout,*)
+        WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
+        WRITE(numout,*) ' qsr       : ', qsr(ki,kj)
+        WRITE(numout,*) ' qns       : ', qns(ki,kj)
+        WRITE(numout,*)
+        WRITE(numout,*) ' - Salt fluxes at bottom interface ***'
+        WRITE(numout,*) ' emp       : ', emp    (ki,kj)
+        WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ki,kj)
+        WRITE(numout,*) ' sfx       : ', sfx    (ki,kj)
+        WRITE(numout,*) ' sfx_res   : ', sfx_res(ki,kj)
+        WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ki,kj)
+        WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
+        WRITE(numout,*) ' fheat_res : ', fheat_res(ki,kj)
+        WRITE(numout,*)
+        WRITE(numout,*) ' - Momentum fluxes '
+        WRITE(numout,*) ' utau      : ', utau(ki,kj) 
+        WRITE(numout,*) ' vtau      : ', vtau(ki,kj)
+      ENDIF
+      WRITE(numout,*) ' '
+      !
+   END SUBROUTINE lim_prt_state
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE sbc_ice_lim ( kt, kblk )     ! Dummy routine
+      WRITE(*,*) 'sbc_ice_lim: You should not have seen this print! error?', kt, kblk
+   END SUBROUTINE sbc_ice_lim
+#endif
+
+   !!======================================================================
+END MODULE sbcice_lim
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcice_lim_2.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcice_lim_2.F90	(revision 2113)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcice_lim_2.F90	(revision 2113)
@@ -0,0 +1,277 @@
+MODULE sbcice_lim_2
+   !!======================================================================
+   !!                       ***  MODULE  sbcice_lim_2  ***
+   !! Surface module :  update surface ocean boundary condition over ice covered area using LIM sea-ice model
+   !! Sea-Ice model  :  LIM-2 Sea ice model time-stepping
+   !!======================================================================
+   !! History :  1.0   !  06-2006  (G. Madec)  from icestp_2.F90
+   !!            3.0   !  08-2008  (S. Masson, E. .... ) coupled interface
+   !!            3.3   !  05-2009  (G.Garric) addition of the lim2_evp case
+   !!----------------------------------------------------------------------
+#if defined key_lim2
+   !!----------------------------------------------------------------------
+   !!   'key_lim2' :                                    LIM-2 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   sbc_ice_lim_2   : sea-ice model time-stepping and update ocean sbc over ice-covered area
+   !!----------------------------------------------------------------------
+   USE oce              ! ocean dynamics and tracers
+   USE dom_oce          ! ocean space and time domain
+   USE ice_2
+   USE par_ice_2
+   USE iceini_2
+   USE dom_ice_2
+
+   USE sbc_oce          ! Surface boundary condition: ocean fields
+   USE sbc_ice          ! Surface boundary condition: ice   fields
+   USE sbcblk_core      ! Surface boundary condition: CORE bulk
+   USE sbcblk_clio      ! Surface boundary condition: CLIO bulk
+   USE sbccpl           ! Surface boundary condition: coupled interface
+   USE albedo
+
+   USE phycst           ! Define parameters for the routines
+   USE eosbn2           ! equation of state
+   USE limdyn_2
+   USE limtrp_2
+   USE limdmp_2
+   USE limthd_2
+   USE limsbc_2         ! sea surface boundary condition
+   USE limdia_2
+   USE limwri_2
+   USE limrst_2
+
+   USE c1d              ! 1D vertical configuration
+
+   USE lbclnk           ! lateral boundary condition - MPP link
+   USE lib_mpp          ! MPP library
+   USE wrk_nemo         ! work arrays
+   USE iom              ! I/O manager library
+   USE in_out_manager   ! I/O manager
+   USE prtctl           ! Print control
+
+# if defined key_agrif
+   USE agrif_ice
+   USE agrif_lim2_update
+# endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   REAL (kind=wp), SAVE :: rdt_ice_dyn, rdt_ice_thermo
+
+   PUBLIC sbc_ice_lim_2 ! routine called by sbcmod.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: sbcice_lim_2.F90 3680 2012-11-27 14:42:24Z rblod $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sbc_ice_lim_2( kt, ksbc )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_ice_lim_2  ***
+      !!                   
+      !! ** Purpose :   update the ocean surface boundary condition via the 
+      !!                Louvain la Neuve Sea Ice Model time stepping 
+      !!
+      !! ** Method  :   ice model time stepping
+      !!              - call the ice dynamics routine 
+      !!              - call the ice advection/diffusion routine 
+      !!              - call the ice thermodynamics routine 
+      !!              - call the routine that computes mass and 
+      !!                heat fluxes at the ice/ocean interface
+      !!              - save the outputs 
+      !!              - save the outputs for restart when necessary
+      !!
+      !! ** Action  : - time evolution of the LIM sea-ice model
+      !!              - update all sbc variables below sea-ice:
+      !!                utau, vtau, taum, wndm, qns , qsr, emp , sfx 
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt      ! ocean time step
+      INTEGER, INTENT(in) ::   ksbc    ! type of sbc ( =3 CLIO bulk ; =4 CORE bulk ; =5 coupled )
+      !!
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      REAL(wp), DIMENSION(:,:,:), POINTER :: zalb_ice_os   ! albedo of the ice under overcast sky
+      REAL(wp), DIMENSION(:,:,:), POINTER :: zalb_ice_cs   ! albedo of ice under clear sky
+      REAL(wp), DIMENSION(:,:,:), POINTER :: zsist         ! surface ice temperature (K)
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi,jpj,1, zalb_ice_os, zalb_ice_cs, zsist )
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc_ice_lim_2 : update ocean surface boudary condition' 
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM) time stepping'
+         !
+         CALL ice_init_2
+         !
+         rdt_ice_dyn    = rdt_ice * REAL ( nn_ice_dyn, wp)
+         rdt_ice_thermo = rdt_ice
+         !
+         IF(lwp) WRITE(numout,*) 'nn_ice_dyn     : ', nn_ice_dyn
+         IF(lwp) WRITE(numout,*) 'rdt_ice        : ', rdt_ice
+         IF(lwp) WRITE(numout,*) 'rdt_ice_thermo : ', rdt_ice_thermo
+         IF(lwp) WRITE(numout,*) 'rdt_ice_dyn    : ', rdt_ice_dyn
+         !
+# if defined key_agrif
+         IF( .NOT. Agrif_Root() ) CALL Agrif_InitValues_cont_lim2   ! AGRIF: set the meshes
+# endif
+      ENDIF
+
+      !                                        !----------------------!
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  !
+         !                                     !----------------------!
+
+         IF(lwp) WRITE(numout,*) 'sbc_ice_lim_2 : pas de temps glace : ', kt
+         
+# if defined key_agrif
+         IF( .NOT. Agrif_Root() ) lim_nbstep = MOD(lim_nbstep,Agrif_rhot()&
+         &*Agrif_PArent(nn_fsbc)/REAL(nn_fsbc)) + 1
+# endif
+         !  Bulk Formulea !
+         !----------------!
+         ! ... mean surface ocean current at ice dynamics point
+         SELECT CASE( cp_ice_msh )
+         CASE( 'I' )                  !== B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
+            DO jj = 2, jpj
+               DO ji = 2, jpi   ! NO vector opt. possible
+                  u_oce(ji,jj) = 0.5_wp * ( ssu_m(ji-1,jj  ) + ssu_m(ji-1,jj-1) ) * tmu(ji,jj)
+                  v_oce(ji,jj) = 0.5_wp * ( ssv_m(ji  ,jj-1) + ssv_m(ji-1,jj-1) ) * tmu(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( u_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices)
+            CALL lbc_lnk( v_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices)
+            !
+         CASE( 'C' )                  !== C-grid ice dynamics :   U & V-points (same as ocean)
+            u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point
+            v_oce(:,:) = ssv_m(:,:)
+            !
+         END SELECT
+
+         ! ... masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
+         tfu(:,:) = tfreez( sss_m ) +  rt0 
+
+         zsist (:,:,1) = sist (:,:) + rt0 * ( 1. - tmask(:,:,1) )
+
+         ! ... ice albedo (clear sky and overcast sky)
+         CALL albedo_ice( zsist, reshape( hicif, (/jpi,jpj,1/) ), &
+                                 reshape( hsnif, (/jpi,jpj,1/) ), &
+                          zalb_ice_cs, zalb_ice_os )
+
+         ! ... Sea-ice surface boundary conditions output from bulk formulae :
+         !     - utau_ice   ! surface ice stress i-component (I-point)   [N/m2]
+         !     - vtau_ice   ! surface ice stress j-component (I-point)   [N/m2]
+         !     - qns_ice    ! non solar heat flux over ice   (T-point)   [W/m2]
+         !     - qsr_ice    !     solar heat flux over ice   (T-point)   [W/m2]
+         !     - qla_ice    ! latent    heat flux over ice   (T-point)   [W/m2]
+         !     - dqns_ice   ! non solar heat sensistivity    (T-point)   [W/m2]
+         !     - dqla_ice   ! latent    heat sensistivity    (T-point)   [W/m2]
+         !     - tprecip    ! total precipitation            (T-point)   [Kg/m2/s]
+         !     - sprecip    ! solid precipitation            (T-point)   [Kg/m2/s]
+         !     - fr1_i0     ! 1sr fraction of qsr penetration in ice     [%]
+         !     - fr2_i0     ! 2nd fraction of qsr penetration in ice     [%]
+         !
+         SELECT CASE( ksbc )
+         CASE( 3 )           ! CLIO bulk formulation
+            CALL blk_ice_clio( zsist, zalb_ice_cs, zalb_ice_os,                         &
+               &                      utau_ice   , vtau_ice   , qns_ice    , qsr_ice,   &
+               &                      qla_ice    , dqns_ice   , dqla_ice   ,            &
+               &                      tprecip    , sprecip    ,                         &
+               &                      fr1_i0     , fr2_i0     , cp_ice_msh , jpl  )
+
+         CASE( 4 )           ! CORE bulk formulation
+            CALL blk_ice_core( zsist, u_ice      , v_ice      , zalb_ice_cs,            &
+               &                      utau_ice   , vtau_ice   , qns_ice    , qsr_ice,   &
+               &                      qla_ice    , dqns_ice   , dqla_ice   ,            &
+               &                      tprecip    , sprecip    ,                         &
+               &                      fr1_i0     , fr2_i0     , cp_ice_msh , jpl  )
+         CASE( 5 )           ! Coupled formulation : atmosphere-ice stress only (fluxes provided after ice dynamics)
+            CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )
+         END SELECT
+
+         CALL iom_put( 'utau_ice', utau_ice )     ! Wind stress over ice along i-axis at I-point
+         CALL iom_put( 'vtau_ice', vtau_ice )     ! Wind stress over ice along j-axis at I-point
+
+         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+            CALL prt_ctl_info( 'Ice Forcings ' )
+            CALL prt_ctl( tab2d_1=tprecip ,clinfo1=' sbc_ice_lim: precip  : ', tab2d_2=sprecip , clinfo2=' Snow    : ' )
+            CALL prt_ctl( tab2d_1=utau_ice,clinfo1=' sbc_ice_lim: utau_ice: ', tab2d_2=vtau_ice, clinfo2=' vtau_ice: ' )
+            CALL prt_ctl( tab2d_1=sst_m   ,clinfo1=' sbc_ice_lim: sst     : ', tab2d_2=sss_m   , clinfo2=' sss     : ' )
+            CALL prt_ctl( tab2d_1=u_oce   ,clinfo1=' sbc_ice_lim: u_io    : ', tab2d_2=v_oce   , clinfo2=' v_io    : ' )
+            CALL prt_ctl( tab2d_1=hsnif   ,clinfo1=' sbc_ice_lim: hsnif  1: ', tab2d_2=hicif   , clinfo2=' hicif   : ' )
+            CALL prt_ctl( tab2d_1=frld    ,clinfo1=' sbc_ice_lim: frld   1: ', tab2d_2=sist    , clinfo2=' sist    : ' )
+         ENDIF
+
+         ! ---------------- !
+         !  Ice model step  !
+         ! ---------------- !
+         numit = numit + nn_fsbc                           ! Ice model time step
+
+                           CALL lim_rst_opn_2  ( kt )  ! Open Ice restart file
+         IF( .NOT. lk_c1d ) THEN                       ! Ice dynamics & transport (except in 1D case)
+            ! - Test OM, call ice dynamics every nn_ice_dyn ice time steps
+            IF( MOD( kt-1, nn_fsbc*nn_ice_dyn ) == 0 ) THEN
+                           IF(lwp) WRITE(numout,*) 'sbc_ice_lim_2 : pas de temps glace dynamique : ', kt
+                           rdt_ice = rdt_ice_dyn
+                           CALL lim_dyn_2      ( kt )      ! Ice dynamics    ( rheology/dynamics )
+                           CALL lim_trp_2      ( kt )      ! Ice transport   ( Advection/diffusion )
+                           rdt_ice = rdt_ice_thermo            
+            END IF
+           rdt_ice = rdt_ice_thermo
+           IF( ln_limdmp ) CALL lim_dmp_2      ( kt )      ! Ice damping 
+         END IF
+#if defined key_coupled
+         !                                             ! Ice surface fluxes in coupled mode 
+         IF( ksbc == 5 )   THEN
+            a_i(:,:,1)=fr_i
+            CALL sbc_cpl_ice_flx( frld,                                              &
+            !                                optional arguments, used only in 'mixed oce-ice' case
+            &                                             palbi = zalb_ice_cs, psst = sst_m, pist = zsist )
+            sprecip(:,:) = - emp_ice(:,:)   ! Ugly patch, WARNING, in coupled mode, sublimation included in snow (parsub = 0.)
+         ENDIF
+#endif
+                           rdt_ice = rdt_ice_thermo
+                           CALL lim_thd_2      ( kt )      ! Ice thermodynamics 
+                           CALL lim_sbc_flx_2  ( kt )      ! update surface ocean mass, heat & salt fluxes 
+
+         IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   &
+            &              CALL lim_dia_2      ( kt )      ! Ice Diagnostics
+# if ! defined key_iomput
+                           CALL lim_wri_2      ( kt )      ! Ice outputs
+# endif
+         IF( lrst_ice  )   CALL lim_rst_write_2( kt )      ! Ice restart file
+         !
+# if defined key_agrif && defined key_lim2
+         IF( .NOT. Agrif_Root() )   CALL agrif_update_lim2( kt )
+# endif
+         !
+      ENDIF                                    ! End sea-ice time step only
+      !
+      !                                        !--------------------------!
+      !                                        !  at all ocean time step  !
+      !                                        !--------------------------!
+      !                                               
+      !                                              ! Update surface ocean stresses (only in ice-dynamic case)
+      !                                                   ! otherwise the atm.-ocean stresses are used everywhere
+      IF( ln_limdyn    )   CALL lim_sbc_tau_2( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents
+      !
+      CALL wrk_dealloc( jpi,jpj,1, zalb_ice_os, zalb_ice_cs, zsist )
+      !
+   END SUBROUTINE sbc_ice_lim_2
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option           Dummy module      NO LIM 2.0 sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE sbc_ice_lim_2 ( kt, ksbc )     ! Dummy routine
+      INTEGER, INTENT(in) ::   kt, ksbc    
+      WRITE(*,*) 'sbc_ice_lim_2: You should not have seen this print! error?', kt, ksbc
+   END SUBROUTINE sbc_ice_lim_2
+#endif
+
+   !!======================================================================
+END MODULE sbcice_lim_2
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcmod.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcmod.F90	(revision 2113)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcmod.F90	(revision 2113)
@@ -0,0 +1,438 @@
+MODULE sbcmod
+   !!======================================================================
+   !!                       ***  MODULE  sbcmod  ***
+   !! Surface module :  provide to the ocean its surface boundary condition
+   !!======================================================================
+   !! History :  3.0  ! 2006-07  (G. Madec)  Original code
+   !!            3.1  ! 2008-08  (S. Masson, A. Caubel, E. Maisonnave, G. Madec) coupled interface
+   !!            3.3  ! 2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
+   !!            3.3  ! 2010-10  (S. Masson)  add diurnal cycle
+   !!            3.3  ! 2010-09  (D. Storkey) add ice boundary conditions (BDY)
+   !!             -   ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
+   !!             -   ! 2010-10  (J. Chanut, C. Bricaud, G. Madec)  add the surface pressure forcing
+   !!            3.4  ! 2011-11  (C. Harris) CICE added as an option
+   !!            3.5  ! 2012-11  (A. Coward, G. Madec) Rethink of heat, mass and salt surface fluxes
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_init       : read namsbc namelist
+   !!   sbc            : surface ocean momentum, heat and freshwater boundary conditions
+   !!----------------------------------------------------------------------
+   USE oce              ! ocean dynamics and tracers
+   USE dom_oce          ! ocean space and time domain
+   USE phycst           ! physical constants
+   USE sbc_oce          ! Surface boundary condition: ocean fields
+   USE sbc_ice          ! Surface boundary condition: ice fields
+   USE sbcdcy           ! surface boundary condition: diurnal cycle
+   USE sbcssm           ! surface boundary condition: sea-surface mean variables
+   USE sbcapr           ! surface boundary condition: atmospheric pressure
+   USE sbcana           ! surface boundary condition: analytical formulation
+   USE sbcflx           ! surface boundary condition: flux formulation
+   USE sbcblk_clio      ! surface boundary condition: bulk formulation : CLIO
+   USE sbcblk_core      ! surface boundary condition: bulk formulation : CORE
+   USE sbcblk_mfs       ! surface boundary condition: bulk formulation : MFS
+   USE sbcice_if        ! surface boundary condition: ice-if sea-ice model
+   USE sbcice_lim       ! surface boundary condition: LIM 3.0 sea-ice model
+   USE sbcice_lim_2     ! surface boundary condition: LIM 2.0 sea-ice model
+   USE sbcice_cice      ! surface boundary condition: CICE    sea-ice model
+   USE sbccpl           ! surface boundary condition: coupled florulation
+   USE cpl_oasis3, ONLY:lk_cpl      ! are we in coupled mode?
+   USE sbcssr           ! surface boundary condition: sea surface restoring
+   USE sbcrnf           ! surface boundary condition: runoffs
+   USE sbcfwb           ! surface boundary condition: freshwater budget
+   USE closea           ! closed sea
+   USE bdy_par          ! for lk_bdy
+   USE bdyice_lim2      ! unstructured open boundary data  (bdy_ice_lim_2 routine)
+   USE icbstp           ! Icebergs!
+
+   USE prtctl           ! Print control                    (prt_ctl routine)
+   USE iom              ! IOM library
+   USE in_out_manager   ! I/O manager
+   USE lib_mpp          ! MPP library
+   USE timing           ! Timing
+   USE sbcwave          ! Wave module
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc        ! routine called by step.F90
+   PUBLIC   sbc_init   ! routine called by opa.F90
+   
+   INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations)
+      
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO-consortium (2011) 
+   !! $Id: sbcmod.F90 3905 2013-05-24 13:40:39Z cetlod $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sbc_init
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc_init ***
+      !!
+      !! ** Purpose :   Initialisation of the ocean surface boundary computation
+      !!
+      !! ** Method  :   Read the namsbc namelist and set derived parameters
+      !!                Call init routines for all other SBC modules that have one
+      !!
+      !! ** Action  : - read namsbc parameters
+      !!              - nsbc: type of sbc
+      !!----------------------------------------------------------------------
+      INTEGER ::   icpt   ! local integer
+      !!
+      NAMELIST/namsbc/ nn_fsbc, nn_ice_dyn, ln_ana   , ln_flx,  ln_blk_clio, ln_blk_core, ln_cpl,   &
+         &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   &
+         &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw, cn_iceflx
+      !!----------------------------------------------------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_init : surface boundary condition setting'
+         WRITE(numout,*) '~~~~~~~~ '
+      ENDIF
+
+      REWIND( numnam )           ! Read Namelist namsbc
+      READ  ( numnam, namsbc )
+
+      !                          ! overwrite namelist parameter using CPP key information
+      IF( Agrif_Root() ) THEN                ! AGRIF zoom
+        IF( lk_lim2 )   nn_ice      = 2
+        IF( lk_lim3 )   nn_ice      = 3
+        IF( lk_cice )   nn_ice      = 4
+      ENDIF
+      IF( cp_cfg == 'gyre' ) THEN            ! GYRE configuration
+          ln_ana      = .TRUE.   
+          nn_ice      =   0
+      ENDIF
+      
+      IF(lwp) THEN               ! Control print
+         WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)'
+         WRITE(numout,*) '           frequency update of sbc (and ice)             nn_fsbc     = ', nn_fsbc
+         WRITE(numout,*) '           frequency ice rheology and dynamics           nn_ice_dyn  = ', nn_ice_dyn
+         WRITE(numout,*) '           Type of sbc : '
+         WRITE(numout,*) '              analytical formulation                     ln_ana      = ', ln_ana
+         WRITE(numout,*) '              flux       formulation                     ln_flx      = ', ln_flx
+         WRITE(numout,*) '              CLIO bulk  formulation                     ln_blk_clio = ', ln_blk_clio
+         WRITE(numout,*) '              CORE bulk  formulation                     ln_blk_core = ', ln_blk_core
+         WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs
+         WRITE(numout,*) '              coupled    formulation (T if key_sbc_cpl)  ln_cpl      = ', ln_cpl
+         WRITE(numout,*) '              Flux handling over ice categories          cn_iceflx   = ', TRIM (cn_iceflx)
+
+         WRITE(numout,*) '           lk_lim2 : ', lk_lim2
+         WRITE(numout,*) '           lk_lim3 : ', lk_lim3
+
+         
+         WRITE(numout,*) '           Misc. options of sbc : '
+         WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn
+         WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice      = ', nn_ice 
+         WRITE(numout,*) '              ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd = ', nn_ice_embd
+         WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc 
+         WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf
+         WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr      = ', ln_ssr
+         WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb      = ', nn_fwb
+         WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea   = ', nn_closea
+      ENDIF
+
+      !   Flux handling over ice categories 
+      SELECT CASE ( TRIM (cn_iceflx))
+      CASE ('ave')
+         ln_iceflx_ave    = .TRUE.
+         ln_iceflx_linear = .FALSE.
+      CASE ('linear')
+         ln_iceflx_ave    = .FALSE.
+         ln_iceflx_linear = .TRUE.
+      CASE default
+         ln_iceflx_ave    = .FALSE.
+         ln_iceflx_linear = .FALSE.
+      END SELECT
+      WRITE(numout,*) '              Fluxes averaged over all ice categories         ln_iceflx_ave    = ', ln_iceflx_ave
+      WRITE(numout,*) '              Fluxes distributed linearly over ice categories ln_iceflx_linear = ', ln_iceflx_linear
+      !
+      !                              ! allocate sbc arrays
+      IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_oce arrays' )
+
+      !                          ! Checks:
+      IF( .NOT. ln_rnf ) THEN                      ! no specific treatment in vicinity of river mouths 
+         ln_rnf_mouth  = .false.                      
+         IF( sbc_rnf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_rnf arrays' )
+         nkrnf         = 0
+         rnf     (:,:) = 0.0_wp
+         rnfmsk  (:,:) = 0.0_wp
+         rnfmsk_z(:)   = 0.0_wp
+      ENDIF
+      IF( nn_ice == 0  )   fr_i(:,:) = 0.e0        ! no ice in the domain, ice fraction is always zero
+
+      sfx(:,:) = 0.0_wp                            ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero) 
+                                                   ! only if sea-ice is present
+ 
+      fmmflx(:,:) = 0.0_wp                        ! freezing-melting array initialisation
+
+      !                                            ! restartability   
+      IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   &
+          MOD( nstock             , nn_fsbc) /= 0 ) THEN 
+         WRITE(ctmp1,*) 'experiment length (', nitend - nit000 + 1, ') or nstock (', nstock,   &
+            &           ' is NOT a multiple of nn_fsbc (', nn_fsbc, ')'
+         CALL ctl_stop( ctmp1, 'Impossible to properly do model restart' )
+      ENDIF
+      !
+      IF( MOD( rday, REAL(nn_fsbc, wp) * rdt ) /= 0 )   &
+         &  CALL ctl_warn( 'nn_fsbc is NOT a multiple of the number of time steps in a day' )
+      !
+      IF( ( nn_ice == 2 .OR. nn_ice ==3 ) .AND. .NOT.( ln_blk_clio .OR. ln_blk_core .OR. lk_cpl ) )   &
+         &   CALL ctl_stop( 'LIM sea-ice model requires a bulk formulation or coupled configuration' )
+      IF( nn_ice == 4 .AND. .NOT.( ln_blk_core .OR. lk_cpl ) )   &
+         &   CALL ctl_stop( 'CICE sea-ice model requires ln_blk_core or lk_cpl' )
+      IF( nn_ice == 4 .AND. lk_agrif )   &
+         &   CALL ctl_stop( 'CICE sea-ice model not currently available with AGRIF' )
+      IF( ( nn_ice == 3 .OR. nn_ice == 4 ) .AND. nn_ice_embd == 0 )   &
+         &   CALL ctl_stop( 'LIM3 and CICE sea-ice models require nn_ice_embd = 1 or 2' )
+
+      IF( ln_iceflx_ave .AND. ln_iceflx_linear ) &
+         &   CALL ctl_stop( ' ln_iceflx_ave and ln_iceflx_linear options are not compatible' )
+
+      IF( ( nn_ice ==3 .AND. lk_cpl) .AND. .NOT. ( ln_iceflx_ave .OR. ln_iceflx_linear ) ) &
+         &   CALL ctl_stop( ' With lim3 coupled, either ln_iceflx_ave or ln_iceflx_linear must be set to .TRUE.' )
+      
+      IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag
+
+      IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) )   &
+         &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' )
+      
+      IF( ln_dm2dc .AND. ( ( NINT(rday) / ( nn_fsbc * NINT(rdt) ) )  < 8 ) )   &
+         &   CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
+
+      IF ( ln_wave ) THEN
+      !Activated wave module but neither drag nor stokes drift activated
+         IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN
+            CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )
+      !drag coefficient read from wave model definable only with mfs bulk formulae and core 
+         ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN       
+             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
+         ENDIF
+      ELSE
+      IF ( ln_cdgw .OR. ln_sdw  )                                         & 
+         &   CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but     &
+         & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ')
+      ENDIF 
+      
+      !                          ! Choice of the Surface Boudary Condition (set nsbc)
+      icpt = 0
+      IF( ln_ana          ) THEN   ;   nsbc =  1   ; icpt = icpt + 1   ;   ENDIF       ! analytical      formulation
+      IF( ln_flx          ) THEN   ;   nsbc =  2   ; icpt = icpt + 1   ;   ENDIF       ! flux            formulation
+      IF( ln_blk_clio     ) THEN   ;   nsbc =  3   ; icpt = icpt + 1   ;   ENDIF       ! CLIO bulk       formulation
+      IF( ln_blk_core     ) THEN   ;   nsbc =  4   ; icpt = icpt + 1   ;   ENDIF       ! CORE bulk       formulation
+      IF( ln_blk_mfs      ) THEN   ;   nsbc =  6   ; icpt = icpt + 1   ;   ENDIF       ! MFS  bulk       formulation
+      IF( ln_cpl          ) THEN   ;   nsbc =  5   ; icpt = icpt + 1   ;   ENDIF       ! Coupled         formulation
+      IF( cp_cfg == 'gyre') THEN   ;   nsbc =  0                       ;   ENDIF       ! GYRE analytical formulation
+      IF( lk_esopa        )            nsbc = -1                                       ! esopa test, ALL formulations
+      !
+      IF( icpt /= 1 .AND. .NOT.lk_esopa ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '           E R R O R in setting the sbc, one and only one namelist/CPP key option '
+         WRITE(numout,*) '                     must be choosen. You choose ', icpt, ' option(s)'
+         WRITE(numout,*) '                     We stop'
+         nstop = nstop + 1
+      ENDIF
+      IF(lwp) THEN
+         WRITE(numout,*)
+         IF( nsbc == -1 )   WRITE(numout,*) '              ESOPA test All surface boundary conditions'
+         IF( nsbc ==  0 )   WRITE(numout,*) '              GYRE analytical formulation'
+         IF( nsbc ==  1 )   WRITE(numout,*) '              analytical formulation'
+         IF( nsbc ==  2 )   WRITE(numout,*) '              flux formulation'
+         IF( nsbc ==  3 )   WRITE(numout,*) '              CLIO bulk formulation'
+         IF( nsbc ==  4 )   WRITE(numout,*) '              CORE bulk formulation'
+         IF( nsbc ==  5 )   WRITE(numout,*) '              coupled formulation'
+         IF( nsbc ==  6 )   WRITE(numout,*) '              MFS Bulk formulation'
+      ENDIF
+      !
+      IF( ln_ssr      )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation
+      !
+      IF( nn_ice == 4 )   CALL cice_sbc_init( nsbc )      ! CICE initialisation
+      !
+   END SUBROUTINE sbc_init
+
+
+   SUBROUTINE sbc( kt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc  ***
+      !!              
+      !! ** Purpose :   provide at each time-step the ocean surface boundary
+      !!                condition (momentum, heat and freshwater fluxes)
+      !!
+      !! ** Method  :   blah blah  to be written ????????? 
+      !!                CAUTION : never mask the surface stress field (tke sbc)
+      !!
+      !! ** Action  : - set the ocean surface boundary condition at before and now 
+      !!                time step, i.e.  
+      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, sfx_b, qrp_b, erp_b
+      !!                utau  , vtau  , qns  , qsr  , emp  , sfx  , qrp  , erp
+      !!              - updte the ice fraction : fr_i
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt       ! ocean time step
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('sbc')
+      !
+      !                                            ! ---------------------------------------- !
+      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
+         !                                         ! ---------------------------------------- !
+         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
+         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
+         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! qsr_b (:,:) = qsr (:,:)
+         emp_b(:,:) = emp(:,:)
+         sfx_b(:,:) = sfx(:,:)
+      ENDIF
+      !                                            ! ---------------------------------------- !
+      !                                            !        forcing field computation         !
+      !                                            ! ---------------------------------------- !
+      !
+      IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
+                                                         ! (caution called before sbc_ssm)
+      !
+      CALL sbc_ssm( kt )                                 ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
+      !                                                  ! averaged over nf_sbc time-step
+
+      IF (ln_wave) CALL sbc_wave( kt )
+                                                   !==  sbc formulation  ==!
+                                                            
+      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
+      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
+      CASE(  0 )   ;   CALL sbc_gyre    ( kt )                    ! analytical formulation : GYRE configuration
+      CASE(  1 )   ;   CALL sbc_ana     ( kt )                    ! analytical formulation : uniform sbc
+      CASE(  2 )   ;   CALL sbc_flx     ( kt )                    ! flux formulation
+      CASE(  3 )   ;   CALL sbc_blk_clio( kt )                    ! bulk formulation : CLIO for the ocean
+      CASE(  4 )   ;   CALL sbc_blk_core( kt )                    ! bulk formulation : CORE for the ocean
+      CASE(  5 )   ;   CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! coupled formulation
+      CASE(  6 )   ;   CALL sbc_blk_mfs ( kt )                    ! bulk formulation : MFS for the ocean
+      CASE( -1 )                                
+                       CALL sbc_ana     ( kt )                    ! ESOPA, test ALL the formulations
+                       CALL sbc_gyre    ( kt )                    !
+                       CALL sbc_flx     ( kt )                    !
+                       CALL sbc_blk_clio( kt )                    !
+                       CALL sbc_blk_core( kt )                    !
+                       CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   !
+      END SELECT
+
+      !                                            !==  Misc. Options  ==!
+      
+      SELECT CASE( nn_ice )                                       ! Update heat and freshwater fluxes over sea-ice areas
+      CASE(  1 )   ;         CALL sbc_ice_if   ( kt )                ! Ice-cover climatology ("Ice-if" model)
+      CASE(  2 )   ;         CALL sbc_ice_lim_2( kt, nsbc )          ! LIM-2 ice model
+              IF( lk_bdy )   CALL bdy_ice_lim_2( kt )                ! BDY boundary condition
+      CASE(  3 )   ;         CALL sbc_ice_lim  ( kt, nsbc )          ! LIM-3 ice model
+      CASE(  4 )   ;         CALL sbc_ice_cice ( kt, nsbc )          ! CICE ice model
+      END SELECT                                              
+
+      IF( ln_icebergs    )   CALL icb_stp( kt )                   ! compute icebergs
+
+      IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes
+ 
+      IF( ln_ssr         )   CALL sbc_ssr( kt )                   ! add SST/SSS damping term
+
+      IF( nn_fwb    /= 0 )   CALL sbc_fwb( kt, nn_fwb, nn_fsbc )  ! control the freshwater budget
+
+      IF( nn_closea == 1 )   CALL sbc_clo( kt )                   ! treatment of closed sea in the model domain 
+      !                                                           ! (update freshwater fluxes)
+!RBbug do not understand why see ticket 667
+      CALL lbc_lnk( emp, 'T', 1. )
+      !
+      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         !                                             ! ---------------------------------------- !
+         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
+            & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN 
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
+            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before i-stress  (U-point)
+            CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b )   ! before j-stress  (V-point)
+            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point)
+            ! The 3D heat content due to qsr forcing is treated in traqsr
+            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  ) ! before     solar heat flux (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'emp_b', emp_b  )    ! before     freshwater flux (T-point)
+            ! To ensure restart capability with 3.3x/3.4 restart files    !! to be removed in v3.6
+            IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN
+               CALL iom_get( numror, jpdom_autoglo, 'sfx_b', sfx_b )  ! before salt flux (T-point)
+            ELSE
+               sfx_b (:,:) = sfx(:,:)
+            ENDIF
+         ELSE                                                   !* no restart: set from nit000 values
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
+            utau_b(:,:) = utau(:,:) 
+            vtau_b(:,:) = vtau(:,:)
+            qns_b (:,:) = qns (:,:)
+            emp_b (:,:) = emp(:,:)
+            sfx_b (:,:) = sfx(:,:)
+         ENDIF
+      ENDIF
+      !                                                ! ---------------------------------------- !
+      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
+         !                                             ! ---------------------------------------- !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )
+         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau )
+         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  )
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
+         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
+         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b' , sfx )
+      ENDIF
+
+      !                                                ! ---------------------------------------- !
+      !                                                !        Outputs and control print         !
+      !                                                ! ---------------------------------------- !
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
+         CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux
+         CALL iom_put( "saltflx", sfx  )                        ! downward salt flux  
+                                                                ! (includes virtual salt flux beneath ice 
+                                                                ! in linear free surface case)
+         CALL iom_put( "fmmflx", fmmflx  )                      ! Freezing-melting water flux
+         CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux 
+         CALL iom_put( "qns"   , qns        )                   ! solar heat flux
+         CALL iom_put( "qsr"   ,       qsr  )                   ! solar heat flux
+         IF( nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction 
+      ENDIF
+      !
+      CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at 
+      CALL iom_put( "vtau", vtau )   ! j-wind stress    each time step in sea-ice)
+      CALL iom_put( "taum", taum )   ! wind stress module 
+      CALL iom_put( "wspd", wndm )   ! wind speed  module 
+      !
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(emp-rnf)        , clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(sfx-rnf)        , clinfo1=' sfx-rnf  - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab3d_1=tmask            , clinfo1=' tmask    - : ', mask1=tmask, ovlap=1, kdim=jpk )
+         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' sst      - : ', mask1=tmask, ovlap=1, kdim=1   )
+         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sss      - : ', mask1=tmask, ovlap=1, kdim=1   )
+         CALL prt_ctl(tab2d_1=utau             , clinfo1=' utau     - : ', mask1=umask,                      &
+            &         tab2d_2=vtau             , clinfo2=' vtau     - : ', mask2=vmask, ovlap=1 )
+      ENDIF
+
+      IF( kt == nitend )   CALL sbc_final         ! Close down surface module if necessary
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('sbc')
+      !
+   END SUBROUTINE sbc
+
+
+   SUBROUTINE sbc_final
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc_final  ***
+      !!
+      !! ** Purpose :   Finalize CICE (if used)
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_ice == 4 )   CALL cice_sbc_final
+      !
+   END SUBROUTINE sbc_final
+
+   !!======================================================================
+END MODULE sbcmod
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/stpctl.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/stpctl.F90	(revision 2113)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/stpctl.F90	(revision 2113)
@@ -0,0 +1,206 @@
+MODULE stpctl
+   !!======================================================================
+   !!                       ***  MODULE  stpctl  ***
+   !! Ocean run control :  gross check of the ocean time stepping
+   !!======================================================================
+   !! History :  OPA  ! 1991-03  (G. Madec) Original code
+   !!            6.0  ! 1992-06  (M. Imbard)
+   !!            8.0  ! 1997-06  (A.M. Treguier)
+   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
+   !!            2.0  ! 2009-07  (G. Madec)  Add statistic for time-spliting
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stp_ctl      : Control the run
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE sol_oce         ! ocean space and time domain variables 
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! distributed memory computing
+   USE dynspg_oce      ! pressure gradient schemes 
+   USE c1d             ! 1D vertical configuration
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC stp_ctl           ! routine called by step.F90
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: stpctl.F90 3294 2012-01-28 16:44:18Z rblod $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE stp_ctl( kt, kindic )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE stp_ctl  ***
+      !!                     
+      !! ** Purpose :   Control the run
+      !!
+      !! ** Method  : - Save the time step in numstp
+      !!              - Print it each 50 time steps
+      !!              - Print solver statistics in numsol 
+      !!              - Stop the run IF problem for the solver ( indec < 0 )
+      !!
+      !! ** Actions :   'time.step' file containing the last ocean time-step
+      !!                
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
+      INTEGER, INTENT( inout ) ::   kindic  ! indicator of solver convergence
+      !!
+      INTEGER  ::   ji, jj, jk              ! dummy loop indices
+      INTEGER  ::   ii, ij, ik              ! temporary integers
+      REAL(wp) ::   zumax, zsmin, zssh2     ! temporary scalars
+      INTEGER, DIMENSION(3) ::   ilocu      ! 
+      INTEGER, DIMENSION(2) ::   ilocs      ! 
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'stp_ctl : time-stepping control'
+         WRITE(numout,*) '~~~~~~~'
+         ! open time.step file
+         CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+      ENDIF
+
+      IF(lwp) WRITE ( numstp, '(1x, i8)' )   kt      !* save the current time step in numstp
+      IF(lwp) REWIND( numstp )                       !  --------------------------
+
+      !                                              !* Test maximum of velocity (zonal only)
+      !                                              !  ------------------------
+      !! zumax = MAXVAL( ABS( un(:,:,:) ) )                ! slower than the following loop on NEC SX5
+      zumax = 0.e0
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zumax = MAX(zumax,ABS(un(ji,jj,jk)))
+          END DO 
+        END DO 
+      END DO        
+      IF( lk_mpp )   CALL mpp_max( zumax )                 ! max over the global domain
+      !
+      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,' abs(U) max: ', zumax
+      !
+      IF( zumax > 20.e0 ) THEN
+         IF( lk_mpp ) THEN
+            CALL mpp_maxloc(ABS(un),umask,zumax,ii,ij,ik)
+         ELSE
+            ilocu = MAXLOC( ABS( un(:,:,:) ) )
+            ii = ilocu(1) + nimpp - 1
+            ij = ilocu(2) + njmpp - 1
+            ik = ilocu(3)
+         ENDIF
+         IF(lwp) THEN
+            WRITE(numout,cform_err)
+            WRITE(numout,*) ' stpctl: the zonal velocity is larger than 20 m/s'
+            WRITE(numout,*) ' ====== '
+            WRITE(numout,9400) kt, zumax, ii, ij, ik
+            WRITE(numout,*)
+            WRITE(numout,*) '          output of last fields in numwso'
+         ENDIF
+         kindic = -3
+      ENDIF
+9400  FORMAT (' kt=',i6,' max abs(U): ',1pg11.4,', i j k: ',3i5)
+
+      !                                              !* Test minimum of salinity
+      !                                              !  ------------------------
+      !! zsmin = MINVAL( tsn(:,:,1,jp_sal), mask = tmask(:,:,1) == 1.e0 )  slower than the following loop on NEC SX5
+      zsmin = 100.e0
+      DO jj = 2, jpjm1
+         DO ji = 1, jpi
+            IF( tmask(ji,jj,1) == 1) zsmin = MIN(zsmin,tsn(ji,jj,1,jp_sal))
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_min( zsmin )                ! min over the global domain
+      !
+      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,' SSS min:', zsmin
+      !
+      !!
+      ! Ajout OM.
+      IF( zsmin < 0.1) THEN 
+         IF (lk_mpp) THEN
+            CALL mpp_minloc ( tsn(:,:,1,jp_sal),tmask(:,:,1), zsmin, ii,ij )
+         ELSE
+            ilocs = MINLOC( tsn(:,:,1,jp_sal), mask = tmask(:,:,1) == 1.e0 )
+            ii = ilocs(1) + nimpp - 1
+            ij = ilocs(2) + njmpp - 1
+         ENDIF
+         !
+         IF(lwp) THEN
+            WRITE(numout,cform_err)
+            WRITE(numout,*) 'stp_ctl : sea surface salinity below 0.1'
+            WRITE(numout,*) '======= '
+            WRITE(numout,9500) kt, zsmin, ii, ij
+            WRITE(numout,*)
+         ENDIF
+         !!kindic = -3
+         !!
+         ! Correct minimum salinity
+         tsn(:,:,1,jp_sal) = MAX (  tsn(:,:,1,jp_sal), 0.1 )
+      ENDIF
+      !! Fin ajout O.M. 
+
+      !!
+      IF( zsmin < 0.) THEN 
+         IF (lk_mpp) THEN
+            CALL mpp_minloc ( tsn(:,:,1,jp_sal),tmask(:,:,1), zsmin, ii,ij )
+         ELSE
+            ilocs = MINLOC( tsn(:,:,1,jp_sal), mask = tmask(:,:,1) == 1.e0 )
+            ii = ilocs(1) + nimpp - 1
+            ij = ilocs(2) + njmpp - 1
+         ENDIF
+         !
+         IF(lwp) THEN
+            WRITE(numout,cform_err)
+            WRITE(numout,*) 'stp_ctl : NEGATIVE sea surface salinity'
+            WRITE(numout,*) '======= '
+            WRITE(numout,9500) kt, zsmin, ii, ij
+            WRITE(numout,*)
+            WRITE(numout,*) '          output of last fields in numwso'
+         ENDIF
+         kindic = -3
+      ENDIF
+9500  FORMAT (' kt=',i6,' min SSS: ',1pg11.4,', i j: ',2i5)
+
+      
+      IF( lk_c1d )  RETURN          ! No log file in case of 1D vertical configuration
+
+      ! log file (solver or ssh statistics)
+      ! --------
+      IF( lk_dynspg_flt ) THEN      ! elliptic solver statistics (if required)
+         !
+         IF(lwp) WRITE(numsol,9200) kt, niter, res, SQRT(epsr)/eps       ! Solver
+         !
+         IF( kindic < 0 .AND. zsmin > 0.e0 .AND. zumax <= 20.e0 ) THEN   ! create a abort file if problem found 
+            IF(lwp) THEN
+               WRITE(numout,*) ' stpctl: the elliptic solver DO not converge or explode'
+               WRITE(numout,*) ' ====== '
+               WRITE(numout,9200) kt, niter, res, sqrt(epsr)/eps
+               WRITE(numout,*)
+               WRITE(numout,*) ' stpctl: output of last fields'
+               WRITE(numout,*) ' ======  '
+            ENDIF
+         ENDIF
+         !
+      ELSE                                   !* ssh statistics (and others...)
+         IF( kt == nit000 .AND. lwp ) THEN   ! open ssh statistics file (put in solver.stat file)
+            CALL ctl_opn( numsol, 'solver.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+         ENDIF
+         !
+         zssh2 = SUM( sshn(:,:) * sshn(:,:) * tmask_i(:,:) )
+         IF( lk_mpp )   CALL mpp_sum( zssh2 )      ! sum over the global domain
+         !
+         IF(lwp) WRITE(numsol,9300) kt, zssh2, zumax, zsmin      ! ssh statistics
+         !
+      ENDIF
+
+9200  FORMAT('it:', i8, ' iter:', i4, ' r: ',e16.10, ' b: ',e16.10 )
+9300  FORMAT(' it :', i8, ' ssh2: ', e16.10, ' Umax: ',e16.10,' Smin: ',e16.10)
+      !
+   END SUBROUTINE stp_ctl
+
+   !!======================================================================
+END MODULE stpctl
