Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8503)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8504)
@@ -36,6 +36,4 @@
    PUBLIC   ice_adv    ! called by icestp
 
-   INTEGER ::   ncfl   ! number of ice time step with CFL>1/2  
-
    !! * Substitution
 #  include "vectopt_loop_substitute.h90"
@@ -51,5 +49,5 @@
       !!                   ***  ROUTINE ice_adv ***
       !!                    
-      !! ** purpose : advection/diffusion process of sea ice
+      !! ** purpose : advection of sea ice
       !!
       !! ** method  : variables included in the process are scalar,   
@@ -72,30 +70,14 @@
       REAL(wp), DIMENSION(jpi,jpj)           ::   zatold, zeiold, zesold, zsmvold 
       REAL(wp), DIMENSION(jpi,jpj,jpl)       ::   zhimax, zviold, zvsold
-      ! --- ultimate macho only --- !
-      REAL(wp) ::   zdt
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)      ::   zudy, zvdx, zcu_box, zcv_box
-      ! --- prather only --- !
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zarea
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0opw
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
-      ! MV MP 2016
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp
-      ! END MV MP 2016
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei
       !!---------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('iceadv')
  
       IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)''
-         WRITE(numout,*)'iceadv : sea-ice advection'
-         WRITE(numout,*)'~~~~~~~'
-         ncfl = 0                ! nb of time step with CFL > 1/2
+         WRITE(numout,*)
+         WRITE(numout,*) 'iceadv: sea-ice advection'
+         WRITE(numout,*) '~~~~~~'
       ENDIF
       
       CALL ice_var_agg( 1 ) ! integrated values + ato_i
-
-      !-------------------------------------!
-      !   Advection of sea ice properties   !
-      !-------------------------------------!
 
       ! conservation test
@@ -103,11 +85,11 @@
       
       ! store old values for diag
-      zviold = v_i
-      zvsold = v_s
-      zsmvold(:,:) = SUM( smv_i(:,:,:), dim=3 )
-      zeiold (:,:) = et_i
-      zesold (:,:) = et_s 
-
-      !--- Thickness correction init. --- !
+      zviold (:,:,:) = v_i(:,:,:)
+      zvsold (:,:,:) = v_s(:,:,:)
+      zsmvold(:,:)   = SUM( smv_i(:,:,:), dim=3 )
+      zeiold (:,:)   = et_i(:,:)
+      zesold (:,:)   = et_s(:,:)
+
+      ! Thickness correction init.
       zatold(:,:) = at_i
       WHERE( a_i(:,:,:) >= epsi20 )
@@ -119,5 +101,5 @@
       END WHERE
          
-      ! --- Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick --- !
+      ! Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick
       zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
       DO jl = 1, jpl
@@ -130,263 +112,28 @@
       END DO
       CALL lbc_lnk( zhimax(:,:,:), 'T', 1. )
-         
-      ! --- If ice drift field is too fast, use an appropriate time step for advection --- !        
-      zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )    ! CFL test for stability
-      zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
-      IF( lk_mpp )   CALL mpp_max( zcfl )
-      
-      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
-      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
-      ENDIF
-      
-!!      IF( zcfl > 0.5_wp .AND. lwp ) THEN
-!!         ncfl = ncfl + 1
-!!         IF( ncfl > 0 ) THEN   
-!!            WRITE(cltmp,'(i6.1)') ncfl
-!!            CALL ctl_warn( 'ice_adv: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
-!!         ENDIF
-!!      ENDIF
-
+      
+      !----------
+      ! Advection
+      !----------
       SELECT CASE ( nn_limadv )
-         
-                       !=============================!
-      CASE ( 0 )       !==  Ultimate-MACHO scheme  ==!                   
-                       !=============================!
-      
-         ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) )
-      
-         IF( kt == nit000 .AND. lwp ) THEN
-            WRITE(numout,*)''
-            WRITE(numout,*)'ice_adv_umx : Ultimate-MACHO advection scheme'
-            WRITE(numout,*)'~~~~~~~~~~~'
-         ENDIF
-         !
-         zdt = rdt_ice / REAL(initad)
-         
-         ! transport
-         zudy(:,:) = u_ice(:,:) * e2u(:,:)
-         zvdx(:,:) = v_ice(:,:) * e1v(:,:)
-         
-         ! define velocity for advection: u*grad(H)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               IF    ( u_ice(ji,jj) * u_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
-               ELSEIF( u_ice(ji,jj)                  >  0._wp ) THEN   ;   zcu_box(ji,jj) = u_ice(ji-1,jj)
-               ELSE                                                    ;   zcu_box(ji,jj) = u_ice(ji  ,jj)
-               ENDIF
-               
-               IF    ( v_ice(ji,jj) * v_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
-               ELSEIF( v_ice(ji,jj)                  >  0._wp ) THEN   ;   zcv_box(ji,jj) = v_ice(ji,jj-1)
-               ELSE                                                    ;   zcv_box(ji,jj) = v_ice(ji,jj  )
-               ENDIF
-            END DO
-         END DO
-         
-         ! advection
-         DO jt = 1, initad
-            CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, ato_i(:,:) )       ! Open water area 
-            DO jl = 1, jpl
-               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_i(:,:,jl) )      ! Ice area
-               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_i(:,:,jl) )      ! Ice  volume
-               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, smv_i(:,:,jl) )    ! Salt content
-               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, oa_i (:,:,jl) )    ! Age content
-               DO jk = 1, nlay_i
-                  CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_i(:,:,jk,jl) )   ! Ice  heat content
-               END DO
-               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_s(:,:,jl) )      ! Snow volume
-               CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_s(:,:,1,jl) )    ! Snow heat content
-               ! MV MP 2016
-               IF ( nn_pnd_scheme > 0 ) THEN
-                  CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_ip(:,:,jl) )  ! Melt pond fraction
-                  CALL ice_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_ip(:,:,jl) )  ! Melt pond volume
-               ENDIF
-               ! END MV MP 2016
-            END DO
-         END DO
-         !
-         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
-         DO jl = 2, jpl
-            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-         END DO
-         !
-         DEALLOCATE( zudy, zvdx, zcu_box, zcv_box )
-         
-                       !=============================!
-      CASE ( -1 )      !==     Prather scheme      ==!                   
-                       !=============================!
-
-         ALLOCATE( zarea(jpi,jpj)     , z0opw(jpi,jpj, 1 ) ,                                           &
-            &      z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) , z0ai(jpi,jpj,jpl) , z0es(jpi,jpj,jpl) ,   &
-            &      z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap(jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) ,   &
-            &      z0ei (jpi,jpj,nlay_i,jpl)                                                           )
-         
-         IF( kt == nit000 .AND. lwp ) THEN
-            WRITE(numout,*)''
-            WRITE(numout,*)'ice_adv_xy : Prather advection scheme'
-            WRITE(numout,*)'~~~~~~~~~~~'
-         ENDIF
-         
-         zarea(:,:) = e1e2t(:,:)
-         
-         !-------------------------
-         ! transported fields                                        
-         !-------------------------
-         z0opw(:,:,1) = ato_i(:,:) * e1e2t(:,:)             ! Open water area 
-         DO jl = 1, jpl
-            z0snw(:,:,jl) = v_s  (:,:,  jl) * e1e2t(:,:)  ! Snow volume
-            z0ice(:,:,jl) = v_i  (:,:,  jl) * e1e2t(:,:)  ! Ice  volume
-            z0ai (:,:,jl) = a_i  (:,:,  jl) * e1e2t(:,:)  ! Ice area
-            z0smi(:,:,jl) = smv_i(:,:,  jl) * e1e2t(:,:)  ! Salt content
-            z0oi (:,:,jl) = oa_i (:,:,  jl) * e1e2t(:,:)  ! Age content
-            z0es (:,:,jl) = e_s  (:,:,1,jl) * e1e2t(:,:)  ! Snow heat content
-            DO jk = 1, nlay_i
-               z0ei(:,:,jk,jl) = e_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
-            END DO
-            ! MV MP 2016
-            IF ( nn_pnd_scheme > 0 ) THEN
-               z0ap(:,:,jl)  = a_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction
-               z0vp(:,:,jl)  = v_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume
-            ENDIF
-            ! END MV MP 2016
-         END DO
-
-         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
-            DO jt = 1, initad
-               CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
-                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-               CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
-                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-               DO jl = 1, jpl
-                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
-                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
-                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
-                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
-                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
-                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
-                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
-                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
-                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
-                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
-                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
-                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
-                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-                  DO jk = 1, nlay_i                                                                !--- ice heat contents ---
-                     CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                     CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                  END DO
-                  ! MV MP 2016
-                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction --
-                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
-                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   --
-                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
-                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-                  ! END MV MP 2016
-               END DO
-            END DO
-         ELSE
-            DO jt = 1, initad
-               CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
-                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-               CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
-                  &                                         sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-               DO jl = 1, jpl
-                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
-                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
-                     &                                         sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
-                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
-                     &                                         sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
-                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
-                     &                                         sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
-                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
-                     &                                         sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
-                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
-                     &                                         sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
-                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-                  CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
-                     &                                         sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
-                     CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                     CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
-                        &                                         sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
-                        &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                  END DO
-                  ! MV MP 2016
-                  IF ( nn_pnd_scheme > 0 ) THEN
-                     CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction ---
-                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                     CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &
-                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                     CALL ice_adv_y( zusnit, v_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   ---
-                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-                     CALL ice_adv_x( zusnit, u_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &
-                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-                  ENDIF
-                  ! END MV MP 2016
-               END DO
-            END DO
-         ENDIF
-
-         !-------------------------------------------
-         ! Recover the properties from their contents
-         !-------------------------------------------
-         ato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:)
-         DO jl = 1, jpl
-            v_i  (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:)
-            v_s  (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:)
-            smv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:)
-            oa_i (:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:)
-            a_i  (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:)
-            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:)
-            DO jk = 1, nlay_i
-               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:)
-            END DO
-            ! MV MP 2016
-            IF ( nn_pnd_scheme > 0 ) THEN
-               a_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:)
-               v_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:)
-            ENDIF
-            ! END MV MP 2016
-         END DO
-         !
-         at_i(:,:) = a_i(:,:,1)      ! total ice fraction
-         DO jl = 2, jpl
-            at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-         END DO
-         !
-         DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0es , z0smi , z0oi , z0ap , z0vp , z0ei )
-         !
+      CASE ( 0 )                    !--  Ultimate-MACHO scheme
+         CALL ice_adv_umx( kt, u_ice, v_ice,  &
+            &              ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+         
+      CASE ( -1 )                   !-- Prather scheme
+         CALL ice_adv_prather( kt, u_ice, v_ice,  &
+            &                  ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+         
       END SELECT
-      
-      ! --- diags ---
+
+      ! total ice fraction
+      at_i(:,:) = a_i(:,:,1)
+      DO jl = 2, jpl
+         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
+      END DO
+      
+      !------------
+      ! diagnostics
+      !------------
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -404,12 +151,15 @@
       IF( iom_use('destrp') )   CALL iom_put( "destrp" , diag_trp_es         )         ! advected snw enthalpy (W/m2)
       
+      !--------------------------------------
+      ! Thickness correction in case too high
+      !--------------------------------------
       IF( nn_limdyn == 2 ) THEN
          !
-         CALL ice_var_zapsmall      !--- zap small areas ---!
+         CALL ice_var_zapsmall                       !-- zap small areas
          !
-         DO jl = 1, jpl             !--- Thickness correction in case too high --- !
+         DO jl = 1, jpl
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN
+                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to zhimax
                      !
                      ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl)
@@ -428,10 +178,10 @@
          END DO
                   
-         WHERE( ht_i(:,:,jpl) > hi_max(jpl) )             !--- bound ht_i to hi_max (99 m)
+         WHERE( ht_i(:,:,jpl) > hi_max(jpl) )        !-- bound ht_i to hi_max (99 m)
             ht_i(:,:,jpl) = hi_max(jpl)
             a_i (:,:,jpl) = v_i(:,:,jpl) / hi_max(jpl)
          END WHERE
 
-         IF ( nn_pnd_scheme > 0 ) THEN                    !--- correct pond fraction to avoid a_ip > a_i
+         IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i
             WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
          ENDIF
@@ -442,6 +192,5 @@
       ! Impose a_i < amax if no ridging/rafting or in mono-category
       !------------------------------------------------------------
-      !
-      IF( l_piling ) THEN ! simple conservative piling, comparable with 1-cat models
+      IF( l_piling ) THEN                            !-- simple conservative piling, comparable with 1-cat models
          at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
          DO jl = 1, jpl
@@ -452,5 +201,5 @@
       ENDIF
       
-      ! --- agglomerate variables -----------------
+      ! agglomerate variables 
       vt_i(:,:) = SUM( v_i(:,:,:), dim=3 )
       vt_s(:,:) = SUM( v_s(:,:,:), dim=3 )
@@ -464,5 +213,5 @@
       ! END MP 2016
       
-      ! --- open water = 1 if at_i=0 --------------------------------
+      ! open water = 1 if at_i=0 
       WHERE( at_i == 0._wp ) ato_i = 1._wp 
       
@@ -470,7 +219,7 @@
       IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
         
-      ! -------------------------------------------------
+      ! --------------
       ! control prints
-      ! -------------------------------------------------
+      ! --------------
       IF( ln_limctl )   CALL ice_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
       !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90	(revision 8503)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90	(revision 8504)
@@ -12,9 +12,9 @@
    !!   'key_lim3'                                       LIM3 sea-ice model
    !!----------------------------------------------------------------------
-   !!   ice_adv_x     : advection of sea ice on x axis
-   !!   ice_adv_y     : advection of sea ice on y axis
+   !!   ice_adv_prather     : advection of sea ice using Prather scheme
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean domain
    USE ice            ! sea-ice variables
+   USE sbc_oce , ONLY : nn_fsbc   ! frequency of sea-ice call
    !
    USE lbclnk         ! lateral boundary condition - MPP exchanges
@@ -27,6 +27,5 @@
    PRIVATE
 
-   PUBLIC   ice_adv_x   ! called by iceadv
-   PUBLIC   ice_adv_y   ! called by iceadv
+   PUBLIC   ice_adv_prather   ! called by iceadv
 
    !! * Substitutions
@@ -39,20 +38,226 @@
 CONTAINS
 
-   SUBROUTINE ice_adv_x( pdf, put , pcrh, psm , ps0 ,   &
-      &                  psx, psxx, psy , psyy, psxy )
+   SUBROUTINE ice_adv_prather( kt, pu_ice, pv_ice,  &
+      &                        pato_i, pv_i, pv_s, psmv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
       !!----------------------------------------------------------------------
-      !!                **  routine ice_adv_x  **
+      !!                **  routine ice_adv_prather  **
       !!  
       !! ** purpose :   Computes and adds the advection trend to sea-ice
-      !!              variable on i-axis
       !!
       !! ** method  :   Uses Prather second order scheme that advects tracers
-      !!              but also theirquadratic forms. The method preserves
-      !!              tracer structures by conserving second order moments.
+      !!                but also their quadratic forms. The method preserves
+      !!                tracer structures by conserving second order moments.
       !!
       !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
       !!----------------------------------------------------------------------
+      INTEGER                     , INTENT(in   ) ::   kt         ! time step
+      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
+      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
+      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psmv_i     ! salt content
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
+      !
+      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, z0es , z0smi , z0oi
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_adv_prather: Prather advection scheme'
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+      ENDIF
+      !
+      ALLOCATE( zarea(jpi,jpj)     , z0opw(jpi,jpj, 1 ) ,                                           &
+         &      z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) , z0ai(jpi,jpj,jpl) , z0es(jpi,jpj,jpl) ,   &
+         &      z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap(jpi,jpj,jpl) , z0vp(jpi,jpj,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(:,:) ) )
+      IF( lk_mpp )   CALL mpp_max( zcfl )
+      
+      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
+      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
+      ENDIF
+      
+      zarea(:,:) = e1e2t(:,:)
+      !-------------------------
+      ! transported fields                                        
+      !-------------------------
+      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) = psmv_i(:,:,  jl) * e1e2t(:,:)     ! Salt content
+         z0oi (:,:,jl) = poa_i (:,:,  jl) * e1e2t(:,:)     ! Age content
+         z0es (:,:,jl) = pe_s  (:,:,1,jl) * e1e2t(:,:)     ! Snow heat content
+         DO jk = 1, nlay_i
+            z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
+         END DO
+         IF ( nn_pnd_scheme > 0 ) THEN
+            z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
+            z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
+         ENDIF
+      END DO
+
+      !                                                    !--------------------------------------------!
+      IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
+         !                                                 !--------------------------------------------!
+         DO jt = 1, initad
+            CALL 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)  )
+               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
+                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
+               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
+                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
+               DO jk = 1, nlay_i                                                               !--- ice heat contents ---
+                  CALL 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 ( nn_pnd_scheme > 0 ) 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
+         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)  )
+               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
+                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
+               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es  (:,:,jl), sxc0 (:,:,jl),   &
+                  &                                      sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
+               DO jk = 1, nlay_i                                                             !--- ice heat contents ---
+                  CALL 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 ( nn_pnd_scheme > 0 ) 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
+         END DO
+      ENDIF
+
+      !-------------------------------------------
+      ! Recover the properties from their contents
+      !-------------------------------------------
+      pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:)
+      DO jl = 1, jpl
+         pv_i  (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:)
+         pv_s  (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:)
+         psmv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:)
+         poa_i (:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:)
+         pa_i  (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:)
+         pe_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e1e2t(:,:)
+         DO jk = 1, nlay_i
+            pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:)
+         END DO
+         ! MV MP 2016
+         IF ( nn_pnd_scheme > 0 ) THEN
+            pa_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:)
+            pv_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:)
+         ENDIF
+         ! END MV MP 2016
+      END DO
+      !
+      DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0es , z0smi , z0oi , z0ap , z0vp , z0ei )
+      !
+   END SUBROUTINE ice_adv_prather
+   
+   SUBROUTINE adv_x( pdf, put , pcrh, psm , ps0 ,   &
+      &              psx, psxx, psy , psyy, psxy )
+      !!----------------------------------------------------------------------
+      !!                **  routine adv_x  **
+      !!  
+      !! ** purpose :   Computes and adds the advection trend to sea-ice
+      !!                variable on x axis
+      !!----------------------------------------------------------------------
       REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
-      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call ice_adv_x then ice_adv_y (=1) or the opposite (=0)
+      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
@@ -209,29 +414,23 @@
 
       IF(ln_ctl) THEN
-         CALL prt_ctl(tab2d_1=psm  , clinfo1=' ice_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
-         CALL prt_ctl(tab2d_1=psx  , clinfo1=' ice_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
-         CALL prt_ctl(tab2d_1=psy  , clinfo1=' ice_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
-         CALL prt_ctl(tab2d_1=psxy , clinfo1=' ice_adv_x: psxy :')
+         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
       !
-   END SUBROUTINE ice_adv_x
-
-
-   SUBROUTINE ice_adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
-      &                  psx, psxx, psy , psyy, psxy )
+   END SUBROUTINE adv_x
+
+
+   SUBROUTINE adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
+      &              psx, psxx, psy , psyy, psxy )
       !!---------------------------------------------------------------------
-      !!                **  routine ice_adv_y  **
+      !!                **  routine adv_y  **
       !!            
       !! ** purpose :   Computes and adds the advection trend to sea-ice 
-      !!              variable on y axis
-      !!
-      !! ** method  :   Uses Prather second order scheme that advects tracers
-      !!              but also their quadratic forms. The method preserves 
-      !!              tracer structures by conserving second order moments.
-      !! 
-      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
+      !!                variable on y axis
       !!---------------------------------------------------------------------
       REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
-      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call ice_adv_x then ice_adv_y (=1) or the opposite (=0)
+      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
@@ -389,11 +588,11 @@
 
       IF(ln_ctl) THEN
-         CALL prt_ctl(tab2d_1=psm  , clinfo1=' ice_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
-         CALL prt_ctl(tab2d_1=psx  , clinfo1=' ice_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
-         CALL prt_ctl(tab2d_1=psy  , clinfo1=' ice_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
-         CALL prt_ctl(tab2d_1=psxy , clinfo1=' ice_adv_y: psxy :')
+         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
       !
-   END SUBROUTINE ice_adv_y
+   END SUBROUTINE adv_y
 
 #else
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90	(revision 8503)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90	(revision 8504)
@@ -43,7 +43,102 @@
 CONTAINS
 
-   SUBROUTINE ice_adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc )
+   SUBROUTINE ice_adv_umx( kt, pu_ice, pv_ice,  &
+      &                    pato_i, pv_i, pv_s, psmv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ice_adv_umx  ***
+      !! 
+      !! **  Purpose :   Compute the now trend due to total advection of 
+      !!                 tracers and add it to the general trend of tracer equations
+      !!                 using an "Ultimate-Macho" scheme
+      !!
+      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
+      !!----------------------------------------------------------------------
+      INTEGER                     , INTENT(in   ) ::   kt         ! time step
+      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
+      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
+      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psmv_i     ! salt content
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
+      !
+      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
+      INTEGER  ::   initad                  ! number of sub-timestep for the advection
+      REAL(wp) ::   zcfl , zusnit, zdt      !   -      -
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zudy, zvdx, zcu_box, zcv_box
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*) 
+         WRITE(numout,*) 'ice_adv_umx : Ultimate-MACHO advection scheme'
+         WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) )
+      !
+      ! --- 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(:,:) ) )
+      IF( lk_mpp )   CALL mpp_max( zcfl )
+
+      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
+      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
+      ENDIF
+
+      zdt = rdt_ice / REAL(initad)
+
+      ! --- transport --- !
+      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
+      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
+
+      ! --- define velocity for advection: u*grad(H) --- !
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
+            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
+            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
+            ENDIF
+
+            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
+            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
+            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
+            ENDIF
+         END DO
+      END DO
+
+      !---------------!
+      !== advection ==!
+      !---------------!
+      DO jt = 1, initad
+         CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area 
+         DO jl = 1, jpl
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, psmv_i(:,:,jl) )       ! Salt content
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i (:,:,jl) )       ! Age content
+            DO jk = 1, nlay_i
+               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content
+            END DO
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content
+            IF ( nn_pnd_scheme > 0 ) THEN
+               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction
+               CALL adv_umx( 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 )
+      !
+   END SUBROUTINE ice_adv_umx
+   
+   SUBROUTINE adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE adv_umx  ***
       !! 
       !! **  Purpose :   Compute the now trend due to total advection of 
@@ -146,5 +241,5 @@
       IF( nn_timing == 1 )  CALL timing_stop('ice_adv_umx')
       !
-   END SUBROUTINE ice_adv_umx
+   END SUBROUTINE adv_umx
 
 
