Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 6988)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 6989)
@@ -71,4 +71,7 @@
 &namicedyn     !   Ice dynamics and transport
 !------------------------------------------------------------------------------
+                  ! -- limtrp & limadv -- !
+   nn_limadv      =    0            !  choose the advection scheme (-1=Prather ; 0=Ultimate-Macho)
+   nn_limadv_ord  =    5            !  choose the order of the advection scheme (if nn_limadv=0)
                   ! -- limitd_me -- !
    nn_icestr      =    0            !  ice strength parameteriztaion                      
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6988)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6989)
@@ -200,4 +200,7 @@
 
    !                                     !!** ice-dynamics namelist (namicedyn) **
+                                          ! -- limtrp & limadv -- !
+   INTEGER , PUBLIC ::   nn_limadv        !: choose the advection scheme (-1=Prather ; 0=Ultimate-Macho)
+   INTEGER , PUBLIC ::   nn_limadv_ord    !: choose the order of the advection scheme (if Ultimate-Macho)   
                                           ! -- limitd_me -- !
    INTEGER , PUBLIC ::   nn_icestr        !: ice strength parameterization (0=Hibler79 1=Rothrock75)
@@ -281,6 +284,4 @@
    REAL(wp), PUBLIC ::   rdt_ice          !: ice time step
    REAL(wp), PUBLIC ::   r1_rdtice        !: = 1. / rdt_ice
-   REAL(wp), PUBLIC ::   usecc2           !:  = 1.0 / ( rn_ecc * rn_ecc )
-   REAL(wp), PUBLIC ::   rhoco            !: = rau0 * cio
    REAL(wp), PUBLIC ::   r1_nlay_i        !: 1 / nlay_i
    REAL(wp), PUBLIC ::   r1_nlay_s        !: 1 / nlay_s 
Index: anches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90	(revision 6988)
+++ 	(revision )
@@ -1,418 +1,0 @@
-MODULE limadv 
-   !!======================================================================
-   !!                       ***  MODULE limadv   ***
-   !! LIM sea-ice model : sea-ice advection
-   !!======================================================================
-   !! History :  LIM  ! 2008-03 (M. Vancoppenolle)  LIM-3 from LIM-2 code
-   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c.
-   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
-   !!--------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                     LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   lim_adv_x  : advection of sea ice on x axis
-   !!   lim_adv_y  : advection of sea ice on y axis
-   !!----------------------------------------------------------------------
-   USE dom_oce          ! ocean domain
-   USE ice              ! LIM-3 variables
-   USE lbclnk           ! lateral boundary condition - MPP exchanges
-   USE in_out_manager   ! I/O manager
-   USE prtctl           ! Print control
-   USE lib_mpp          ! MPP library
-   USE wrk_nemo         ! work arrays
-   USE lib_fortran      ! to use key_nosignedzero
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   lim_adv_x   ! called by lim_trp
-   PUBLIC   lim_adv_y   ! called by lim_trp
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE lim_adv_x( pdf, put , pcrh, psm , ps0 ,   &
-      &                  psx, psxx, psy , psyy, psxy )
-      !!---------------------------------------------------------------------
-      !!                **  routine lim_adv_x  **
-      !!  
-      !! ** purpose :   Computes and adds the advection trend to sea-ice
-      !!              variable on i-axis
-      !!
-      !! ** method  :   Uses Prather second order scheme that advects tracers
-      !!              but also theirquadratic forms. The method preserves
-      !!              tracer structures by conserving second order moments.
-      !!
-      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
-      !!--------------------------------------------------------------------
-      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
-      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
-      !! 
-      INTEGER  ::   ji, jj                               ! dummy loop indices
-      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars
-      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
-      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zalg, zalg1, zalg1q         !  -      -
-      !---------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
-      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
-
-      ! Limitation of moments.                                           
-
-      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zslpmax = MAX( 0._wp, ps0(ji,jj) )
-            zs1max  = 1.5 * zslpmax
-            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
-            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
-               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
-            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
-
-            ps0 (ji,jj) = zslpmax  
-            psx (ji,jj) = zs1new      * rswitch
-            psxx(ji,jj) = zs2new      * rswitch
-            psy (ji,jj) = psy (ji,jj) * rswitch
-            psyy(ji,jj) = psyy(ji,jj) * rswitch
-            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
-         END DO
-      END DO
-
-      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
-      psm (:,:)  = MAX( pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
-
-      !  Calculate fluxes and moments between boxes i<-->i+1              
-      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0 
-         DO ji = 1, jpi
-            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
-            zalf         =  MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
-            zalfq        =  zalf * zalf
-            zalf1        =  1.0 - zalf
-            zalf1q       =  zalf1 * zalf1
-            !
-            zfm (ji,jj)  =  zalf  *   psm (ji,jj)
-            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) )  )
-            zfx (ji,jj)  =  zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
-            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj) * zalfq
-            zfy (ji,jj)  =  zalf  * ( psy (ji,jj) + zalf1 * psxy(ji,jj) )
-            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj)
-            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj)
-
-            !  Readjust moments remaining in the box.
-            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
-            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
-            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
-            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
-            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
-            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
-            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
-         END DO
-      END DO
-
-      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
-         DO ji = 1, fs_jpim1
-            zalf          = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
-            zalg  (ji,jj) = zalf
-            zalfq         = zalf * zalf
-            zalf1         = 1.0 - zalf
-            zalg1 (ji,jj) = zalf1
-            zalf1q        = zalf1 * zalf1
-            zalg1q(ji,jj) = zalf1q
-            !
-            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji+1,jj)
-            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
-            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
-            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji+1,jj) * zalfq
-            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) )
-            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji+1,jj)
-            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji+1,jj)
-         END DO
-      END DO
-
-      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box. 
-         DO ji = fs_2, fs_jpim1
-            zbt  =       zbet(ji-1,jj)
-            zbt1 = 1.0 - zbet(ji-1,jj)
-            !
-            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
-            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
-            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
-            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
-            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
-            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
-            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
-         END DO
-      END DO
-
-      !   Put the temporary moments into appropriate neighboring boxes.    
-      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
-         DO ji = fs_2, fs_jpim1
-            zbt  =       zbet(ji-1,jj)
-            zbt1 = 1.0 - zbet(ji-1,jj)
-            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
-            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
-            zalf1       = 1.0 - zalf
-            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
-            !
-            ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj)
-            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
-            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               &
-               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
-               &                                                + zbt1 * psxx(ji,jj)
-            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
-               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
-               &                                                + zbt1 * psxy(ji,jj)
-            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
-            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
-         END DO
-      END DO
-
-      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
-         DO ji = fs_2, fs_jpim1
-            zbt  =       zbet(ji,jj)
-            zbt1 = 1.0 - zbet(ji,jj)
-            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
-            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
-            zalf1       = 1.0 - zalf
-            ztemp       = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
-            !
-            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
-            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
-            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  &
-               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )      &
-               &                                      + ( zalf1 - zalf ) * ztemp ) )
-            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  &
-               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
-            psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
-            psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
-         END DO
-      END DO
-
-      !-- Lateral boundary conditions
-      CALL lbc_lnk_multi( psm , 'T',  1., ps0 , 'T',  1.   &
-         &              , psx , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
-         &              , psxx, 'T',  1., psyy, 'T',  1.   &
-         &              , psxy, 'T',  1. )
-
-      IF(ln_ctl) THEN
-         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
-         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
-         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
-         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
-      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
-      !
-   END SUBROUTINE lim_adv_x
-
-
-   SUBROUTINE lim_adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
-      &                  psx, psxx, psy , psyy, psxy )
-      !!---------------------------------------------------------------------
-      !!                **  routine lim_adv_y  **
-      !!            
-      !! ** purpose :   Computes and adds the advection trend to sea-ice 
-      !!              variable on y axis
-      !!
-      !! ** method  :   Uses Prather second order scheme that advects tracers
-      !!              but also their quadratic forms. The method preserves 
-      !!              tracer structures by conserving second order moments.
-      !! 
-      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
-      !!---------------------------------------------------------------------
-      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
-      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
-      !!
-      INTEGER  ::   ji, jj                               ! dummy loop indices
-      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars
-      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
-      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zf0, zfx , zfy , zbet   ! 2D workspace
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zfm, zfxx, zfyy, zfxy   !  -      -
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zalg, zalg1, zalg1q     !  -      -
-      !---------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
-      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
-
-      ! Limitation of moments.
-
-      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zslpmax = MAX( 0._wp, ps0(ji,jj) )
-            zs1max  = 1.5 * zslpmax
-            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
-            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
-               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
-            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
-            !
-            ps0 (ji,jj) = zslpmax  
-            psx (ji,jj) = psx (ji,jj) * rswitch
-            psxx(ji,jj) = psxx(ji,jj) * rswitch
-            psy (ji,jj) = zs1new * rswitch
-            psyy(ji,jj) = zs2new * rswitch
-            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
-         END DO
-      END DO
-
-      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
-      psm(:,:)  = MAX(  pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
-
-      !  Calculate fluxes and moments between boxes j<-->j+1              
-      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
-         DO ji = 1, jpi
-            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
-            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
-            zalfq        =  zalf * zalf
-            zalf1        =  1.0 - zalf
-            zalf1q       =  zalf1 * zalf1
-            !
-            zfm (ji,jj)  =  zalf  * psm(ji,jj)
-            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
-            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
-            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
-            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
-            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
-            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
-            !
-            !  Readjust moments remaining in the box.
-            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
-            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
-            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
-            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
-            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
-            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
-            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
-         END DO
-      END DO
-      !
-      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
-         DO ji = 1, jpi
-            zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
-            zalg  (ji,jj) = zalf
-            zalfq         = zalf * zalf
-            zalf1         = 1.0 - zalf
-            zalg1 (ji,jj) = zalf1
-            zalf1q        = zalf1 * zalf1
-            zalg1q(ji,jj) = zalf1q
-            !
-            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1)
-            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
-            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
-            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq
-            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) )
-            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1)
-            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1)
-         END DO
-      END DO
-
-      !  Readjust moments remaining in the box. 
-      DO jj = 2, jpj
-         DO ji = 1, jpi
-            zbt  =         zbet(ji,jj-1)
-            zbt1 = ( 1.0 - zbet(ji,jj-1) )
-            !
-            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
-            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
-            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
-            psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj)
-            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
-            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
-            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
-         END DO
-      END DO
-
-      !   Put the temporary moments into appropriate neighboring boxes.    
-      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
-         DO ji = 1, jpi
-            zbt  =         zbet(ji,jj-1)
-            zbt1 = ( 1.0 - zbet(ji,jj-1) )
-            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
-            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
-            zalf1       = 1.0 - zalf
-            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
-            !
-            ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj)
-            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
-               &                                               + zbt1 * psy(ji,jj)  
-            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
-               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
-               &                                               + zbt1 * psyy(ji,jj)
-            psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               &
-               &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   &
-               &                                                + zbt1 * psxy(ji,jj)
-            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
-            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
-         END DO
-      END DO
-
-      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
-         DO ji = 1, jpi
-            zbt  =         zbet(ji,jj)
-            zbt1 = ( 1.0 - zbet(ji,jj) )
-            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
-            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
-            zalf1       = 1.0 - zalf
-            ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj)
-            ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) )
-            psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
-            psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   &
-               &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          &
-               &                                         + ( zalf1 - zalf ) * ztemp )                                )
-            psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       &
-               &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  )
-            psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
-            psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
-         END DO
-      END DO
-
-      !-- Lateral boundary conditions
-      CALL lbc_lnk_multi( psm , 'T',  1.,  ps0 , 'T',  1.   &
-         &              , psx , 'T', -1.,  psy , 'T', -1.   &   ! caution gradient ==> the sign changes
-         &              , psxx, 'T',  1.,  psyy, 'T',  1.   &
-         &              , psxy, 'T',  1. )
-
-      IF(ln_ctl) THEN
-         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
-         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
-         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
-         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
-      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
-      !
-   END SUBROUTINE lim_adv_y
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option            Dummy module         NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-#endif
-
-   !!======================================================================
-END MODULE limadv
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_prather.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_prather.F90	(revision 6989)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_prather.F90	(revision 6989)
@@ -0,0 +1,418 @@
+MODULE limadv_prather 
+   !!======================================================================
+   !!                       ***  MODULE limadv_prather   ***
+   !! LIM sea-ice model : sea-ice advection => Prather scheme
+   !!======================================================================
+   !! History :  LIM  ! 2008-03 (M. Vancoppenolle)  LIM-3 from LIM-2 code
+   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c.
+   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!--------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                     LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   lim_adv_x  : advection of sea ice on x axis
+   !!   lim_adv_y  : advection of sea ice on y axis
+   !!----------------------------------------------------------------------
+   USE dom_oce          ! ocean domain
+   USE ice              ! LIM-3 variables
+   USE lbclnk           ! lateral boundary condition - MPP exchanges
+   USE in_out_manager   ! I/O manager
+   USE prtctl           ! Print control
+   USE lib_mpp          ! MPP library
+   USE wrk_nemo         ! work arrays
+   USE lib_fortran      ! to use key_nosignedzero
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   lim_adv_x   ! called by lim_trp
+   PUBLIC   lim_adv_y   ! called by lim_trp
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
+   !! $Id: limadv.F90 6746 2016-06-27 17:20:57Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE lim_adv_x( pdf, put , pcrh, psm , ps0 ,   &
+      &                  psx, psxx, psy , psyy, psxy )
+      !!---------------------------------------------------------------------
+      !!                **  routine lim_adv_x  **
+      !!  
+      !! ** purpose :   Computes and adds the advection trend to sea-ice
+      !!              variable on i-axis
+      !!
+      !! ** method  :   Uses Prather second order scheme that advects tracers
+      !!              but also theirquadratic forms. The method preserves
+      !!              tracer structures by conserving second order moments.
+      !!
+      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
+      !!--------------------------------------------------------------------
+      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
+      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
+      !! 
+      INTEGER  ::   ji, jj                               ! dummy loop indices
+      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars
+      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
+      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zalg, zalg1, zalg1q         !  -      -
+      !---------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
+      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
+
+      ! Limitation of moments.                                           
+
+      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zslpmax = MAX( 0._wp, ps0(ji,jj) )
+            zs1max  = 1.5 * zslpmax
+            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
+            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
+               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
+            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
+
+            ps0 (ji,jj) = zslpmax  
+            psx (ji,jj) = zs1new      * rswitch
+            psxx(ji,jj) = zs2new      * rswitch
+            psy (ji,jj) = psy (ji,jj) * rswitch
+            psyy(ji,jj) = psyy(ji,jj) * rswitch
+            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
+         END DO
+      END DO
+
+      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
+      psm (:,:)  = MAX( pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
+
+      !  Calculate fluxes and moments between boxes i<-->i+1              
+      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0 
+         DO ji = 1, jpi
+            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
+            zalf         =  MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
+            zalfq        =  zalf * zalf
+            zalf1        =  1.0 - zalf
+            zalf1q       =  zalf1 * zalf1
+            !
+            zfm (ji,jj)  =  zalf  *   psm (ji,jj)
+            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) )  )
+            zfx (ji,jj)  =  zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
+            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj) * zalfq
+            zfy (ji,jj)  =  zalf  * ( psy (ji,jj) + zalf1 * psxy(ji,jj) )
+            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj)
+            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj)
+
+            !  Readjust moments remaining in the box.
+            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
+            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
+            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
+            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
+            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
+            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
+            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
+         END DO
+      END DO
+
+      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
+         DO ji = 1, fs_jpim1
+            zalf          = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
+            zalg  (ji,jj) = zalf
+            zalfq         = zalf * zalf
+            zalf1         = 1.0 - zalf
+            zalg1 (ji,jj) = zalf1
+            zalf1q        = zalf1 * zalf1
+            zalg1q(ji,jj) = zalf1q
+            !
+            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji+1,jj)
+            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
+            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
+            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji+1,jj) * zalfq
+            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) )
+            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji+1,jj)
+            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji+1,jj)
+         END DO
+      END DO
+
+      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box. 
+         DO ji = fs_2, fs_jpim1
+            zbt  =       zbet(ji-1,jj)
+            zbt1 = 1.0 - zbet(ji-1,jj)
+            !
+            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
+            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
+            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
+            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
+            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
+            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
+            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
+         END DO
+      END DO
+
+      !   Put the temporary moments into appropriate neighboring boxes.    
+      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
+         DO ji = fs_2, fs_jpim1
+            zbt  =       zbet(ji-1,jj)
+            zbt1 = 1.0 - zbet(ji-1,jj)
+            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
+            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
+            zalf1       = 1.0 - zalf
+            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
+            !
+            ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj)
+            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
+            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               &
+               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
+               &                                                + zbt1 * psxx(ji,jj)
+            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
+               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
+               &                                                + zbt1 * psxy(ji,jj)
+            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
+            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
+         END DO
+      END DO
+
+      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
+         DO ji = fs_2, fs_jpim1
+            zbt  =       zbet(ji,jj)
+            zbt1 = 1.0 - zbet(ji,jj)
+            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
+            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
+            zalf1       = 1.0 - zalf
+            ztemp       = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
+            !
+            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
+            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
+            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  &
+               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )      &
+               &                                      + ( zalf1 - zalf ) * ztemp ) )
+            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  &
+               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
+            psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
+            psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
+         END DO
+      END DO
+
+      !-- Lateral boundary conditions
+      CALL lbc_lnk_multi( psm , 'T',  1., ps0 , 'T',  1.   &
+         &              , psx , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
+         &              , psxx, 'T',  1., psyy, 'T',  1.   &
+         &              , psxy, 'T',  1. )
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
+         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
+         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
+         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
+      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
+      !
+   END SUBROUTINE lim_adv_x
+
+
+   SUBROUTINE lim_adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
+      &                  psx, psxx, psy , psyy, psxy )
+      !!---------------------------------------------------------------------
+      !!                **  routine lim_adv_y  **
+      !!            
+      !! ** purpose :   Computes and adds the advection trend to sea-ice 
+      !!              variable on y axis
+      !!
+      !! ** method  :   Uses Prather second order scheme that advects tracers
+      !!              but also their quadratic forms. The method preserves 
+      !!              tracer structures by conserving second order moments.
+      !! 
+      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
+      !!---------------------------------------------------------------------
+      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
+      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_adv_y (=1) or the opposite (=0)
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
+      !!
+      INTEGER  ::   ji, jj                               ! dummy loop indices
+      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars
+      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
+      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zf0, zfx , zfy , zbet   ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zfm, zfxx, zfyy, zfxy   !  -      -
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zalg, zalg1, zalg1q     !  -      -
+      !---------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
+      CALL wrk_alloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
+
+      ! Limitation of moments.
+
+      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zslpmax = MAX( 0._wp, ps0(ji,jj) )
+            zs1max  = 1.5 * zslpmax
+            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
+            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
+               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
+            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
+            !
+            ps0 (ji,jj) = zslpmax  
+            psx (ji,jj) = psx (ji,jj) * rswitch
+            psxx(ji,jj) = psxx(ji,jj) * rswitch
+            psy (ji,jj) = zs1new * rswitch
+            psyy(ji,jj) = zs2new * rswitch
+            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
+         END DO
+      END DO
+
+      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
+      psm(:,:)  = MAX(  pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
+
+      !  Calculate fluxes and moments between boxes j<-->j+1              
+      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
+         DO ji = 1, jpi
+            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
+            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
+            zalfq        =  zalf * zalf
+            zalf1        =  1.0 - zalf
+            zalf1q       =  zalf1 * zalf1
+            !
+            zfm (ji,jj)  =  zalf  * psm(ji,jj)
+            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
+            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
+            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
+            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
+            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
+            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
+            !
+            !  Readjust moments remaining in the box.
+            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
+            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
+            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
+            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
+            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
+            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
+            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
+         END DO
+      END DO
+      !
+      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
+         DO ji = 1, jpi
+            zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
+            zalg  (ji,jj) = zalf
+            zalfq         = zalf * zalf
+            zalf1         = 1.0 - zalf
+            zalg1 (ji,jj) = zalf1
+            zalf1q        = zalf1 * zalf1
+            zalg1q(ji,jj) = zalf1q
+            !
+            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1)
+            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
+            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
+            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq
+            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) )
+            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1)
+            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1)
+         END DO
+      END DO
+
+      !  Readjust moments remaining in the box. 
+      DO jj = 2, jpj
+         DO ji = 1, jpi
+            zbt  =         zbet(ji,jj-1)
+            zbt1 = ( 1.0 - zbet(ji,jj-1) )
+            !
+            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
+            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
+            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
+            psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj)
+            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
+            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
+            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
+         END DO
+      END DO
+
+      !   Put the temporary moments into appropriate neighboring boxes.    
+      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
+         DO ji = 1, jpi
+            zbt  =         zbet(ji,jj-1)
+            zbt1 = ( 1.0 - zbet(ji,jj-1) )
+            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
+            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
+            zalf1       = 1.0 - zalf
+            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
+            !
+            ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj)
+            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
+               &                                               + zbt1 * psy(ji,jj)  
+            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
+               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
+               &                                               + zbt1 * psyy(ji,jj)
+            psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               &
+               &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   &
+               &                                                + zbt1 * psxy(ji,jj)
+            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
+            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
+         END DO
+      END DO
+
+      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
+         DO ji = 1, jpi
+            zbt  =         zbet(ji,jj)
+            zbt1 = ( 1.0 - zbet(ji,jj) )
+            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
+            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
+            zalf1       = 1.0 - zalf
+            ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj)
+            ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) )
+            psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
+            psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   &
+               &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          &
+               &                                         + ( zalf1 - zalf ) * ztemp )                                )
+            psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       &
+               &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  )
+            psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
+            psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
+         END DO
+      END DO
+
+      !-- Lateral boundary conditions
+      CALL lbc_lnk_multi( psm , 'T',  1.,  ps0 , 'T',  1.   &
+         &              , psx , 'T', -1.,  psy , 'T', -1.   &   ! caution gradient ==> the sign changes
+         &              , psxx, 'T',  1.,  psyy, 'T',  1.   &
+         &              , psxy, 'T',  1. )
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
+         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
+         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
+         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, zf0 , zfx , zfy , zbet, zfm )
+      CALL wrk_dealloc( jpi, jpj, zfxx, zfyy, zfxy, zalg, zalg1, zalg1q )
+      !
+   END SUBROUTINE lim_adv_y
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option            Dummy module         NO LIM sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+
+   !!======================================================================
+END MODULE limadv_prather
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90	(revision 6988)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90	(revision 6989)
@@ -29,8 +29,7 @@
 
    PUBLIC   lim_adv_umx    ! routine called by limtrp.F90
-   
-   INTEGER  ::   nn_ult_order = 4       ! order of the ULTIMATE scheme
-   
-   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
+      
+   REAL(wp) ::   z1_6   = 1._wp / 6._wp   ! =1/6
+   REAL(wp) ::   z1_120 = 1._wp / 120._wp ! =1/120
 
    !! * Substitutions
@@ -43,5 +42,5 @@
 CONTAINS
 
-   SUBROUTINE lim_adv_umx( lcon, kt, pdt, pu, pv, puc, pvc, pt, ptc, pfu, pfv )
+   SUBROUTINE lim_adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE lim_adv_umx  ***
@@ -56,13 +55,9 @@
       !! ** Action : - pt  the after advective tracer
       !!----------------------------------------------------------------------
-      LOGICAL                     , INTENT(in   )           ::   lcon       ! "continuity" equations for a and H 
       INTEGER                     , INTENT(in   )           ::   kt         ! number of iteration
       REAL(wp)                    , INTENT(in   )           ::   pdt        ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pu, pv     ! 2 ice velocity components
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc, pvc   ! 2 ice velocity components  (for a or H)
-                                                                            ! 2 ice transport components (for tracers q)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pt         ! tracer field
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc, pvc   ! 2 ice velocity components => u*e2
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc        ! tracer content field
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out), OPTIONAL ::   pfu, pfv   ! advective fluxes at u- and v-points 
       !
       INTEGER  ::   ji, jj           ! dummy loop indices  
@@ -77,6 +72,4 @@
       CALL wrk_alloc( jpi,jpj,   zt_ups, zfu_ups, zfv_ups, ztrd, zfu_ho, zfv_ho, zt_u, zt_v )
       !
-      !
-!clem      zt_ups(:,:) = 0._wp
       !
       !  upstream advection with initial mass fluxes & intermediate update
@@ -88,6 +81,6 @@
             zfp_vj = pvc(ji,jj) + ABS( pvc(ji,jj) )
             zfm_vj = pvc(ji,jj) - ABS( pvc(ji,jj) )
-            zfu_ups(ji,jj) = 0.5_wp * ( zfp_ui * pt(ji,jj) + zfm_ui * pt(ji+1,jj  ) )
-            zfv_ups(ji,jj) = 0.5_wp * ( zfp_vj * pt(ji,jj) + zfm_vj * pt(ji  ,jj+1) )
+            zfu_ups(ji,jj) = 0.5_wp * ( zfp_ui * ptc(ji,jj) + zfm_ui * ptc(ji+1,jj  ) )
+            zfv_ups(ji,jj) = 0.5_wp * ( zfp_vj * ptc(ji,jj) + zfm_vj * ptc(ji  ,jj+1) )
          END DO
       END DO
@@ -106,16 +99,15 @@
       ! High order (_ho) fluxes 
       ! -----------------------
-      SELECT CASE( nn_ult_order )
+      SELECT CASE( nn_limadv_ord )
       CASE ( 20 )                          ! centered second order
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) )
-               zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) )
-            END DO
-         END DO
-         !
-      CASE ( 1 , 4 )                      ! 1st to 4th order ULTIMATE-MACHO scheme
-         CALL macho( lcon, kt, nn_ult_order, pdt, pt, pu, pv, zt_u, zt_v )
-!!         CALL macho( lcon, kt, nn_ult_order, pdt, ptc, pu, pv, zt_u, zt_v )
+               zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( ptc(ji,jj) + ptc(ji+1,jj) )
+               zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj) + ptc(ji,jj+1) )
+            END DO
+         END DO
+         !
+      CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme
+         CALL macho( kt, nn_limadv_ord, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v )
          !
          DO jj = 2, jpjm1
@@ -153,13 +145,4 @@
       CALL lbc_lnk( ptc(:,:) , 'T',  1. )
       !
-      IF( PRESENT( pfu ) ) THEN
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               pfu(ji,jj) = zfu_ups(ji,jj) + zfu_ho(ji,jj)
-               pfv(ji,jj) = zfv_ups(ji,jj) + zfv_ho(ji,jj)
-            END DO
-         END DO
-         CALL lbc_lnk_multi( pfu, 'U', -1., pfv, 'V', -1. )         ! Lateral bondary conditions
-      ENDIF
       !
       CALL wrk_dealloc( jpi,jpj,   zt_ups, zfu_ups, zfv_ups, ztrd, zfu_ho, zfv_ho, zt_u, zt_v )
@@ -170,5 +153,5 @@
 
 
-   SUBROUTINE macho( lcon, kt, k_order, pdt, pt, pu, pv, pt_u, pt_v )
+   SUBROUTINE macho( kt, k_order, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_x  ***
@@ -181,10 +164,10 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
-      LOGICAL                     , INTENT(in   ) ::   lcon       ! "continuity" equations for a and ah 
       INTEGER                     , INTENT(in   ) ::   kt         ! number of iteration
       INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme
       REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu, pv     ! 2 ice velocity components
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt         ! tracer fields
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc        ! tracer fields
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc, pvc   ! 2 ice velocity components
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points 
       !
@@ -201,48 +184,36 @@
          !
          !                                                           !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( lcon, k_order, pdt, pt, pu, pt_u )
-         !
-         !                                                           !--  advective form update in zzt  --!
-         DO jj = 1, jpj
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               !
-               IF( pu(ji,jj) * pu(ji-1,jj) <= 0._wp ) THEN   ;   zc_box = 0._wp
-               ELSEIF(         pu(ji  ,jj) >  0._wp ) THEN   ;   zc_box = pu(ji-1,jj)
-               ELSE                                          ;   zc_box = pu(ji  ,jj)
-               ENDIF
-               !
-               zzt(ji,jj) = pt(ji,jj) - zc_box * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e12t(ji,jj)
-               IF( lcon .eqv. .TRUE. ) THEN ! clem => for a.div(u) ??
-                  zzt(ji,jj) = zzt(ji,jj) - pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e12t(ji,jj)
-               END IF
-            END DO
-         END DO
-         !
-         !                                                           !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( lcon, k_order, pdt, zzt, pv, pt_v )
-         !
-      ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==!
-         !
-         !                                                           !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( lcon, k_order, pdt, pt, pv, pt_v )
+         CALL ultimate_x( k_order, pdt, ptc, puc, pt_u )
          !
          !                                                           !--  advective form update in zzt  --!
          DO jj = 2, jpjm1
-            DO ji = 1, jpi
-               !
-               IF( pv(ji,jj) * pv(ji,jj-1) <= 0._wp ) THEN   ;   zc_box = 0._wp
-               ELSEIF(         pv(ji,jj  ) >  0._wp ) THEN   ;   zc_box = pv(ji,jj-1)
-               ELSE                                          ;   zc_box = pv(ji,jj  )
-               ENDIF
-               !
-               zzt(ji,jj) = pt(ji,jj) - zc_box * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e12t(ji,jj)
-               IF( lcon .eqv. .TRUE. ) THEN ! clem => for a.div(u) ??
-                  zzt(ji,jj) = zzt(ji,jj) - pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e12t(ji,jj)
-               END IF
-            END DO
-         END DO
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj)  &
+                  &                    - ptc  (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e12t(ji,jj)
+               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
+            END DO
+         END DO
+         CALL lbc_lnk( zzt, 'T', 1. )
+         !
+         !                                                           !--  ultimate interpolation of pt at v-point  --!
+         CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v )
+         !
+      ELSE                                                  !==  even ice time step:  adv_y then adv_x  ==!
+         !
+         !                                                           !--  ultimate interpolation of pt at v-point  --!
+         CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v )
+         !
+         !                                                           !--  advective form update in zzt  --!
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj)  &
+                  &                    - ptc  (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e12t(ji,jj)
+               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
+            END DO
+         END DO
+         CALL lbc_lnk( zzt, 'T', 1. )
          !
          !                                                           !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( lcon, k_order, pdt, zzt, pu, pt_u )
+         CALL ultimate_x( k_order, pdt, zzt, puc, pt_u )
          !      
       ENDIF      
@@ -255,5 +226,5 @@
 
 
-   SUBROUTINE ultimate_x( lcon, k_order, pdt, pt, pu, pt_u )
+   SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_x  ***
@@ -266,19 +237,43 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
-      LOGICAL                     , INTENT(in   ) ::   lcon      ! "continuity" equations for a and ah 
       INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index
       REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pu        ! ice i-velocity component
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc       ! ice i-velocity component
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_u      ! tracer at u-point 
       !
       INTEGER  ::   ji, jj       ! dummy loop indices
-      REAL(wp) ::   zcu, zdx2    !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:) :: ztu, zltu
+      REAL(wp) ::   zcu, zdx2, zdx4    !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:) :: ztu1, ztu2, ztu3, ztu4
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('ultimate_x')
       !
-      CALL wrk_alloc( jpi,jpj,   ztu, zltu )
+      CALL wrk_alloc( jpi,jpj,   ztu1, ztu2, ztu3, ztu4 )
+      !
+      !                                                     !--  Laplacian in i-direction  --!
+      DO jj = 2, jpjm1         ! First derivative (gradient)
+         DO ji = 1, fs_jpim1
+            ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
+         END DO
+         !                     ! Second derivative (Laplacian)
+         DO ji = fs_2, fs_jpim1
+            ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj)
+         END DO
+      END DO
+      CALL lbc_lnk( ztu2, 'T', 1. )
+      !
+      !                                                     !--  BiLaplacian in i-direction  --!
+      DO jj = 2, jpjm1         ! Third derivative
+         DO ji = 1, fs_jpim1
+            ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
+         END DO
+         !                     ! Fourth derivative
+         DO ji = fs_2, fs_jpim1
+            ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj)
+         END DO
+      END DO
+      CALL lbc_lnk( ztu4, 'T', 1. )
+      !
       !
       SELECT CASE (k_order )
@@ -288,15 +283,16 @@
          DO jj = 1, jpj
             DO ji = 1, fs_jpim1   ! vector opt.
-               pt_u(ji,jj) = 0.5_wp * (                             ( pt(ji+1,jj) + pt(ji,jj) )    &
-                  &                    - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) )  ) * umask(ji,jj,1)
+               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj) + pt(ji,jj)   &
+                  &                                    - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) )
             END DO
          END DO
          !
       CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
+         !
          DO jj = 1, jpj
             DO ji = 1, fs_jpim1   ! vector opt.
-               pt_u(ji,jj) = 0.5_wp * (                    ( pt(ji+1,jj) + pt(ji,jj) )                 &
-!!                  &                    - pdt * pu(ji,jj) * ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj)  ) * umask(ji,jj,1)
-                  &                    - pdt * pu(ji,jj) * ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e12u(ji,jj) ) * umask(ji,jj,1)
+               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
+               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (                                   pt(ji+1,jj) + pt(ji,jj)   &
+                  &                                               -              zcu   * ( pt(ji+1,jj) - pt(ji,jj) ) ) 
             END DO
          END DO
@@ -305,51 +301,44 @@
       CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
          !
-         !                                                                 !--  Laplacian in i- and in j-directions  --!
-         DO jj = 2, jpjm1         ! First derivative (gradient)
-            DO ji = 1, fs_jpim1   ! vector opt.
-               ztu(ji,jj) = ( pt(ji+1,jj  ) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
-            END DO
-            !                     ! Second derivative (divergence)
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zltu(ji,jj) = (  ztu(ji,jj) - ztu(ji-1,jj)  ) * r1_e1t(ji,jj)
-            END DO
-         END DO
-         CALL lbc_lnk( zltu, 'T', 1. )   ! Lateral boundary cond. (unchanged sign)
-         !    
          DO jj = 1, jpj
             DO ji = 1, fs_jpim1   ! vector opt.
-!!               zcu  = pu(ji,jj) * pdt * r1_e1u(ji,jj)
-               zcu  = pu(ji,jj) * pdt * r1_e12u(ji,jj)
+               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
+!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
                pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                         pt  (ji+1,jj) + pt  (ji,jj)        &
                   &                                               -              zcu   * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
-                  &        - r1_6 * zdx2 * ( 1._wp - zcu*zcu ) * (                         zltu(ji+1,jj) + zltu(ji,jj)        &
-                  &                                               - SIGN( 1._wp, zcu ) * ( zltu(ji+1,jj) - zltu(ji,jj) )  )   )
+                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                         ztu2(ji+1,jj) + ztu2(ji,jj)        &
+                  &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
             END DO
          END DO
          !
       CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
-         !
-         !                                                                 !--  Laplacian in i- and in j-directions  --!
-         DO jj = 1, jpjm1         ! First derivative (gradient)
-            DO ji = 1, fs_jpim1   ! vector opt.
-               ztu(ji,jj) = ( pt(ji+1,jj  ) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
-            END DO
-            !                     ! Second derivative (divergence)
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zltu(ji,jj) = (  ztu(ji,jj) - ztu(ji-1,jj)  ) * r1_e1t(ji,jj)
-            END DO
-         END DO
-         CALL lbc_lnk( zltu, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
          !
          DO jj = 1, jpj
             DO ji = 1, fs_jpim1   ! vector opt.
-!!               zcu  = pu(ji,jj) * pdt * r1_e1u(ji,jj)
-               zcu  = pu(ji,jj) * pdt * r1_e12u(ji,jj)
+               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
+!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
                pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (         (                   pt  (ji+1,jj) + pt  (ji,jj)        &
                   &                                               -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) )  )   &
-                  &        - r1_6 * zdx2 * ( 1._wp - zcu*zcu ) * (                   zltu(ji+1,jj) + zltu(ji,jj)        &
-                  &                                               - 0.5_wp * zcu * ( zltu(ji+1,jj) - zltu(ji,jj) )  )   )
+                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                   ztu2(ji+1,jj) + ztu2(ji,jj)        &
+                  &                                               - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
+            END DO
+         END DO
+         !
+      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
+         !
+         DO jj = 1, jpj
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
+               zdx2 = e1u(ji,jj) * e1u(ji,jj)
+!!rachid       zdx2 = e1u(ji,jj) * e1t(ji,jj)
+               zdx4 = zdx2 * zdx2
+               pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * (               (                   pt  (ji+1,jj) + pt  (ji,jj)       &
+                  &                                                     -          zcu * ( pt  (ji+1,jj) - pt  (ji,jj) ) )   &
+                  &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *     (                   ztu2(ji+1,jj) + ztu2(ji,jj)       &
+                  &                                                     - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) )   &
+                  &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj)       &
+                  &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) )
             END DO
          END DO
@@ -357,5 +346,5 @@
       END SELECT
       !
-      CALL wrk_dealloc( jpi,jpj,   ztu, zltu )
+      CALL wrk_dealloc( jpi,jpj,   ztu1, ztu2, ztu3, ztu4 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('ultimate_x')
@@ -364,5 +353,5 @@
    
  
-   SUBROUTINE ultimate_y( lcon, k_order, pdt, pt, pv, pt_v )
+   SUBROUTINE ultimate_y( k_order, pdt, pt, pvc, pt_v )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_y  ***
@@ -375,19 +364,45 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
-      LOGICAL                     , INTENT(in   ) ::   lcon      ! "continuity" equations for a and ah 
       INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index
       REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pv        ! ice j-velocity component
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvc       ! ice j-velocity component
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pt        ! tracer fields
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pt_v      ! tracer at v-point 
       !
       INTEGER  ::   ji, jj       ! dummy loop indices
-      REAL(wp) ::   zcv, zdy2    !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:) :: ztv, zltv
+      REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:) :: ztv1, ztv2, ztv3, ztv4
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('ultimate_y')
       !
-      CALL wrk_alloc( jpi,jpj,   ztv, zltv )
+      CALL wrk_alloc( jpi,jpj,   ztv1, ztv2, ztv3, ztv4 )
+      !
+      !                                                     !--  Laplacian in j-direction  --!
+      DO jj = 1, jpjm1         ! First derivative (gradient)
+         DO ji = fs_2, fs_jpim1
+            ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
+         END DO
+      END DO
+      DO jj = 2, jpjm1         ! Second derivative (Laplacian)
+         DO ji = fs_2, fs_jpim1
+            ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj)
+         END DO
+      END DO
+      CALL lbc_lnk( ztv2, 'T', 1. )
+      !
+      !                                                     !--  BiLaplacian in j-direction  --!
+      DO jj = 1, jpjm1         ! First derivative
+         DO ji = fs_2, fs_jpim1
+            ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
+         END DO
+      END DO
+      DO jj = 2, jpjm1         ! Second derivative
+         DO ji = fs_2, fs_jpim1
+            ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj)
+         END DO
+      END DO
+      CALL lbc_lnk( ztv4, 'T', 1. )
+      !
       !
       SELECT CASE (k_order )
@@ -397,6 +412,6 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpi
-               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                             ( pt(ji,jj+1) + pt(ji,jj) )  &
-                  &                                     - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) )  )
+               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1) + pt(ji,jj) )  &
+                  &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) )
             END DO
          END DO
@@ -405,7 +420,7 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpi
-               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                    ( pt(ji,jj+1) + pt(ji,jj) )               &
-!!                  &                                     - pdt * pv(ji,jj) * ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj)  )
-                  &                                     - pdt * pv(ji,jj) * ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e12v(ji,jj) )
+               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (        ( pt(ji,jj+1) + pt(ji,jj) )  &
+                  &                                     - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) )
             END DO
          END DO
@@ -413,54 +428,45 @@
          !
       CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
-         !
-         !                                                                 !--  Laplacian in i- and in j-directions  --!
-         DO jj = 1, jpjm1            ! First derivative (gradient)
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ztv(ji,jj) = ( pt(ji  ,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
-            END DO
-         END DO
-         DO jj = 2, jpjm1            ! Second derivative (divergence)
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zltv(ji,jj) = (  ztv(ji,jj) - ztv(ji,jj-1)  ) * r1_e2t(ji,jj)
-            END DO
-         END DO
-         CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sign)
          !
          DO jj = 1, jpjm1
             DO ji = 1, jpi
-!!               zcv  = pv(ji,jj) * pdt * r1_e2v(ji,jj)
-               zcv  = pv(ji,jj) * pdt * r1_e12v(ji,jj)
+               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
-               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                (  pt  (ji,jj+1) + pt  (ji,jj)        &
-                  &                                     -                        zcv   * ( pt  (ji,jj+1) - pt  (ji,jj) )  )   &
-                  &        - r1_6 * zdy2 * ( 1._wp - zcv*zcv ) * (                         zltv(ji,jj+1) + zltv(ji,jj)        &
-                  &                                               - SIGN( 1._wp, zcv ) * ( zltv(ji,jj+1) - zltv(ji,jj) )  )   )
+!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
+               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)       &
+                  &                                     -                        zcv   * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
+                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1) + ztv2(ji,jj)       &
+                  &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
             END DO
          END DO
          !
       CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
-         !
-         !                                                                 !--  Laplacian in i- and in j-directions  --!
-         DO jj = 1, jpjm1            ! First derivative (gradient)
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ztv(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
-            END DO
-         END DO
-         DO jj = 2, jpjm1            ! Second derivative (divergence)
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zltv(ji,jj) = (  ztv(ji,jj) - ztv(ji,jj-1)  ) * r1_e2t(ji,jj)
-            END DO
-         END DO
-         CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
          !
          DO jj = 1, jpjm1
             DO ji = 1, jpi
-!!               zcv  = pv(ji,jj) * pdt * r1_e2v(ji,jj)
-               zcv  = pv(ji,jj) * pdt * r1_e12v(ji,jj)
+               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
-               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                          (  pt  (ji,jj+1) + pt  (ji,jj)        &
-                  &                                               -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) )  )   &
-                  &        - r1_6 * zdy2 * ( 1._wp - zcv*zcv ) * (                   zltv(ji,jj+1) + zltv(ji,jj)        &
-                  &                                               - 0.5_wp * zcv * ( zltv(ji,jj+1) - zltv(ji,jj) )  )   )
+!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
+               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt  (ji,jj+1) + pt  (ji,jj)       &
+                  &                                               -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )   &
+                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1) + ztv2(ji,jj)       &
+                  &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) )
+            END DO
+         END DO
+         !
+      CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
+         !
+         DO jj = 1, jpjm1
+            DO ji = 1, jpi
+               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+               zdy2 = e2v(ji,jj) * e2v(ji,jj)
+!!rachid       zdy2 = e2v(ji,jj) * e2t(ji,jj)
+               zdy4 = zdy2 * zdy2
+               pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * (                                 ( pt  (ji,jj+1) + pt  (ji,jj)      &
+                  &                                                     -          zcv * ( pt  (ji,jj+1) - pt  (ji,jj) ) )  &
+                  &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1) + ztv2(ji,jj)      &
+                  &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) )  &
+                  &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj)      &
+                  &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) )
             END DO
          END DO
@@ -468,5 +474,5 @@
       END SELECT
       !
-      CALL wrk_dealloc( jpi,jpj,   ztv, zltv )
+      CALL wrk_dealloc( jpi,jpj,   ztv1, ztv2, ztv3, ztv4 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('ultimate_y')
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 6988)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 6989)
@@ -26,6 +26,6 @@
    USE prtctl           ! Print control
    USE lib_fortran      ! glob_sum
-   USE timing          ! Timing
-   USE limcons        ! conservation tests
+   USE timing           ! Timing
+   USE limcons          ! conservation tests
    USE limvar
 
@@ -33,5 +33,6 @@
    PRIVATE
 
-   PUBLIC   lim_dyn   ! routine called by ice_step
+   PUBLIC   lim_dyn        ! routine called by sbcice_lim.F90
+   PUBLIC   lim_dyn_init   ! routine called by sbcice_lim.F90
 
    !! * Substitutions
@@ -64,6 +65,4 @@
 
       CALL lim_var_agg(1)                      ! aggregate ice categories
-
-      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
       !
       ! conservation test
@@ -147,5 +146,6 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, &
+      NAMELIST/namicedyn/ nn_limadv, nn_limadv_ord,  &
+         &                nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, &
          &                nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr, rn_lfrelax
       !!-------------------------------------------------------------------
@@ -164,4 +164,8 @@
          WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
          WRITE(numout,*) '~~~~~~~~~~~~'
+         ! limtrp
+         WRITE(numout,*)'    choose the advection scheme (-1=Prather, 0=Ulimate-Macho)   nn_limadv     = ', nn_limadv 
+         WRITE(numout,*)'    choose the order of the scheme (if ultimate)                nn_limadv_ord = ', nn_limadv_ord  
+         ! limrhg
          WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)         nn_icestr     = ', nn_icestr 
          WRITE(numout,*)'    Including brine volume in ice strength comp.                ln_icestr_bvf = ', ln_icestr_bvf
@@ -180,7 +184,4 @@
       ENDIF
       !
-      usecc2 = 1._wp / ( rn_ecc * rn_ecc )
-      rhoco  = rau0  * rn_cio
-      !
    END SUBROUTINE lim_dyn_init
 
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 6988)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 6989)
@@ -114,4 +114,5 @@
       CHARACTER (len=50) ::   charout
 
+      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio
       REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling
       REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity
@@ -203,4 +204,6 @@
       ! 1) define some variables and initialize arrays
       !------------------------------------------------------------------------------!
+      zrhoco = rau0 * rn_cio 
+
       ! ecc2: square of yield ellipse eccenticrity
       ecc2    = rn_ecc * rn_ecc
@@ -357,5 +360,5 @@
                
                ! delta at T points
-               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )  
+               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
 
                ! P/delta at T points
@@ -424,6 +427,6 @@
 
                   ! tau_io/(v_oce - v_ice)
-                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
-                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
 
                   ! tau_bottom/v_ice
@@ -469,6 +472,6 @@
                                
                   ! tau_io/(u_oce - u_ice)
-                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
-                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
 
                   ! tau_bottom/u_ice
@@ -515,6 +518,6 @@
 
                   ! tau_io/(u_oce - u_ice)
-                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
-                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
 
                   ! tau_bottom/u_ice
@@ -559,6 +562,6 @@
          
                   ! tau_io/(v_oce - v_ice)
-                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
-                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
 
                   ! tau_bottom/v_ice
@@ -651,5 +654,5 @@
             
             ! delta at T points
-            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )  
+            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
             rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
             delta_i(ji,jj) = zdelta + rn_creepl * rswitch
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 6988)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 6989)
@@ -53,6 +53,6 @@
       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
+      CHARACTER(len=20)   ::   clkt     ! ocean time-step define as a character
+      CHARACTER(len=50)   ::   clname   ! ice output restart file name
       CHARACTER(len=256)  ::   clpath   ! full path to ice output restart file 
       !!----------------------------------------------------------------------
@@ -107,6 +107,6 @@
       INTEGER ::   ji, jj, jk ,jl   ! dummy loop indices
       INTEGER ::   iter
-      CHARACTER(len=15) ::   znam
-      CHARACTER(len=1)  ::   zchar, zchar1
+      CHARACTER(len=25) ::   znam
+      CHARACTER(len=2)  ::   zchar, zchar1
       REAL(wp), POINTER, DIMENSION(:,:) :: z2d
       !!----------------------------------------------------------------------
@@ -130,5 +130,5 @@
       ! Prognostic variables 
       DO jl = 1, jpl 
-         WRITE(zchar,'(I1)') jl
+         WRITE(zchar,'(I2.2)') jl
          znam = 'v_i'//'_htc'//zchar
          z2d(:,:) = v_i(:,:,jl)
@@ -152,5 +152,5 @@
 
       DO jl = 1, jpl 
-         WRITE(zchar,'(I1)') jl
+         WRITE(zchar,'(I2.2)') jl
          znam = 'tempt_sl1'//'_htc'//zchar
          z2d(:,:) = e_s(:,:,1,jl)
@@ -159,7 +159,7 @@
 
       DO jl = 1, jpl 
-         WRITE(zchar,'(I1)') jl
+         WRITE(zchar,'(I2.2)') jl
          DO jk = 1, nlay_i 
-            WRITE(zchar1,'(I1)') jk
+            WRITE(zchar1,'(I2.2)') jk
             znam = 'tempt'//'_il'//zchar1//'_htc'//zchar
             z2d(:,:) = e_i(:,:,jk,jl)
@@ -176,128 +176,136 @@
       CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b )
 
-      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
-
+      ! In case Prather scheme is used for advection, write second order moments
+      ! ------------------------------------------------------------------------
+      IF( nn_limadv == -1 ) THEN
+         
+         DO jl = 1, jpl 
+            WRITE(zchar,'(I2.2)') 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,'(I2.2)') jl
+            DO jk = 1, nlay_i 
+               WRITE(zchar1,'(I2.2)') 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
+
+      ENDIF
+      
+      ! close restart file
+      ! ------------------
       IF( iter == nitrst ) THEN
-         CALL iom_close( numriw )                         ! close the restart file
+         CALL iom_close( numriw )
          lrst_ice = .FALSE.
       ENDIF
@@ -317,6 +325,6 @@
       REAL(wp) ::   zfice, ziter
       REAL(wp), POINTER, DIMENSION(:,:) ::   z2d
-      CHARACTER(len=15) ::   znam
-      CHARACTER(len=1)  ::   zchar, zchar1
+      CHARACTER(len=25) ::   znam
+      CHARACTER(len=2)  ::   zchar, zchar1
       INTEGER           ::   jlibalt = jprstlib
       LOGICAL           ::   llok
@@ -356,6 +364,7 @@
          &                   '   control of time parameter  nrstdt' )
 
+      ! Prognostic variables 
       DO jl = 1, jpl 
-         WRITE(zchar,'(I1)') jl
+         WRITE(zchar,'(I2.2)') jl
          znam = 'v_i'//'_htc'//zchar
          CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
@@ -379,5 +388,5 @@
 
       DO jl = 1, jpl 
-         WRITE(zchar,'(I1)') jl
+         WRITE(zchar,'(I2.2)') jl
          znam = 'tempt_sl1'//'_htc'//zchar
          CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
@@ -386,7 +395,7 @@
 
       DO jl = 1, jpl 
-         WRITE(zchar,'(I1)') jl
+         WRITE(zchar,'(I2.2)') jl
          DO jk = 1, nlay_i 
-            WRITE(zchar1,'(I1)') jk
+            WRITE(zchar1,'(I2.2)') jk
             znam = 'tempt'//'_il'//zchar1//'_htc'//zchar
             CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
@@ -403,126 +412,132 @@
       CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b )
 
-      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
-      !
+      ! In case Prather scheme is used for advection, read second order moments
+      ! ------------------------------------------------------------------------
+      IF( nn_limadv == -1 ) THEN
+
+         DO jl = 1, jpl 
+            WRITE(zchar,'(I2.2)') 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,'(I2.2)') jl
+            DO jk = 1, nlay_i 
+               WRITE(zchar1,'(I2.2)') 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
+         !
+      END IF
+      
       ! clem: I do not understand why the following IF is needed
       !       I suspect something inconsistent in the main code with option nn_icesal=1
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6988)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6989)
@@ -260,6 +260,7 @@
       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          !   -      -
+      REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
       !!---------------------------------------------------------------------
+      zrhoco = rau0 * rn_cio
       !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
@@ -272,6 +273,6 @@
                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
+               taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * zrhoco * zmodt
+               tmod_io(ji,jj) = zrhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point
             END DO
          END DO
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 6988)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 6989)
@@ -18,7 +18,8 @@
    USE sbc_oce        ! ocean surface boundary condition
    USE ice            ! ice variables
-   USE limadv         ! ice advection
    USE limhdf         ! ice horizontal diffusion
    USE limvar         ! 
+   USE limadv_prather ! advection scheme (Prather)
+   USE limadv_umx     ! advection scheme (ultimate-macho)
    !
    USE in_out_manager ! I/O manager
@@ -31,5 +32,4 @@
    USE limcons        ! conservation tests
    USE limctl         ! control prints
-   USE limadv_umx     ! advection scheme
 
    IMPLICIT NONE
@@ -57,5 +57,7 @@
       !! ** method  : variables included in the process are scalar,   
       !!     other values are considered as second order. 
-      !!     For advection, a second order Prather scheme is used.  
+      !!     For advection, one can choose between
+      !!     a) an Ultimate-Macho scheme (whose order is defined by nn_limadv_ord) => nn_limadv=0
+      !!     b) and a second order Prather scheme => nn_limadv=-1
       !!
       !! ** action :
@@ -79,6 +81,5 @@
       ! --- ultimate macho only --- !
       REAL(wp)                               ::   zdt
-      LOGICAL                                ::   lcon
-      REAL(wp), POINTER, DIMENSION(:,:)      ::   ze, zu_trp, zv_trp, z1_v, zudy, zvdx
+      REAL(wp), POINTER, DIMENSION(:,:)      ::   zudy, zvdx, zcu_box, zcv_box
       ! --- prather only --- !
       REAL(wp), POINTER, DIMENSION(:,:)      ::   zarea
@@ -142,5 +143,5 @@
       zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )    ! CFL test for stability
       zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
-      IF(lk_mpp )   CALL mpp_max( zcfl )
+      IF( lk_mpp )   CALL mpp_max( zcfl )
       
       IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
@@ -156,222 +157,208 @@
 !!      ENDIF
 
-#if defined key_limumx
-      !=============================!
-      !==  Ultimate-MACHO scheme  ==!
-      !=============================!
-      
-      CALL wrk_alloc( jpi,jpj, ze, zu_trp, zv_trp, z1_v, zudy, zvdx )
-      
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)''
-         WRITE(numout,*)'lim_adv_umx : Ultimate-MACHO advection scheme'
-         WRITE(numout,*)'~~~~~~~~~~~'
-      ENDIF
-      !
-      zdt = rdt_ice / REAL(initad)
-      
-      ! transport
-      zudy(:,:) = u_ice(:,:) * e2u(:,:)
-      zvdx(:,:) = v_ice(:,:) * e1v(:,:)
-      !
-      DO jt = 1, initad
-         lcon = .TRUE. 
-!!!            lcon = .false. 
-         CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, ato_i(:,:), ato_i(:,:) )                                       ! Open water area 
+      SELECT CASE ( nn_limadv )
+         
+                       !=============================!
+      CASE ( 0 )       !==  Ultimate-MACHO scheme  ==!                   
+                       !=============================!
+      
+         CALL wrk_alloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box )
+      
+         IF( kt == nit000 .AND. lwp ) THEN
+            WRITE(numout,*)''
+            WRITE(numout,*)'lim_adv_umx : Ultimate-MACHO advection scheme'
+            WRITE(numout,*)'~~~~~~~~~~~'
+         ENDIF
          !
+         zdt = rdt_ice / REAL(initad)
+         
+         ! transport
+         zudy(:,:) = u_ice(:,:) * e2u(:,:)
+         zvdx(:,:) = v_ice(:,:) * e1v(:,:)
+         
+         ! define velocity for advection: u*grad(H)
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               IF    ( u_ice(ji,jj) * u_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
+               ELSEIF( u_ice(ji,jj)                  >  0._wp ) THEN   ;   zcu_box(ji,jj) = u_ice(ji-1,jj)
+               ELSE                                                    ;   zcu_box(ji,jj) = u_ice(ji  ,jj)
+               ENDIF
+               
+               IF    ( v_ice(ji,jj) * v_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
+               ELSEIF( v_ice(ji,jj)                  >  0._wp ) THEN   ;   zcv_box(ji,jj) = v_ice(ji,jj-1)
+               ELSE                                                    ;   zcv_box(ji,jj) = v_ice(ji,jj  )
+               ENDIF
+            END DO
+         END DO
+         
+         ! advection
+         DO jt = 1, initad
+            CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, ato_i(:,:) )       ! Open water area 
+            DO jl = 1, jpl
+               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_i(:,:,jl) )      ! Ice area
+               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_i(:,:,jl) )      ! Ice  volume
+               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, smv_i(:,:,jl) )    ! Salt content
+               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, oa_i (:,:,jl) )    ! Age content
+               DO jk = 1, nlay_i
+                  CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_i(:,:,jk,jl) )   ! Ice  heat content
+               END DO
+               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_s(:,:,jl) )      ! Snow volume
+               CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_s(:,:,1,jl) )    ! Snow heat content
+            END DO
+         END DO
+         !
+         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
+         DO jl = 2, jpl
+            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
+         END DO
+         !
+         CALL wrk_dealloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box )
+         
+                       !=============================!
+      CASE ( -1 )      !==     Prather scheme      ==!                   
+                       !=============================!
+
+         CALL wrk_alloc( jpi,jpj,            zarea )
+         CALL wrk_alloc( jpi,jpj,1,          z0opw )
+         CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
+         CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
+         
+         IF( kt == nit000 .AND. lwp ) THEN
+            WRITE(numout,*)''
+            WRITE(numout,*)'lim_adv_xy : Prather advection scheme'
+            WRITE(numout,*)'~~~~~~~~~~~'
+         ENDIF
+         
+         zarea(:,:) = e12t(:,:)
+         
+         !-------------------------
+         ! transported fields                                        
+         !-------------------------
+         z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area 
          DO jl = 1, jpl
-            WHERE( v_i(:,:,jl) /= 0._wp )   ;   z1_v(:,:) = 1._wp / v_i(:,:,jl)
-            ELSEWHERE                       ;   z1_v(:,:) = 0._wp
-            END WHERE
-            !
-            lcon = .TRUE. 
-!!!               lcon = .false. 
-            CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, a_i(:,:,jl), a_i(:,:,jl) )                      ! Ice area
-            CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, v_i(:,:,jl), v_i(:,:,jl), zu_trp, zv_trp )      ! Ice  volume
-            !
-            lcon = .FALSE. 
-            ze(:,:) = smv_i(:,:,jl) * z1_v(:,:)
-            CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, smv_i(:,:,jl) )                           ! Salt content
-            !
-!!!check that               ze(:,:) = oa_i (:,:,jl) * z1_v(:,:)
-!!!check that               CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, oa_i (:,:,jl) )                           ! Age content
-            !
-            zu_trp(:,:) = zu_trp(:,:) * r1_nlay_i
-            zv_trp(:,:) = zv_trp(:,:) * r1_nlay_i
-            z1_v  (:,:) = z1_v  (:,:) * REAL( nlay_i, wp )
+            z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
+            z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
+            z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
+            z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
+            z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
+            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
             DO jk = 1, nlay_i
-               ze (:,:) = e_i(:,:,jk,jl) * z1_v(:,:)
-               CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, e_i(:,:,jk,jl) )                                  ! Ice  heat content
-            END DO
-            !
-            WHERE( v_s(:,:,jl) /= 0._wp )   ;   z1_v(:,:) = 1._wp / v_s(:,:,jl)
-            ELSEWHERE                       ;   z1_v(:,:) = 0._wp
-            END WHERE
-            !
-            lcon = .TRUE. 
-!!!               lcon = .false. 
-            CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zudy, zvdx, v_s(:,:,jl), v_s(:,:,jl), zu_trp, zv_trp )      ! Snow volume
-            !
-            lcon = .FALSE. 
-            ze (:,:) = e_s(:,:,1,jl) * z1_v(:,:)
-            CALL lim_adv_umx( lcon, kt, zdt, zudy, zvdx, zu_trp, zv_trp, ze, e_s(:,:,1,jl) )                                     ! Snow heat content
-            !
-         END DO
-      END DO
-      !
-      at_i(:,:) = a_i(:,:,1)      ! total ice fraction
-      DO jl = 2, jpl
-         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-      END DO
-      !
-      CALL wrk_dealloc( jpi,jpj, ze, zu_trp, zv_trp, z1_v, zudy, zvdx )
-        
-#else
-      !=============================!
-      !==      Prather scheme     ==!
-      !=============================!
-
-      CALL wrk_alloc( jpi,jpj,            zarea )
-      CALL wrk_alloc( jpi,jpj,1,          z0opw )
-      CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
-      CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
-      
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)''
-         WRITE(numout,*)'lim_adv_xy : Prather advection scheme'
-         WRITE(numout,*)'~~~~~~~~~~~'
-      ENDIF
-
-      zarea(:,:) = e12t(:,:)
-      
-      !-------------------------
-      ! transported fields                                        
-      !-------------------------
-      z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area 
-      DO jl = 1, jpl
-         z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
-         z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
-         z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
-         z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
-         z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
-         z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
-         DO jk = 1, nlay_i
-            z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
-         END DO
-      END DO
-
-      
-      IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
-         DO jt = 1, initad
-            CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
-               &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
-               &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            DO jl = 1, jpl
-               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
-                  &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
-                  &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
-                  &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
-                  &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
-                  &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
-                  &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
-                  &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
-                  &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
-                  &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
-                  &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
-                  &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
-                  &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-               DO jk = 1, nlay_i                                                                !--- ice heat contents ---
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                     &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                     &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                     &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                     &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
+            END DO
+         END DO
+         
+         
+         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
+            DO jt = 1, initad
+               CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
+                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
+               CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
+                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
+               DO jl = 1, jpl
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
+                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
+                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
+                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
+                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
+                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
+                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
+                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
+                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
+                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
+                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
+                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
+                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
+                  DO jk = 1, nlay_i                                                                !--- ice heat contents ---
+                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
+                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
+                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+                  END DO
                END DO
             END DO
-         END DO
-      ELSE
-         DO jt = 1, initad
-            CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
-               &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
-               &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            DO jl = 1, jpl
-               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
-                  &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
-                  &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
-                  &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
-                  &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
-                  &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
-                  &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
-                  &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
-                  &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
-                  &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
-                  &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
-                  &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
-                  &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-               DO jk = 1, nlay_i                                                           !--- ice heat contents ---
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                     &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                     &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                     &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                     &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+         ELSE
+            DO jt = 1, initad
+               CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
+                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
+               CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
+                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
+               DO jl = 1, jpl
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
+                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
+                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
+                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
+                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
+                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
+                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
+                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
+                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
+                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
+                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
+                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
+                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
+                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
+                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
+                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
+                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+                  END DO
                END DO
             END DO
-         END DO
-      ENDIF
-      
-      !-------------------------------------------
-      ! Recover the properties from their contents
-      !-------------------------------------------
-      ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
-      DO jl = 1, jpl
-         v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:)
-         v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:)
-         smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:)
-         oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:)
-         a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:)
-         e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
-         DO jk = 1, nlay_i
-            e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
-         END DO
-      END DO
-      
-      at_i(:,:) = a_i(:,:,1)      ! total ice fraction
-      DO jl = 2, jpl
-         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-      END DO
-
-      CALL wrk_dealloc( jpi,jpj,            zarea )
-      CALL wrk_dealloc( jpi,jpj,1,          z0opw )
-      CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
-      CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
-      
-#endif
-
+         ENDIF
+         
+         !-------------------------------------------
+         ! Recover the properties from their contents
+         !-------------------------------------------
+         ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
+         DO jl = 1, jpl
+            v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:)
+            v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:)
+            smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:)
+            oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:)
+            a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:)
+            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
+            DO jk = 1, nlay_i
+               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
+            END DO
+         END DO
+         
+         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
+         DO jl = 2, jpl
+            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
+         END DO
+         
+         CALL wrk_dealloc( jpi,jpj,            zarea )
+         CALL wrk_dealloc( jpi,jpj,1,          z0opw )
+         CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
+         CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
+
+      END SELECT
+      
       !------------------------------!
       ! Diffusion of Ice fields                  
Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6988)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6989)
@@ -321,4 +321,6 @@
          &     'use more ocean levels or less ice/snow layers/categories.' )
       !
+      CALL lim_dyn_init                ! set ice dynamics parameters
+      !
       CALL lim_itd_init                ! ice thickness distribution initialization
       !
