Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90	(revision 10170)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90	(revision 10180)
@@ -71,6 +71,9 @@
       INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
       INTEGER  ::   initad                  ! number of sub-timestep for the advection
+      INTEGER  ::   ipl                     ! third dimention of tracer array
+
       REAL(wp) ::   zcfl , zusnit, zdt      !   -      -
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zudy, zvdx, zcu_box, zcv_box
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zudy, zvdx, zcu_box, zcv_box
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zpato
       !!----------------------------------------------------------------------
       !
@@ -78,4 +81,5 @@
       !
       ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) )
+      ALLOCATE( zpato(jpi,jpj,1) )
       !
       ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !        
@@ -109,34 +113,36 @@
       END DO
 
+      zpato (:,:,1) = pato_i(:,:)
+
       !---------------!
       !== advection ==!
       !---------------!
       DO jt = 1, initad
-         CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area 
-         DO jl = 1, jpl
-            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area
-            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume
-            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,jl) )        ! Salt content
-            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,jl) )        ! Age content
-            DO jk = 1, nlay_i
-               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content
-            END DO
-            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume
-            DO jk = 1, nlay_s
-               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,jl) )   ! Snow heat content
-            END DO
-            IF ( ln_pnd_H12 ) THEN
-               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction
-               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume
-            ENDIF
-         END DO
-      END DO
-      !
-      DEALLOCATE( zudy, zvdx, zcu_box, zcv_box )
+         CALL adv_umx( k_order, kt,   1, zdt, zudy, zvdx, zcu_box, zcv_box, zpato(:,:,1) )        ! Open water area 
+         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,:) )         ! Ice area
+         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,:) )         ! Ice  volume
+         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,:) )        ! Salt content
+         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,:) )        ! Age content
+         DO jk = 1, nlay_i
+            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,:) )   ! Ice  heat content
+         END DO
+         CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,:) )         ! Snow volume
+         DO jk = 1, nlay_s
+            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,:) )   ! Snow heat content
+         END DO
+         IF ( ln_pnd_H12 ) THEN
+            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,:) )     ! Melt pond fraction
+            CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,:) )     ! Melt pond volume
+         ENDIF
+      END DO
+      !
+      pato_i(:,:) = zpato (:,:,1)
+      !
+      DEALLOCATE( zudy, zvdx, zcu_box, zcv_box, zpato )
       !
    END SUBROUTINE ice_dyn_adv_umx
 
    
-   SUBROUTINE adv_umx( k_order, kt, pdt, puc, pvc, pubox, pvbox, ptc )
+   SUBROUTINE adv_umx( k_order, kt, ipl, pdt, puc, pvc, pubox, pvbox, ptc )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE adv_umx  ***
@@ -151,15 +157,19 @@
       !! ** Action : - pt  the after advective tracer
       !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme
-      INTEGER                     , INTENT(in   ) ::   kt             ! number of iteration
-      REAL(wp)                    , INTENT(in   ) ::   pdt            ! tracer time-step
-      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
-      !
-      INTEGER  ::   ji, jj           ! dummy loop indices  
+      INTEGER                         , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme
+      INTEGER                         , INTENT(in   ) ::   kt             ! number of iteration
+      INTEGER                         , INTENT(in   ) ::   ipl            ! third dimension of tracer array
+      REAL(wp)                        , INTENT(in   ) ::   pdt            ! tracer time-step
+      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,ipl), INTENT(inout) ::   ptc            ! tracer content field
+      !
+      INTEGER  ::   ji, jj, jl       ! dummy loop indices  
+
       REAL(wp) ::   ztra             ! local scalar
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfu_ho, zt_u, zt_ups
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfv_ups, zfv_ho, zt_v, ztrd
+      REAL(wp), DIMENSION(jpi,jpj,ipl) ::   zfu_ups, zfu_ho, zt_u, zt_ups
+      REAL(wp), DIMENSION(jpi,jpj,ipl) ::   zfv_ups, zfv_ho, zt_v, ztrd
+
+      DO jl = 1, ipl
       !!----------------------------------------------------------------------
       !
@@ -168,19 +178,19 @@
       DO jj = 1, jpjm1         ! upstream tracer flux in the i and j direction
          DO ji = 1, fs_jpim1   ! vector opt.
-            zfu_ups(ji,jj) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj)
-            zfv_ups(ji,jj) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1)
+            zfu_ups(ji,jj,jl) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj,jl)
+            zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1,jl)
          END DO
       END DO
       
-      DO jj = 2, jpjm1            ! total intermediate advective trends
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            ztra = - (   zfu_ups(ji,jj) - zfu_ups(ji-1,jj  )   &
-               &       + zfv_ups(ji,jj) - zfv_ups(ji  ,jj-1)   ) * r1_e1e2t(ji,jj)
+         DO jj = 2, jpjm1            ! total intermediate advective trends
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ztra = - (   zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj  ,jl)   &
+                  &       + zfv_ups(ji,jj,jl) - zfv_ups(ji  ,jj-1,jl)   ) * r1_e1e2t(ji,jj)
             !
-            ztrd(ji,jj) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ]  
-            zt_ups (ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)      ! guess after content field with monotonic scheme
-         END DO
-      END DO
-      CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
+               ztrd(ji,jj,jl) =                         ztra                         ! upstream trend [ -div(uh) or -div(uhT) ]  
+               zt_ups (ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1)   ! guess after content field with monotonic scheme
+            END DO
+         END DO
+      END DO
       
       ! High order (_ho) fluxes 
@@ -188,18 +198,22 @@
       SELECT CASE( k_order )
       CASE ( 20 )                          ! centered second order
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               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) )
+         DO jl = 1, ipl
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji+1,jj,jl) )
+                  zfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji,jj+1,jl) )
+               END DO
             END DO
          END DO
          !
       CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme
-         CALL macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v )
-         !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zfu_ho(ji,jj) = puc(ji,jj) * zt_u(ji,jj)
-               zfv_ho(ji,jj) = pvc(ji,jj) * zt_v(ji,jj)
+         CALL macho( k_order, kt, ipl, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v )
+         !
+         DO jl = 1, ipl
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zfu_ho(ji,jj,jl) = puc(ji,jj) * zt_u(ji,jj,jl)
+                  zfv_ho(ji,jj,jl) = pvc(ji,jj) * zt_v(ji,jj,jl)
+               END DO
             END DO
          END DO
@@ -209,22 +223,28 @@
       ! antidiffusive flux : high order minus low order
       ! --------------------------------------------------
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector opt.
-            zfu_ho(ji,jj) = zfu_ho(ji,jj) - zfu_ups(ji,jj)
-            zfv_ho(ji,jj) = zfv_ho(ji,jj) - zfv_ups(ji,jj)
-         END DO
-      END DO
+      DO jl = 1, ipl
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zfu_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) - zfu_ups(ji,jj,jl)
+               zfv_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) - zfv_ups(ji,jj,jl)
+            END DO
+         END DO
+      END DO
+
+      CALL lbc_lnk("icedyn_adv_umx",zt_ups, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
       
       ! monotonicity algorithm
       ! -------------------------
-      CALL nonosc_2d( ptc, zfu_ho, zfv_ho, zt_ups, pdt )
+      CALL nonosc_2d( ipl, ptc, zfu_ho, zfv_ho, zt_ups, pdt )
       
       ! final trend with corrected fluxes
       ! ------------------------------------
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.  
-            ztra       = ztrd(ji,jj)  - (  zfu_ho(ji,jj) - zfu_ho(ji-1,jj  )   &
-               &                         + zfv_ho(ji,jj) - zfv_ho(ji  ,jj-1) ) * r1_e1e2t(ji,jj)  
-            ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)
+      DO jl = 1, ipl
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.  
+               ztra       = ztrd(ji,jj,jl)  - (  zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj  ,jl)   &
+                  &                         +    zfv_ho(ji,jj,jl) - zfv_ho(ji  ,jj-1,jl) ) * r1_e1e2t(ji,jj)  
+               ptc(ji,jj,jl) = ptc(ji,jj,jl) + pdt * ztra * tmask(ji,jj,1)
+            END DO
          END DO
       END DO
@@ -233,6 +253,5 @@
    END SUBROUTINE adv_umx
 
-
-   SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v )
+   SUBROUTINE macho( k_order, kt, ipl, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_x  ***
@@ -245,15 +264,16 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme
-      INTEGER                     , INTENT(in   ) ::   kt         ! number of iteration
-      REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
-      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 
-      !
-      INTEGER  ::   ji, jj    ! dummy loop indices
+      INTEGER                         , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme
+      INTEGER                         , INTENT(in   ) ::   kt         ! number of iteration
+      INTEGER                         , INTENT(in   ) ::   ipl        ! third dimension of tracer array
+      REAL(wp)                        , INTENT(in   ) ::   pdt        ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,ipl), 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,ipl), INTENT(  out) ::   pt_u, pt_v ! tracer at u- and v-points 
+      !
+      INTEGER  ::   ji, jj, jl    ! dummy loop indices
       REAL(wp) ::   zc_box    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: zzt
+      REAL(wp), DIMENSION(jpi,jpj,ipl) :: zzt
       !!----------------------------------------------------------------------
       !
@@ -261,12 +281,14 @@
          !
          !                                                           !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( k_order, pdt, ptc, puc, pt_u )
+         CALL ultimate_x( k_order, ipl, pdt, ptc, puc, pt_u )
          !
          !                                                           !--  advective form update in zzt  --!
-         DO jj = 2, jpjm1
-            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_e1e2t(ji,jj)
-               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
+         DO jl = 1, ipl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zzt(ji,jj,jl) = ptc(ji,jj,jl) - pubox(ji,jj) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj)  &
+                   &            - ptc(ji,jj,jl) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj)
+                  zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1)
+               END DO
             END DO
          END DO
@@ -274,17 +296,19 @@
          !
          !                                                           !--  ultimate interpolation of pt at v-point  --!
-         CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v )
+         CALL ultimate_y( k_order, ipl, 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 )
+         CALL ultimate_y( k_order, ipl, 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_e1e2t(ji,jj)
-               zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1)
+         DO jl = 1, ipl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zzt(ji,jj,jl) = ptc(ji,jj,jl) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj)  &
+                     &                    - ptc  (ji,jj,jl) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj)
+                  zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1)
+               END DO
             END DO
          END DO
@@ -292,5 +316,5 @@
          !
          !                                                           !--  ultimate interpolation of pt at u-point  --!
-         CALL ultimate_x( k_order, pdt, zzt, puc, pt_u )
+         CALL ultimate_x( k_order, ipl, pdt, zzt, puc, pt_u )
          !      
       ENDIF      
@@ -299,5 +323,5 @@
 
 
-   SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u )
+   SUBROUTINE ultimate_x( k_order, ipl, pdt, pt, puc, pt_u )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_x  ***
@@ -310,23 +334,27 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index
-      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
-      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
+      INTEGER                         , INTENT(in   ) ::   k_order   ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   ipl        ! third dimension of tracer array
+      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   puc       ! ice i-velocity component
+      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in   ) ::   pt        ! tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_u      ! tracer at u-point 
+      !
+      INTEGER  ::   ji, jj, jl       ! dummy loop indices
       REAL(wp) ::   zcu, zdx2, zdx4    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4
+      REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu3
+      REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztu2, 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)
+      DO jl = 1, ipl
+         DO jj = 2, jpjm1         ! First derivative (gradient)
+            DO ji = 1, fs_jpim1
+               ztu1(ji,jj) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
+            END DO
+            !                     ! Second derivative (Laplacian)
+            DO ji = fs_2, fs_jpim1
+               ztu2(ji,jj,jl) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj)
+            END DO
          END DO
       END DO
@@ -334,11 +362,13 @@
       !
       !                                                     !--  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)
+      DO jl = 1, ipl
+         DO jj = 2, jpjm1         ! Third derivative
+            DO ji = 1, fs_jpim1
+               ztu3(ji,jj) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1)
+            END DO
+            !                     ! Fourth derivative
+            DO ji = fs_2, fs_jpim1
+               ztu4(ji,jj,jl) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj)
+            END DO
          END DO
       END DO
@@ -350,8 +380,10 @@
       CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21)
          !        
-         DO jj = 2, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               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) ) )
+         DO jl = 1, ipl
+            DO jj = 2, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
+                     &                                       - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
+               END DO
             END DO
          END DO
@@ -359,9 +391,11 @@
       CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23)
          !
-         DO jj = 2, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               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) ) ) 
+         DO jl = 1, ipl
+            DO jj = 2, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zcu  = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
+                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                               pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
+                     &                                              -              zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
+               END DO
             END DO
          END DO
@@ -369,13 +403,15 @@
       CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24)
          !
-         DO jj = 2, jpjm1
-            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)
+         DO jl = 1, ipl
+            DO jj = 2, jpjm1
+               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)
-               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)        &
-                  &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) )  )   )
+                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                     pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        &
+                     &                                              -              zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   &
+                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                        ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        &
+                     &                                              - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   )
+               END DO
             END DO
          END DO
@@ -383,13 +419,15 @@
       CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27)
          !
-         DO jj = 2, jpjm1
-            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)
+         DO jl = 1, ipl
+            DO jj = 2, jpjm1
+               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)
-               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) )  )   )
+                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                     pt  (ji+1,jj,jl) + pt  (ji,jj,jl)        &
+                     &                                                    -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) )  )   &
+                     &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * (                        ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)        &
+                     &                                                    - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) )  )   )
+               END DO
             END DO
          END DO
@@ -397,16 +435,18 @@
       CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29)
          !
-         DO jj = 2, jpjm1
-            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)
+         DO jl = 1, ipl
+            DO jj = 2, jpjm1
+               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) ) ) )
+                  zdx4 = zdx2 * zdx2
+                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (           (                   pt  (ji+1,jj,jl) + pt  (ji,jj,jl)       &
+                     &                                                    -          zcu * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) )   &
+                     &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) *    (                   ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)       &
+                     &                                                    - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) )   &
+                     &       + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)       &
+                     &                                              - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
+               END DO
             END DO
          END DO
@@ -417,5 +457,5 @@
    
  
-   SUBROUTINE ultimate_y( k_order, pdt, pt, pvc, pt_v )
+   SUBROUTINE ultimate_y( k_order, ipl, pdt, pt, pvc, pt_v )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_y  ***
@@ -428,24 +468,28 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in   ) ::   k_order   ! ocean time-step index
-      REAL(wp)                    , INTENT(in   ) ::   pdt       ! tracer time-step
-      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
+      INTEGER                         , INTENT(in   ) ::   k_order   ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   ipl       ! third dimension of tracer array
+      REAL(wp)                        , INTENT(in   ) ::   pdt       ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pvc       ! ice j-velocity component
+      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in   ) ::   pt        ! tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(  out) ::   pt_v      ! tracer at v-point 
+      !
+      INTEGER  ::   ji, jj, jl       ! dummy loop indices
       REAL(wp) ::   zcv, zdy2, zdy4    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4
+      REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv3
+      REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztv2, 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)
+      DO jl = 1, ipl
+         DO jj = 1, jpjm1         ! First derivative (gradient)
+            DO ji = fs_2, fs_jpim1
+               ztv1(ji,jj) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * 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,jl) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj)
+            END DO
          END DO
       END DO
@@ -453,12 +497,14 @@
       !
       !                                                     !--  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)
+      DO jl = 1, ipl
+         DO jj = 1, jpjm1         ! First derivative
+            DO ji = fs_2, fs_jpim1
+            ztv3(ji,jj) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * 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,jl) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj)
+            END DO
          END DO
       END DO
@@ -469,17 +515,21 @@
       !
       CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
-         DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               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) ) )
+         DO jl = 1, ipl
+            DO jj = 1, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                           ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
+                     &                                     - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
+               END DO
             END DO
          END DO
          !
       CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
-         DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               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) ) )
+         DO jl = 1, ipl
+            DO jj = 1, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
+                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (     ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
+                     &                                     - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
+               END DO
             END DO
          END DO
@@ -487,42 +537,48 @@
          !
       CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
-         DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
-               zdy2 = e2v(ji,jj) * e2v(ji,jj)
+         DO jl = 1, ipl
+            DO jj = 1, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  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)
-               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) ) ) )
+                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)       &
+                     &                                     -                        zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )   &
+                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)       &
+                     &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
+               END DO
             END DO
          END DO
          !
       CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
-         DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
-               zdy2 = e2v(ji,jj) * e2v(ji,jj)
+         DO jl = 1, ipl
+            DO jj = 1, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  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)
-               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) ) ) )
+                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                        ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)       &
+                     &                                               -          zcv * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )   &
+                     &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                   ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)       &
+                     &                                               - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
+               END DO
             END DO
          END DO
          !
       CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
-         DO jj = 1, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zcv  = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
-               zdy2 = e2v(ji,jj) * e2v(ji,jj)
+         DO jl = 1, ipl
+            DO jj = 1, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  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) ) ) )
+                  zdy4 = zdy2 * zdy2
+                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)      &
+                     &                                                     -          zcv * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) )  &
+                     &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) *     (                   ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)      &
+                     &                                                     - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) )  &
+                     &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)      &
+                     &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
+               END DO
             END DO
          END DO
@@ -533,5 +589,5 @@
    
   
-   SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, pdt )
+   SUBROUTINE nonosc_2d( ipl, pbef, paa, pbb, paft, pdt )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nonosc  ***
@@ -546,13 +602,15 @@
       !!       in-space based differencing for fluid
       !!----------------------------------------------------------------------
-      REAL(wp)                     , INTENT(in   ) ::   pdt          ! tracer time-step
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(in   ) ::   pbef, paft   ! before & after field
-      REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) ::   paa, pbb     ! monotonic fluxes in the 2 directions
-      !
-      INTEGER  ::   ji, jj    ! dummy loop indices
+      INTEGER                          , INTENT(in   ) ::   ipl          ! third dimension of tracer array
+      REAL(wp)                         , INTENT(in   ) ::   pdt          ! tracer time-step
+      REAL(wp), DIMENSION (jpi,jpj,ipl), INTENT(in   ) ::   pbef, paft   ! before & after field
+      REAL(wp), DIMENSION (jpi,jpj,ipl), INTENT(inout) ::   paa, pbb     ! monotonic fluxes in the 2 directions
+      !
+      INTEGER  ::   ji, jj, jl    ! dummy loop indices
       INTEGER  ::   ikm1      ! local integer
       REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt   ! local scalars
       REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zmsk, zdiv
+      REAL(wp), DIMENSION(jpi,jpj) :: zbup, zbdo, zmsk
+      REAL(wp), DIMENSION(jpi,jpj,ipl) :: zbetup, zbetdo, zdiv
       !!----------------------------------------------------------------------
       !
@@ -561,28 +619,33 @@
 
       ! test on divergence
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.  
-            zdiv(ji,jj) =  - (  paa(ji,jj) - paa(ji-1,jj  )   &
-               &              + pbb(ji,jj) - pbb(ji  ,jj-1) )  
+      DO jl = 1, ipl
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.  
+               zdiv(ji,jj,jl) =  - (  paa(ji,jj,jl) - paa(ji-1,jj  ,jl)   &
+                  &              +    pbb(ji,jj,jl) - pbb(ji  ,jj-1,jl) )  
+            END DO
          END DO
       END DO
       CALL lbc_lnk( 'icedyn_adv_umx', zdiv, 'T', 1. )        ! Lateral boundary conditions   (unchanged sign)
 
-      ! Determine ice masks for before and after tracers 
-      WHERE( pbef(:,:) == 0._wp .AND. paft(:,:) == 0._wp .AND. zdiv(:,:) == 0._wp )   ;   zmsk(:,:) = 0._wp
-      ELSEWHERE                                                                       ;   zmsk(:,:) = 1._wp * tmask(:,:,1)
-      END WHERE
-
-      ! Search local extrema
-      ! --------------------
-      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
-!      zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ),   &
-!         &             paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) )  )
-!      zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ),   &
-!         &             paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) )  )
-      zbup(:,:) = MAX( pbef(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   &
-         &             paft(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  )
-      zbdo(:,:) = MIN( pbef(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   &
-         &             paft(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  )
+      DO jl = 1, ipl
+         ! Determine ice masks for before and after tracers 
+         WHERE( pbef(:,:,jl) == 0._wp .AND. paft(:,:,jl) == 0._wp .AND. zdiv(:,:,jl) == 0._wp ) 
+            zmsk(:,:) = 0._wp
+         ELSEWHERE                                                                                   
+            zmsk(:,:) = 1._wp * tmask(:,:,1)
+         END WHERE
+
+         ! Search local extrema
+         ! --------------------
+         ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
+!         zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ),   &
+!            &             paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) )  )
+!         zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ),   &
+!            &             paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) )  )
+         zbup(:,:) = MAX( pbef(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ),   &
+            &             paft(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) )  )
+         zbdo(:,:) = MIN( pbef(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ),   &
+            &             paft(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) )  )
 
       z1_dt = 1._wp / pdt
@@ -595,35 +658,38 @@
                &                       zbdo(ji  ,jj-1), zbdo(ji  ,jj+1)    )
                !
-            zpos = MAX( 0., paa(ji-1,jj  ) ) - MIN( 0., paa(ji  ,jj  ) )   &        ! positive/negative  part of the flux
-               & + MAX( 0., pbb(ji  ,jj-1) ) - MIN( 0., pbb(ji  ,jj  ) )
-            zneg = MAX( 0., paa(ji  ,jj  ) ) - MIN( 0., paa(ji-1,jj  ) )   &
-               & + MAX( 0., pbb(ji  ,jj  ) ) - MIN( 0., pbb(ji  ,jj-1) )
+               zpos = MAX( 0., paa(ji-1,jj  ,jl) ) - MIN( 0., paa(ji  ,jj  ,jl) )   &        ! positive/negative  part of the flux
+                  & + MAX( 0., pbb(ji  ,jj-1,jl) ) - MIN( 0., pbb(ji  ,jj  ,jl) )
+               zneg = MAX( 0., paa(ji  ,jj  ,jl) ) - MIN( 0., paa(ji-1,jj  ,jl) )   &
+                  & + MAX( 0., pbb(ji  ,jj  ,jl) ) - MIN( 0., pbb(ji  ,jj-1,jl) )
                !
-            zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms
-            zbetup(ji,jj) = ( zup         - paft(ji,jj) ) / ( zpos + zsml ) * zbt
-            zbetdo(ji,jj) = ( paft(ji,jj) - zdo         ) / ( zneg + zsml ) * zbt
+               zbt = e1e2t(ji,jj) * z1_dt                                   ! up & down beta terms
+               zbetup(ji,jj,jl) = ( zup            - paft(ji,jj,jl) ) / ( zpos + zsml ) * zbt
+               zbetdo(ji,jj,jl) = ( paft(ji,jj,jl) - zdo         )    / ( zneg + zsml ) * zbt
+            END DO
          END DO
       END DO
       CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
 
-      ! monotonic flux in the i & j direction (paa & pbb)
-      ! -------------------------------------
-      DO jj = 2, jpjm1
-         DO ji = 1, fs_jpim1   ! vector opt.
-            zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) )
-            zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) )
-            zcu = 0.5  + SIGN( 0.5 , paa(ji,jj) )
-            !
-            paa(ji,jj) = paa(ji,jj) * ( zcu * zau + ( 1._wp - zcu) * zbu )
-         END DO
-      END DO
-      !
-      DO jj = 1, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) )
-            zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) )
-            zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj) )
-            !
-            pbb(ji,jj) = pbb(ji,jj) * ( zcv * zav + ( 1._wp - zcv) * zbv )
+      DO jl = 1, ipl
+         ! monotonic flux in the i & j direction (paa & pbb)
+         ! -------------------------------------
+         DO jj = 2, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
+               zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
+               zcu = 0.5  + SIGN( 0.5 , paa(ji,jj,jl) )
+               !
+               paa(ji,jj,jl) = paa(ji,jj,jl) * ( zcu * zau + ( 1._wp - zcu) * zbu )
+            END DO
+         END DO
+         !
+         DO jj = 1, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
+               zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
+               zcv = 0.5  + SIGN( 0.5 , pbb(ji,jj,jl) )
+               !
+               pbb(ji,jj,jl) = pbb(ji,jj,jl) * ( zcv * zav + ( 1._wp - zcv) * zbv )
+            END DO
          END DO
       END DO
