Index: trunk/NEMO/LIM_SRC_2/limadv_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limadv_2.F90	(revision 1529)
+++ trunk/NEMO/LIM_SRC_2/limadv_2.F90	(revision 1530)
@@ -4,4 +4,10 @@
    !! LIM 2.0 sea-ice model : sea-ice advection
    !!======================================================================
+   !! History :  OPA  ! 2000-01 (LIM)  Original code
+   !!                 ! 2001-05 (G. Madec, R. Hordoir) Doctor norm
+   !!   NEMO     1.0  ! 2003-10 (C. Ethe) F90, module
+   !!             -   ! 2003-12 (R. Hordoir, G. Madec) mpp
+   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c.
+   !!--------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
@@ -11,31 +17,27 @@
    !!   lim_adv_y_2  : advection of sea ice on y axis
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE dom_oce
    USE dom_ice_2
-   USE in_out_manager  ! I/O manager
    USE ice_2
    USE lbclnk
-   USE prtctl          ! Print control
+   USE in_out_manager     ! I/O manager
+   USE prtctl             ! Print control
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_adv_x_2    ! called by lim_trp
-   PUBLIC lim_adv_y_2    ! called by lim_trp
-
+   PUBLIC   lim_adv_x_2   ! called by lim_trp
+   PUBLIC   lim_adv_y_2   ! called by lim_trp
+
+   REAL(wp)  ::   epsi20 = 1.e-20   ! constant values
+   REAL(wp)  ::   rzero  = 0.e0     !    -       -
+   REAL(wp)  ::   rone   = 1.e0     !    -       -
+   
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
-
-   !! * Module variables
-   REAL(wp)  ::            &  ! constant values
-      epsi20 = 1e-20    ,  &
-      rzero  = 0.e0     ,  &
-      rone   = 1.e0
-   !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2009) 
    !! $Id$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -48,46 +50,32 @@
       !!  
       !! ** purpose :   Computes and adds the advection trend to sea-ice
-      !!      variable on x axis
+      !!              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.
-      !!      Ref.: "Numerical Advection by Conservation of Second Order
-      !!      Moments", JGR, VOL. 91. NO. D6. PAGES 6671-6681. MAY 20, 1986
-      !!     
-      !! History :
-      !!        !  00-01 (LIM)
-      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
-      !!        !  03-10 (C. Ethe) F90, module
-      !!        !  03-12 (R. Hordoir, G. Madec) mpp
+      !! ** method  :   Uses Prather second order scheme that advects tracers
+      !!              but also theirquadratic forms. The method preserves
+      !!              tracer structures by conserving second order moments.
+      !!
+      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
       !!--------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp)                    , INTENT(in)  ::  &
-         pdf ,       &  ! ???
-         pcrh           ! = 1. : lim_adv_x is called before lim_adv_y
-         !              ! = 0. : lim_adv_x is called after  lim_adv_y
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  ::  &
-         put            ! i-direction ice velocity at ocean U-point (m/s)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::  & 
-         ps0 , psm , &  ! ???
-         psx , psy , &  ! ??? 
-         psxx, psyy, psxy
-
-      !! * Local declarations
-      INTEGER ::   ji, jj      ! dummy loop indices
-      REAL(wp)       ::  &
-         zrdt, zslpmax, ztemp, zin0,     &  ! temporary scalars
-         zs1max, zs1new, zs2new,         &  !    "         "
-         zalf, zalfq, zalf1, zalf1q,     &  !    "         "
-         zbt , zbt1                         !    "         "
-      REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace
-         zf0 , zfx , zfy , zbet,         &  !    "           "
-         zfxx, zfyy, zfxy,               &  !    "           "
-         zfm, zalg, zalg1, zalg1q           !    "           "
+      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
+      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_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
+      !! 
+      INTEGER  ::   ji, jj                               ! dummy loop indices
+      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
+      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
+      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
+      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.
+      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
 
       DO jj = 1, jpj
@@ -99,5 +87,5 @@
                &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
             zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
-
+            !
             ps0 (ji,jj) = zslpmax  
             psx (ji,jj) = zs1new      * zin0
@@ -120,4 +108,5 @@
             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)  ) )
@@ -127,5 +116,5 @@
             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)
@@ -181,7 +170,8 @@
             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)   &
+            !
+            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)
@@ -202,25 +192,21 @@
             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) )
+            !
+            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( psm , 'T',  1. )
-      CALL lbc_lnk( ps0 , 'T',  1. )
-      CALL lbc_lnk( psx , 'T', -1. )
-      CALL lbc_lnk( psxx, 'T',  1. )
-      CALL lbc_lnk( psy , 'T', -1. )
-      CALL lbc_lnk( psyy, 'T',  1. )
+      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
+      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
+      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
       CALL lbc_lnk( psxy, 'T',  1. )
 
@@ -231,5 +217,5 @@
          CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
       ENDIF
-
+      !
    END SUBROUTINE lim_adv_x_2
 
@@ -241,39 +227,27 @@
       !!            
       !! ** purpose :   Computes and adds the advection trend to sea-ice 
-      !!      variable on y axis
+      !!              variable on j-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.
+      !!              but also their quadratic forms. The method preserves 
+      !!              tracer structures by conserving second order moments.
       !!
-      !! History :
-      !!   1.0  !  00-01 (LIM)
-      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
-      !!   2.0  !  03-10 (C. Ethe) F90, module
-      !!        !  03-12 (R. Hordoir, G. Madec) mpp
+      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
       !!---------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp),                     INTENT(in)  :: &
-         pdf,        &  ! ???
-         pcrh           ! = 1. : lim_adv_x is called before lim_adv_y
-         !              ! = 0. : lim_adv_x is called after  lim_adv_y
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: &
-         pvt            ! j-direction ice velocity at ocean V-point (m/s)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: &
-         psm , ps0 , psx , psy,   &
-         psxx, psyy, psxy
-
-      !! * Local Variables
-      INTEGER  ::   ji, jj                ! dummy loop indices
-      REAL(wp) ::   &
-         zrdt, zslpmax, zin0, ztemp,  &   ! temporary scalars
-         zs1max, zs1new, zs2new,      &   !    "         "
-         zalf, zalfq, zalf1, zalf1q,  &   !    "         "
-         zbt , zbt1                       !
-      REAL(wp), DIMENSION(jpi,jpj)  :: &
-         zf0 , zfx , zfy ,             &  ! temporary workspace
-         zfxx, zfyy, zfxy,             &  !    "           "
-         zfm , zbet,                   &  !    "           "
-         zalg, zalg1, zalg1q              !    "           "
+      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
+      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_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
+      !!
+      INTEGER  ::   ji, jj                               ! dummy loop indices
+      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
+      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
+      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
       !---------------------------------------------------------------------
 
@@ -282,141 +256,137 @@
       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( rzero, 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) )  )
-             zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
-             ps0 (ji,jj) = zslpmax  
-             psx (ji,jj)  = psx (ji,jj) * zin0
-             psxx(ji,jj)  = psxx(ji,jj) * zin0
-             psy (ji,jj) = zs1new * zin0
-             psyy(ji,jj) = zs2new * zin0
-             psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
-          END DO
-       END DO
-
-       !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
-       psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
-
-       !  Calculate fluxes and moments between boxes j<-->j+1              
-       DO jj = 1, jpj
-          DO ji = 1, jpi
-             !  Flux from j to j+1 WHEN v GT 0   
-             zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) )
-             zalf         =  MAX( rzero, 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(rzero, -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  * zalfq * psyy(ji,jj+1)
-             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
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zslpmax = MAX( rzero, 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) )  )
+            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
+            !
+            ps0 (ji,jj) = zslpmax  
+            psx (ji,jj) = psx (ji,jj) * zin0
+            psxx(ji,jj) = psxx(ji,jj) * zin0
+            psy (ji,jj) = zs1new * zin0
+            psyy(ji,jj) = zs2new * zin0
+            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0
+         END DO
+      END DO
+
+      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
+      psm(:,:)  = MAX(  pcrh * area(:,:) + ( 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( rzero, SIGN( rone, pvt(ji,jj) ) )
+            zalf         =  MAX( rzero, 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(rzero, -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  * zalfq * psyy(ji,jj+1)
+            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
+      !  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( psm , 'T',  1. )
-      CALL lbc_lnk( ps0 , 'T',  1. )
-      CALL lbc_lnk( psx , 'T', -1. )
-      CALL lbc_lnk( psxx, 'T',  1. )
-      CALL lbc_lnk( psy , 'T', -1. )
-      CALL lbc_lnk( psyy, 'T',  1. )
+      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
+      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
+      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
       CALL lbc_lnk( psxy, 'T',  1. )
 
@@ -427,5 +397,5 @@
          CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
       ENDIF
-
+      !
    END SUBROUTINE lim_adv_y_2
 
@@ -434,11 +404,6 @@
    !!   Default option            Dummy module     NO LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_adv_x_2         ! Empty routine
-   END SUBROUTINE lim_adv_x_2
-   SUBROUTINE lim_adv_y_2         ! Empty routine
-   END SUBROUTINE lim_adv_y_2
-
 #endif
 
+   !!======================================================================
 END MODULE limadv_2
Index: trunk/NEMO/LIM_SRC_3/limadv.F90
===================================================================
--- trunk/NEMO/LIM_SRC_3/limadv.F90	(revision 1529)
+++ trunk/NEMO/LIM_SRC_3/limadv.F90	(revision 1530)
@@ -4,4 +4,7 @@
    !! LIM sea-ice model : sea-ice advection
    !!======================================================================
+   !! History :  LIM  ! 2008-03 (M. Vancoppenolle)  LIM-3 from LIM-2 code
+   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c.
+   !!--------------------------------------------------------------------
 #if defined key_lim3
    !!----------------------------------------------------------------------
@@ -11,10 +14,9 @@
    !!   lim_adv_y  : advection of sea ice on y axis
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE dom_oce
    USE dom_ice
    USE ice
+   USE lbclnk
    USE in_out_manager  ! I/O manager
-   USE lbclnk
    USE prtctl          ! Print control
 
@@ -22,20 +24,17 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_adv_x    ! called by lim_trp
-   PUBLIC lim_adv_y    ! called by lim_trp
+   PUBLIC   lim_adv_x   ! called by lim_trp
+   PUBLIC   lim_adv_y   ! called by lim_trp
+
+   REAL(wp)  ::   epsi20 = 1.e-20   ! constant values
+   REAL(wp)  ::   rzero  = 0.e0     !    -       -
+   REAL(wp)  ::   rone   = 1.e0     !    -       -
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
-
-   !! * Module variables
-   REAL(wp)  ::            &  ! constant values
-      epsi20 = 1e-20    ,  &
-      rzero  = 0.e0     ,  &
-      rone   = 1.e0
-   !!----------------------------------------------------------------------
-   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2005) 
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2009) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -48,46 +47,32 @@
       !!  
       !! ** purpose :   Computes and adds the advection trend to sea-ice
-      !!      variable on x axis
+      !!              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.
-      !!      Ref.: "Numerical Advection by Conservation of Second Order
-      !!      Moments", JGR, VOL. 91. NO. D6. PAGES 6671-6681. MAY 20, 1986
-      !!     
-      !! History :
-      !!        !  00-01 (LIM)
-      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
-      !!        !  03-10 (C. Ethe) F90, module
-      !!        !  03-12 (R. Hordoir, G. Madec) mpp
+      !! ** method  :   Uses Prather second order scheme that advects tracers
+      !!              but also theirquadratic forms. The method preserves
+      !!              tracer structures by conserving second order moments.
+      !!
+      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
       !!--------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp)                    , INTENT(in)  ::  &
-         pdf ,       &  ! ???
-         pcrh           ! = 1. : lim_adv_x is called before lim_adv_y
-      !              ! = 0. : lim_adv_x is called after  lim_adv_y
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  ::  &
-         put            ! i-direction ice velocity at ocean U-point (m/s)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::  & 
-         ps0 , psm , &  ! ???
-         psx , psy , &  ! ??? 
-         psxx, psyy, psxy
-
-      !! * Local declarations
-      INTEGER ::   ji, jj      ! dummy loop indices
-      REAL(wp)       ::  &
-         zrdt, zslpmax, ztemp, zin0,     &  ! temporary scalars
-         zs1max, zs1new, zs2new,         &  !    "         "
-         zalf, zalfq, zalf1, zalf1q,     &  !    "         "
-         zbt , zbt1                         !    "         "
-      REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace
-         zf0 , zfx , zfy , zbet,         &  !    "           "
-         zfxx, zfyy, zfxy,               &  !    "           "
-         zfm, zalg, zalg1, zalg1q           !    "           "
+      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
+      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_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
+      !! 
+      INTEGER  ::   ji, jj                               ! dummy loop indices
+      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
+      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
+      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
+      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.
+      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
 
       DO jj = 1, jpj
@@ -120,11 +105,12 @@
             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  * zalfq * psxx(ji,jj)
-            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)
+            !
+            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.
@@ -148,11 +134,12 @@
             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  * zalfq * psxx(ji+1,jj)
-            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)
+            !
+            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
@@ -162,4 +149,5 @@
             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) )
@@ -181,12 +169,13 @@
             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)   &
+            !
+            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)
+               &                                                + 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)
+               &                                                + 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)
@@ -201,14 +190,13 @@
             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) )  )
+            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) )
@@ -217,10 +205,7 @@
 
       !-- Lateral boundary conditions
-      CALL lbc_lnk( psm , 'T',  1. )
-      CALL lbc_lnk( ps0 , 'T',  1. )
-      CALL lbc_lnk( psx , 'T', -1. )
-      CALL lbc_lnk( psxx, 'T',  1. )
-      CALL lbc_lnk( psy , 'T', -1. )
-      CALL lbc_lnk( psyy, 'T',  1. )
+      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
+      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
+      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
       CALL lbc_lnk( psxy, 'T',  1. )
 
@@ -231,5 +216,5 @@
          CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
       ENDIF
-
+      !
    END SUBROUTINE lim_adv_x
 
@@ -241,39 +226,27 @@
       !!            
       !! ** purpose :   Computes and adds the advection trend to sea-ice 
-      !!      variable on y axis
+      !!              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.
+      !!              but also their quadratic forms. The method preserves 
+      !!              tracer structures by conserving second order moments.
+      !! 
+      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681.
+      !!---------------------------------------------------------------------
+      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
+      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call lim_adv_x then lim_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
       !!
-      !! History :
-      !!   1.0  !  00-01 (LIM)
-      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
-      !!   2.0  !  03-10 (C. Ethe) F90, module
-      !!        !  03-12 (R. Hordoir, G. Madec) mpp
-      !!---------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp),                     INTENT(in)  :: &
-         pdf,        &  ! ???
-         pcrh           ! = 1. : lim_adv_x is called before lim_adv_y
-      !              ! = 0. : lim_adv_x is called after  lim_adv_y
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: &
-         pvt            ! j-direction ice velocity at ocean V-point (m/s)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: &
-         psm , ps0 , psx , psy,   &
-         psxx, psyy, psxy
-
-      !! * Local Variables
-      INTEGER  ::   ji, jj                ! dummy loop indices
-      REAL(wp) ::   &
-         zrdt, zslpmax, zin0, ztemp,  &   ! temporary scalars
-         zs1max, zs1new, zs2new,      &   !    "         "
-         zalf, zalfq, zalf1, zalf1q,  &   !    "         "
-         zbt , zbt1                       !
-      REAL(wp), DIMENSION(jpi,jpj)  :: &
-         zf0 , zfx , zfy ,             &  ! temporary workspace
-         zfxx, zfyy, zfxy,             &  !    "           "
-         zfm , zbet,                   &  !    "           "
-         zalg, zalg1, zalg1q              !    "           "
+      INTEGER  ::   ji, jj                               ! dummy loop indices
+      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp, zin0   ! temporary scalars
+      REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
+      REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
       !---------------------------------------------------------------------
 
@@ -290,7 +263,8 @@
                &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
             zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
+            !
             ps0 (ji,jj) = zslpmax  
-            psx (ji,jj)  = psx (ji,jj) * zin0
-            psxx(ji,jj)  = psxx(ji,jj) * zin0
+            psx (ji,jj) = psx (ji,jj) * zin0
+            psxx(ji,jj) = psxx(ji,jj) * zin0
             psy (ji,jj) = zs1new * zin0
             psyy(ji,jj) = zs2new * zin0
@@ -300,10 +274,9 @@
 
       !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
-      psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
+      psm(:,:)  = MAX(  pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
 
       !  Calculate fluxes and moments between boxes j<-->j+1              
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !  Flux from j to j+1 WHEN v GT 0   
+      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
+         DO ji = 1, jpi
             zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) )
             zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
@@ -311,4 +284,5 @@
             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)  ) ) 
@@ -318,5 +292,5 @@
             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)
@@ -329,5 +303,5 @@
          END DO
       END DO
-
+      !
       DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
          DO ji = 1, jpi
@@ -339,11 +313,12 @@
             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  * zalfq * psyy(ji,jj+1)
-            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)
+            !
+            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
@@ -354,8 +329,9 @@
             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)
+            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) )
@@ -373,16 +349,14 @@
             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)
-
+            !
+            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)  
-
+               &                                               + 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)
+               &                                               + 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)
@@ -397,26 +371,21 @@
             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) )
+            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( psm , 'T',  1. )
-      CALL lbc_lnk( ps0 , 'T',  1. )
-      CALL lbc_lnk( psx , 'T', -1. )
-      CALL lbc_lnk( psxx, 'T',  1. )
-      CALL lbc_lnk( psy , 'T', -1. )
-      CALL lbc_lnk( psyy, 'T',  1. )
+      CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. )
+      CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes
+      CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. )
       CALL lbc_lnk( psxy, 'T',  1. )
 
@@ -427,18 +396,14 @@
          CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
       ENDIF
-
+      !
    END SUBROUTINE lim_adv_y
 
+
 #else
    !!----------------------------------------------------------------------
    !!   Default option            Dummy module         NO LIM sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_adv_x         ! Empty routine
-   END SUBROUTINE lim_adv_x
-   SUBROUTINE lim_adv_y           ! Empty routine
-   END SUBROUTINE lim_adv_y
-
 #endif
 
+   !!======================================================================
 END MODULE limadv
