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 6860)
+++ 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')
