Index: NEMO/trunk/src/ICE/icectl.F90
===================================================================
--- NEMO/trunk/src/ICE/icectl.F90	(revision 11601)
+++ NEMO/trunk/src/ICE/icectl.F90	(revision 11612)
@@ -158,7 +158,7 @@
             !    only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc)
             !    so the formulation for conservation is different (and not coded) 
-            !    it does not mean UM is not conservative (it is checked with above prints)
-            IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) &
-               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice
+            !    it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now
+            !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) &
+            !   &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice
          ENDIF
          !
Index: NEMO/trunk/src/ICE/icedyn_adv_pra.F90
===================================================================
--- NEMO/trunk/src/ICE/icedyn_adv_pra.F90	(revision 11601)
+++ NEMO/trunk/src/ICE/icedyn_adv_pra.F90	(revision 11612)
@@ -19,4 +19,5 @@
    USE ice            ! sea-ice variables
    USE sbc_oce , ONLY : nn_fsbc   ! frequency of sea-ice call
+   USE icevar         ! sea-ice: operations
    !
    USE in_out_manager ! I/O manager
@@ -25,5 +26,4 @@
    USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
    USE lbclnk         ! lateral boundary conditions (or mpp links)
-   USE prtctl         ! Print control
 
    IMPLICIT NONE
@@ -39,5 +39,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   ! ice salinity
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   ! ice age
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    ! snow layers heat content
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content
@@ -82,30 +82,31 @@
       !
       INTEGER  ::   jk, jl, jt              ! dummy loop indices
-      INTEGER  ::   initad                  ! number of sub-timestep for the advection
-      REAL(wp) ::   zcfl , zusnit           !   -      -
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zarea
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0opw
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ice, z0snw, z0ai, z0smi, z0oi
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0es
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei
+      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
+      REAL(wp) ::   zdt                     !   -      -
+      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication
+      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea
+      REAL(wp), DIMENSION(jpi,jpj,1)          ::   z0opw
+      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi
+      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp
+      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es
+      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei
       !!----------------------------------------------------------------------
       !
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
       !
-      ALLOCATE( zarea(jpi,jpj)    , z0opw(jpi,jpj, 1 ) , z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) ,                       &
-         &      z0ai(jpi,jpj,jpl) , z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap (jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) ,   &
-         &      z0es (jpi,jpj,nlay_s,jpl), z0ei(jpi,jpj,nlay_i,jpl) )
-      !
-      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !        
-      zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
-      zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
-      CALL mpp_max( 'icedyn_adv_pra', zcfl )
+      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- !
+      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
+      !              this should not affect too much the stability
+      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
+      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
       
-      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
-      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
+      ! non-blocking global communication send zcflnow and receive zcflprv
+      CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
+
+      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2
+      ELSE                         ;   icycle = 1
       ENDIF
+      zdt = rdt_ice / REAL(icycle)
       
-      zarea(:,:) = e1e2t(:,:)
       !-------------------------
       ! transported fields                                        
@@ -113,9 +114,10 @@
       z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:)              ! Open water area 
       DO jl = 1, jpl
-         z0snw(:,:,jl) = pv_s (:,:,  jl) * e1e2t(:,:)     ! Snow volume
-         z0ice(:,:,jl) = pv_i (:,:,  jl) * e1e2t(:,:)     ! Ice  volume
-         z0ai (:,:,jl) = pa_i (:,:,  jl) * e1e2t(:,:)     ! Ice area
-         z0smi(:,:,jl) = psv_i(:,:,  jl) * e1e2t(:,:)     ! Salt content
-         z0oi (:,:,jl) = poa_i(:,:,  jl) * e1e2t(:,:)     ! Age content
+         zarea(:,:,jl) = e1e2t(:,:)
+         z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:)        ! Snow volume
+         z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:)        ! Ice  volume
+         z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:)        ! Ice area
+         z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:)        ! Salt content
+         z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:)        ! Age content
          DO jk = 1, nlay_s
             z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
@@ -133,106 +135,72 @@
       IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
          !                                                 !--------------------------------------------!
-         DO jt = 1, initad
-            CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
-               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
-               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            DO jl = 1, jpl
-               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
-                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
-                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
-                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
-                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
-                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
-                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
-                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
-                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
-                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
-                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               DO jk = 1, nlay_s                                                               !--- snow heat contents ---
-                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
-                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
-                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
-                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
-               END DO
-               DO jk = 1, nlay_i                                                               !--- ice heat contents ---
-                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-               END DO
-               IF ( ln_pnd_H12 ) THEN
-                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction --
-                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
-                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   --
-                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
-                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-               ENDIF
-            END DO
+         DO jt = 1, icycle
+            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) !--- open water
+            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw )
+            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
+            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
+            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
+            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
+            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
+            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
+            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
+            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
+            CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
+            CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
+            !
+            DO jk = 1, nlay_s                                                                             !--- snow heat content
+               CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+               CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+            END DO
+            DO jk = 1, nlay_i                                                                             !--- ice heat content
+               CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+               CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+            END DO
+            !
+            IF ( ln_pnd_H12 ) THEN
+               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
+               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
+               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
+               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
+            ENDIF
          END DO
       !                                                    !--------------------------------------------!
       ELSE                                                 !== even ice time step:  adv_y then adv_x  ==!
          !                                                 !--------------------------------------------!
-         DO jt = 1, initad
-            CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
-               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
-               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            DO jl = 1, jpl
-               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
-                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
-                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
-                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
-                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
-                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
-                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
-                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
-                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
-                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
-                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               DO jk = 1, nlay_s                                                             !--- snow heat contents ---
-                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
-                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
-                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
-                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
-               END DO
-               DO jk = 1, nlay_i                                                             !--- ice heat contents ---
-                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-               END DO
-               IF ( ln_pnd_H12 ) THEN
-                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction ---
-                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &
-                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   ---
-                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &
-                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-               ENDIF
-            END DO
+         DO jt = 1, icycle
+            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) !--- open water
+            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw )
+            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
+            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
+            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
+            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
+            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
+            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
+            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
+            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
+            CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
+            CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
+            DO jk = 1, nlay_s                                                                             !--- snow heat content
+               CALL adv_y( zdt, pv_ice, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+               CALL adv_x( zdt, pu_ice, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+            END DO
+            DO jk = 1, nlay_i                                                                             !--- ice heat content
+               CALL adv_y( zdt, pv_ice, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+               CALL adv_x( zdt, pu_ice, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+            END DO
+            IF ( ln_pnd_H12 ) THEN
+               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
+               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
+               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
+               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
+            ENDIF
          END DO
       ENDIF
@@ -243,9 +211,9 @@
       pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) * tmask(:,:,1)
       DO jl = 1, jpl
-         pv_i (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         pv_s (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         psv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         poa_i(:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         pa_i (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+         pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+         pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+         psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+         poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+         pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
          DO jk = 1, nlay_s
             pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
@@ -255,10 +223,13 @@
          END DO
          IF ( ln_pnd_H12 ) THEN
-            pa_ip  (:,:,jl) = z0ap (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-            pv_ip  (:,:,jl) = z0vp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
          ENDIF
       END DO
       !
-      DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0smi , z0oi , z0ap , z0vp , z0es, z0ei )
+      ! --- Ensure non-negative fields --- !
+      ! Remove negative values (conservation is ensured)
+      !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
+      CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
       !
       IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file
@@ -267,5 +238,5 @@
    
    
-   SUBROUTINE adv_x( pdf, put , pcrh, psm , ps0 ,   &
+   SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 ,   &
       &              psx, psxx, psy , psyy, psxy )
       !!----------------------------------------------------------------------
@@ -275,14 +246,14 @@
       !!                variable on x axis
       !!----------------------------------------------------------------------
-      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
-      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then 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
+      REAL(wp)                  , INTENT(in   ) ::   pdt                ! the time step
+      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
+      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments 
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
       !! 
-      INTEGER  ::   ji, jj                               ! dummy loop indices
-      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars
+      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
+      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars
       REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
       REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
@@ -291,154 +262,150 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   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
+      !
+      jcat = SIZE( ps0 , 3 )   ! size of input arrays
+      !
+      DO jl = 1, jcat   ! loop on categories
+         !
+         ! Limitation of moments.                                           
+         DO jj = 2, jpjm1
+            DO ji = 1, jpi
+               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
+               psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 )
+               !
+               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
+               zs1max  = 1.5 * zslpmax
+               zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) )
+               zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
+                  &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  )
+               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
+
+               ps0 (ji,jj,jl) = zslpmax  
+               psx (ji,jj,jl) = zs1new         * rswitch
+               psxx(ji,jj,jl) = zs2new         * rswitch
+               psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch
+               psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch
+               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
+            END DO
+         END DO
+
+         !  Calculate fluxes and moments between boxes i<-->i+1              
+         DO jj = 2, jpjm1                      !  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) ) * pdt * e2u(ji,jj) / psm(ji,jj,jl)
+               zalfq        =  zalf * zalf
+               zalf1        =  1.0 - zalf
+               zalf1q       =  zalf1 * zalf1
+               !
+               zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl)
+               zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) )
+               zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) )
+               zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq
+               zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
+               zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl)
+               zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl)
+
+               !  Readjust moments remaining in the box.
+               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
+               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
+               psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) )
+               psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl)
+               psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj)
+               psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj)
+               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
+            END DO
+         END DO
+
+         DO jj = 2, jpjm1                      !  Flux from i+1 to i when u LT 0.
+            DO ji = 1, fs_jpim1
+               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji+1,jj,jl) 
+               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,jl)
+               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) &
+                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
+               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
+               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq
+               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
+               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl)
+               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl)
+            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,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) )
+               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) )
+               psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) )
+               psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl)
+               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) )
+               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) )
+               psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl)
+            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,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl)
+               zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl)
+               zalf1         = 1.0 - zalf
+               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj)
+               !
+               ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl)
+               psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl)
+               psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             &
+                  &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) &
+                  &            + zbt1 * psxx(ji,jj,jl)
+               psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             &
+                  &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   &
+                  &            + zbt1 * psxy(ji,jj,jl)
+               psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl)
+               psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl)
+            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,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
+               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
+               zalf1         = 1.0 - zalf
+               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
+               !
+               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) )
+               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp )
+               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) &
+                  &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    &
+                  &                                           + ( zalf1 - zalf ) * ztemp ) )
+               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
+                  &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) )
+               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) )
+               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) )
+            END DO
+         END DO
+
       END DO
 
-      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
-      psm (:,:)  = MAX( pcrh * e1e2t(:,:) + ( 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( 'icedyn_adv_pra', 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=' adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
-         CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
-         CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
-         CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_x: psxy :')
-      ENDIF
+      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , '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. )
       !
    END SUBROUTINE adv_x
 
 
-   SUBROUTINE adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
+   SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 ,   &
       &              psx, psxx, psy , psyy, psxy )
       !!---------------------------------------------------------------------
@@ -448,14 +415,14 @@
       !!                variable on y axis
       !!---------------------------------------------------------------------
-      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
-      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then 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
+      REAL(wp)                  , INTENT(in   ) ::   pdt                ! time step
+      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
+      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments 
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
       !!
-      INTEGER  ::   ji, jj                               ! dummy loop indices
-      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars
+      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
+      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars
       REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
       REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
@@ -464,150 +431,147 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   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
+      !
+      jcat = SIZE( ps0 , 3 )   ! size of input arrays
+      !      
+      DO jl = 1, jcat   ! loop on categories
+         !
+         ! Limitation of moments.
+         DO jj = 1, jpj
+            DO ji = fs_2, fs_jpim1
+               !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
+               psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  )
+               !
+               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
+               zs1max  = 1.5 * zslpmax
+               zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) )
+               zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
+                  &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  )
+               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
+               !
+               ps0 (ji,jj,jl) = zslpmax  
+               psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch
+               psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch
+               psy (ji,jj,jl) = zs1new         * rswitch
+               psyy(ji,jj,jl) = zs2new         * rswitch
+               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
+            END DO
+         END DO
+ 
+         !  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 = fs_2, fs_jpim1
+               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
+               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt * e1v(ji,jj) / psm(ji,jj,jl)
+               zalfq        =  zalf * zalf
+               zalf1        =  1.0 - zalf
+               zalf1q       =  zalf1 * zalf1
+               !
+               zfm (ji,jj)  =  zalf  * psm(ji,jj,jl)
+               zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 
+               zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) )
+               zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl)
+               zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
+               zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl)
+               zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl)
+               !
+               !  Readjust moments remaining in the box.
+               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
+               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
+               psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) )
+               psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl)
+               psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj)
+               psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj)
+               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
+            END DO
+         END DO
+         !
+         DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
+            DO ji = fs_2, fs_jpim1
+               zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * pdt * e1v(ji,jj) ) / psm(ji,jj+1,jl) 
+               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,jl)
+               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji,jj+1,jl) &
+                  &                                   - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
+               zfy   (ji,jj) = zfy (ji,jj) + zalfq * (  psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
+               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji,jj+1,jl) * zalfq
+               zfx   (ji,jj) = zfx (ji,jj) + zalf  * (  psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
+               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji,jj+1,jl)
+               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji,jj+1,jl)
+            END DO
+         END DO
+
+         !  Readjust moments remaining in the box. 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zbt  =         zbet(ji,jj-1)
+               zbt1 = ( 1.0 - zbet(ji,jj-1) )
+               !
+               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) )
+               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) )
+               psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) )
+               psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl)
+               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) )
+               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) )
+               psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl)
+            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 = fs_2, fs_jpim1
+               zbt  =       zbet(ji,jj-1)
+               zbt1 = 1.0 - zbet(ji,jj-1)
+               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 
+               zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 
+               zalf1         = 1.0 - zalf
+               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1)
+               !
+               ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl)
+               psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  &
+                  &             + zbt1 * psy(ji,jj,jl)  
+               psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           &
+                  &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 
+                  &             + zbt1 * psyy(ji,jj,jl)
+               psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            &
+                  &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  &
+                  &             + zbt1 * psxy(ji,jj,jl)
+               psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl)
+               psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl)
+            END DO
+         END DO
+
+         DO jj = 2, jpjm1                      !  Flux from j+1 to j IF v LT 0.
+            DO ji = fs_2, fs_jpim1
+               zbt  =       zbet(ji,jj)
+               zbt1 = 1.0 - zbet(ji,jj)
+               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
+               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
+               zalf1         = 1.0 - zalf
+               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
+               !
+               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) )
+               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )
+               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) &
+                  &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    &
+                  &                                            + ( zalf1 - zalf ) * ztemp ) )
+               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
+                  &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) )
+               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) )
+               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) )
+            END DO
+         END DO
+
       END DO
 
-      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
-      psm(:,:)  = MAX(  pcrh * e1e2t(:,:) + ( 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( 'icedyn_adv_pra', 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=' adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
-         CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
-         CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
-         CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_y: psxy :')
-      ENDIF
+      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , '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. )
       !
    END SUBROUTINE adv_y
@@ -624,5 +588,5 @@
       !
       !                             !* allocate prather fields
-      ALLOCATE( sxopw(jpi,jpj)     , syopw(jpi,jpj)     , sxxopw(jpi,jpj)     , syyopw(jpi,jpj)     , sxyopw(jpi,jpj)     ,   &
+      ALLOCATE( sxopw(jpi,jpj,1)   , syopw(jpi,jpj,1)   , sxxopw(jpi,jpj,1)   , syyopw(jpi,jpj,1)   , sxyopw(jpi,jpj,1)   ,   &
          &      sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
          &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   &
@@ -652,5 +616,5 @@
       !!                   ***  ROUTINE adv_pra_rst  ***
       !!                     
-      !! ** Purpose :   Read or write RHG file in restart file
+      !! ** Purpose :   Read or write file in restart file
       !!
       !! ** Method  :   use of IOM library
Index: NEMO/trunk/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 11601)
+++ NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 11612)
@@ -319,9 +319,7 @@
          !
          !== Ice age ==!
-         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN
-            zamsk = 1._wp
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
-               &                                      poa_i, poa_i )
-         ENDIF
+         zamsk = 1._wp
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+            &                                      poa_i, poa_i )
          !
          !== melt ponds ==!
