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 8486)
+++ 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
