Index: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5234)
+++ branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5443)
@@ -17,9 +17,10 @@
    USE dom_oce        ! ocean domain
    USE sbc_oce        ! ocean surface boundary condition
-   USE par_ice        ! ice parameter
    USE dom_ice        ! ice domain
    USE ice            ! ice variables
    USE limadv         ! ice advection
    USE limhdf         ! ice horizontal diffusion
+   USE limvar         ! 
+   !
    USE in_out_manager ! I/O manager
    USE lbclnk         ! lateral boundary conditions -- MPP exchanges
@@ -28,12 +29,14 @@
    USE prtctl         ! Print control
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-   USE limvar          ! clem for ice thickness correction
-   USE timing          ! Timing
+   USE timing         ! Timing
    USE limcons        ! conservation tests
+   USE limctl         ! control prints
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   lim_trp    ! called by ice_step
+   PUBLIC   lim_trp    ! called by sbcice_lim
+
+   INTEGER  ::   ncfl                 ! number of ice time step with CFL>1/2  
 
    !! * Substitution
@@ -58,27 +61,28 @@
       !! ** action :
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! number of iteration
+      INTEGER, INTENT(in) ::   kt           ! number of iteration
       !
-      INTEGER  ::   ji, jj, jk, jl, jn      ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
       INTEGER  ::   initad                  ! number of sub-timestep for the advection
       REAL(wp) ::   zcfl , zusnit           !   -      -
+      CHARACTER(len=80) ::   cltmp
       !
-      REAL(wp), POINTER, DIMENSION(:,:)      ::   zui_u, zvi_v, zsm, zs0at, zs0ow
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi
-      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold   ! old ice volume...
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zaiold, zhimax   ! old ice concentration and thickness
-      REAL(wp), POINTER, DIMENSION(:,:)      ::   zeiold, zesold   ! old enthalpies
-      REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei
-      !
-      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
+      REAL(wp), POINTER, DIMENSION(:,:)      ::   zsm
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw
+      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold, zsmvold  ! old ice volume...
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness
+      REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold   ! old concentration, enthalpies
+      REAL(wp) ::    zdv, zvi, zvs, zsmv, zes, zei
+      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
       !!---------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limtrp')
 
-      CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )
-      CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
-      CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, zs0e )
-
-      CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold )   ! clem
+      CALL wrk_alloc( jpi,jpj,            zsm, zatold, zeiold, zesold )
+      CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
+      CALL wrk_alloc( jpi,jpj,1,          z0opw )
+      CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
+      CALL wrk_alloc( jpi,jpj,jpl,        zhimax, zviold, zvsold, zsmvold )
 
       IF( numit == nstart .AND. lwp ) THEN
@@ -88,8 +92,9 @@
          ENDIF
          WRITE(numout,*) '~~~~~~~~~~~~'
+         ncfl = 0                ! nb of time step with CFL > 1/2
       ENDIF
+
+      zsm(:,:) = e12t(:,:)
       
-      zsm(:,:) = area(:,:)
-
       !                             !-------------------------------------!
       IF( ln_limdyn ) THEN          !   Advection of sea ice properties   !
@@ -97,27 +102,33 @@
 
          ! conservation test
-         IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-         ! mass and salt flux init (clem)
+         IF( ln_limdiahsb )   CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
+         ! mass and salt flux init
          zviold(:,:,:)  = v_i(:,:,:)
-         zeiold(:,:)  = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 
-         zesold(:,:)  = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 
-
-         !--- Thickness correction init. (clem) -------------------------------
-         CALL lim_var_glo2eqv
-         zaiold(:,:,:) = a_i(:,:,:)
+         zvsold(:,:,:)  = v_s(:,:,:)
+         zsmvold(:,:,:) = smv_i(:,:,:)
+         zeiold(:,:)    = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 
+         zesold(:,:)    = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 
+
+         !--- Thickness correction init. -------------------------------
+         zatold(:,:) = SUM( a_i(:,:,:), dim=3 )
+         DO jl = 1, jpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
+                  ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
+                  ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
+               END DO
+            END DO
+         END DO
          !---------------------------------------------------------------------
-         ! Record max of the surrounding ice thicknesses for correction in limupdate
+         ! Record max of the surrounding ice thicknesses for correction
          ! in case advection creates ice too thick.
          !---------------------------------------------------------------------
-         zhimax(:,:,:) = ht_i(:,:,:)
+         zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
          DO jl = 1, jpl
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) )
-                  !zhimax(ji,jj,jl) = ( ht_i(ji  ,jj  ,jl) * tmask(ji,  jj  ,1) + ht_i(ji-1,jj-1,jl) * tmask(ji-1,jj-1,1) + ht_i(ji+1,jj+1,jl) * tmask(ji+1,jj+1,1) &
-                  !     &             + ht_i(ji-1,jj  ,jl) * tmask(ji-1,jj  ,1) + ht_i(ji  ,jj-1,jl) * tmask(ji  ,jj-1,1) &
-                  !     &             + ht_i(ji+1,jj  ,jl) * tmask(ji+1,jj  ,1) + ht_i(ji  ,jj+1,jl) * tmask(ji  ,jj+1,1) &
-                  !     &             + ht_i(ji-1,jj+1,jl) * tmask(ji-1,jj+1,1) + ht_i(ji+1,jj-1,jl) * tmask(ji+1,jj-1,1) )
+                  zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) )
                END DO
             END DO
@@ -125,74 +136,80 @@
          END DO
          
+         !=============================!
+         !==      Prather scheme     ==!
+         !=============================!
+
+         ! 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 )   ncfl = ncfl + 1
+!!         IF( lwp ) THEN
+!!            IF( ncfl > 0 ) THEN   
+!!               WRITE(cltmp,'(i6.1)') ncfl
+!!               CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
+!!            ELSE
+!!            !  WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 '
+!!            ENDIF
+!!         ENDIF
+
          !-------------------------
          ! transported fields                                        
          !-------------------------
-         ! Snow vol, ice vol, salt and age contents, area
-         zs0ow(:,:) = ato_i(:,:) * area(:,:)               ! Open water area 
-         DO jl = 1, jpl
-            zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume
-            zs0ice(:,:,jl)   = v_i  (:,:,jl) * area(:,:)    ! Ice  volume
-            zs0a  (:,:,jl)   = a_i  (:,:,jl) * area(:,:)    ! Ice area
-            zs0sm (:,:,jl)   = smv_i(:,:,jl) * area(:,:)    ! Salt content
-            zs0oi (:,:,jl)   = oa_i (:,:,jl) * area(:,:)    ! Age content
-            zs0c0 (:,:,jl)   = e_s  (:,:,1,jl)              ! Snow heat content
-            zs0e  (:,:,:,jl) = e_i  (:,:,:,jl)              ! Ice  heat content
-         END DO
-
-         !--------------------------
-         ! Advection of Ice fields  (Prather scheme)                                            
-         !--------------------------
-         ! If ice drift field is too fast, use an appropriate time step for advection.         
-         ! CFL test for stability
-         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) )
-         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) )
-         IF(lk_mpp )   CALL mpp_max( zcfl )
-!!gm more readability:
-!         IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
-!         ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
-!         ENDIF
-!!gm end
-         initad = 1 + NINT( MAX( 0._wp, SIGN( 1._wp, zcfl-0.5 ) ) )
-         zusnit = 1.0 / REAL( initad ) 
-         IF( zcfl > 0.5 .AND. lwp )   &
-            WRITE(numout,*) 'lim_trp   : CFL violation at day ', nday, ', cfl = ', zcfl,   &
-               &                        ': the ice time stepping is split in two'
+         z0opw(:,:,1) = ato_i(:,:) * e12t(:,:)             ! Open water area 
+         DO jl = 1, jpl
+            z0snw (:,:,jl)  = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
+            z0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
+            z0ai  (:,:,jl)  = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
+            z0smi (:,:,jl)  = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
+            z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
+            z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
+            DO jk = 1, nlay_i
+               z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
+            END DO
+         END DO
+
 
          IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
-            DO jn = 1,initad
-               CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area
-                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-               CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:),   &
-                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
+            DO jt = 1, initad
+               CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
+                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
+               CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &
+                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
                DO jl = 1, jpl
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
                      &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &
                      &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---     
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
                      &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &
                      &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
                      &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
                      &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents ---
                      &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &
                      &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
-                     CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                  DO jk = 1, nlay_i                                                                !--- ice heat contents ---
+                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
                         &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
                         &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
                         &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
                         &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
@@ -201,40 +218,39 @@
             END DO
          ELSE
-            DO jn = 1, initad
-               CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area
-                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:),   &
-                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
+            DO jt = 1, initad
+               CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
+                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
+               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:),   &
+                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
                DO jl = 1, jpl
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
                      &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &
                      &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
                      &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &
                      &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
                      &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai  (:,:,jl), sxa  (:,:,jl),   &
                      &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-                  CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
+                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents ---
                      &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es  (:,:,jl), sxc0 (:,:,jl),   &
                      &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
                   DO jk = 1, nlay_i                                                           !--- ice heat contents ---
-                     CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
                         &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
                         &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl),   & 
                         &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
                         &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
@@ -247,37 +263,41 @@
          ! Recover the properties from their contents
          !-------------------------------------------
-         zs0ow(:,:) = zs0ow(:,:) / area(:,:)
-         DO jl = 1, jpl
-            zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:)
-            zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:)
-            zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:)
-            zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:)
-            zs0a  (:,:,jl) = zs0a  (:,:,jl) / area(:,:)
-            !
+         ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
+         DO jl = 1, jpl
+            v_i  (:,:,jl)   = z0ice(:,:,jl) * r1_e12t(:,:)
+            v_s  (:,:,jl)   = z0snw(:,:,jl) * r1_e12t(:,:)
+            smv_i(:,:,jl)   = z0smi(:,:,jl) * r1_e12t(:,:)
+            oa_i (:,:,jl)   = z0oi (:,:,jl) * r1_e12t(:,:)
+            a_i  (:,:,jl)   = z0ai (:,:,jl) * r1_e12t(:,:)
+            e_s  (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
+            DO jk = 1, nlay_i
+               e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
+            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
 
          !------------------------------------------------------------------------------!
-         ! 4) Diffusion of Ice fields                  
+         ! Diffusion of Ice fields                  
          !------------------------------------------------------------------------------!
 
+         !
          !--------------------------------
          !  diffusion of open water area
          !--------------------------------
-         zs0at(:,:) = zs0a(:,:,1)      ! total ice fraction
-         DO jl = 2, jpl
-            zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl)
-         END DO
-         !
          !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
          DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row
             DO ji = 1 , fs_jpim1   ! vector opt.
-               pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji  ,jj) ) ) )   &
-                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj)
-               pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji,jj  ) ) ) )   &
-                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj)
+               pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
+                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
+               pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
+                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
             END DO
          END DO
          !
-         CALL lim_hdf( zs0ow (:,:) )   ! Diffusion
+         CALL lim_hdf( ato_i (:,:) )
 
          !------------------------------------
@@ -288,93 +308,64 @@
             DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
                DO ji = 1 , fs_jpim1   ! vector opt.
-                  pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji  ,jj,jl) ) ) )   &
-                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
-                  pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji,jj  ,jl) ) ) )   &
-                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
-               END DO
-            END DO
-
-            CALL lim_hdf( zs0ice (:,:,jl) )
-            CALL lim_hdf( zs0sn  (:,:,jl) )
-            CALL lim_hdf( zs0sm  (:,:,jl) )
-            CALL lim_hdf( zs0oi  (:,:,jl) )
-            CALL lim_hdf( zs0a   (:,:,jl) )
-            CALL lim_hdf( zs0c0  (:,:,jl) )
+                  pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,jl) ) ) )   &
+                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
+                  pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj  ,jl) ) ) )   &
+                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
+               END DO
+            END DO
+
+            CALL lim_hdf( v_i  (:,:,  jl) )
+            CALL lim_hdf( v_s  (:,:,  jl) )
+            CALL lim_hdf( smv_i(:,:,  jl) )
+            CALL lim_hdf( oa_i (:,:,  jl) )
+            CALL lim_hdf( a_i  (:,:,  jl) )
+            CALL lim_hdf( e_s  (:,:,1,jl) )
             DO jk = 1, nlay_i
-               CALL lim_hdf( zs0e (:,:,jk,jl) )
+               CALL lim_hdf( e_i(:,:,jk,jl) )
             END DO
          END DO
 
          !------------------------------------------------------------------------------!
-         ! 5) Update and limit ice properties after transport                           
+         ! limit ice properties after transport                           
          !------------------------------------------------------------------------------!
-
-         !--------------------------------------------------
-         ! 5.1) Recover mean values over the grid squares.
-         !--------------------------------------------------
-         zs0at(:,:) = 0._wp
+!!gm & cr   :  MAX should not be active if adv scheme is positive !
          DO jl = 1, jpl
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zs0sn (ji,jj,jl) = MAX( 0._wp, zs0sn (ji,jj,jl) )
-                  zs0ice(ji,jj,jl) = MAX( 0._wp, zs0ice(ji,jj,jl) )
-                  zs0sm (ji,jj,jl) = MAX( 0._wp, zs0sm (ji,jj,jl) )
-                  zs0oi (ji,jj,jl) = MAX( 0._wp, zs0oi (ji,jj,jl) )
-                  zs0a  (ji,jj,jl) = MAX( 0._wp, zs0a  (ji,jj,jl) )
-                  zs0c0 (ji,jj,jl) = MAX( 0._wp, zs0c0 (ji,jj,jl) )
-                  zs0at (ji,jj)    = zs0at(ji,jj) + zs0a(ji,jj,jl)
-               END DO
-            END DO
-         END DO
-
-         !---------------------------------------------------------
-         ! 5.2) Update and mask variables
-         !---------------------------------------------------------
-         DO jl = 1, jpl          
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) )
-
-                  zvi  = zs0ice(ji,jj,jl)
-                  zvs  = zs0sn (ji,jj,jl)
-                  zes  = zs0c0 (ji,jj,jl)      
-                  zsmv = zs0sm (ji,jj,jl)
-                  !
-                  ! Remove very small areas
-                  v_s(ji,jj,jl)   = rswitch * zs0sn (ji,jj,jl) 
-                  v_i(ji,jj,jl)   = rswitch * zs0ice(ji,jj,jl)
-                  a_i(ji,jj,jl)   = rswitch * zs0a  (ji,jj,jl)
-                  e_s(ji,jj,1,jl) = rswitch * zs0c0 (ji,jj,jl)      
-                  ! Ice salinity and age
-                  IF(  num_sal == 2  ) THEN
-                     smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) )
-                  ENDIF
-                  oa_i(ji,jj,jl) = MAX( rswitch * zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl)
-
-                 ! Update fluxes
-                  wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 
-                  wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
-                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
-                  hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
-              END DO
-            END DO
-         END DO
-
-         DO jl = 1, jpl
+                  v_s  (ji,jj,jl)   = MAX( 0._wp, v_s  (ji,jj,jl) )
+                  v_i  (ji,jj,jl)   = MAX( 0._wp, v_i  (ji,jj,jl) )
+                  smv_i(ji,jj,jl)   = MAX( 0._wp, smv_i(ji,jj,jl) )
+                  oa_i (ji,jj,jl)   = MAX( 0._wp, oa_i (ji,jj,jl) )
+                  a_i  (ji,jj,jl)   = MAX( 0._wp, a_i  (ji,jj,jl) )
+                  e_s  (ji,jj,1,jl) = MAX( 0._wp, e_s  (ji,jj,1,jl) )
+               END DO
+            END DO
+
             DO jk = 1, nlay_i
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     rswitch          = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) )
-                     zei              = zs0e(ji,jj,jk,jl)      
-                     e_i(ji,jj,jk,jl) = rswitch * MAX( 0._wp, zs0e(ji,jj,jk,jl) )
-                     ! Update fluxes
-                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
-                  END DO !ji
-               END DO ! jj
-            END DO ! jk
-         END DO ! jl
-
-         !--- Thickness correction in case too high (clem) --------------------------------------------------------
-         CALL lim_var_glo2eqv
+                     e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) )
+                  END DO
+               END DO
+            END DO
+         END DO
+!!gm & cr 
+
+         ! --- diags ---
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice
+               diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice
+
+               diag_trp_vi (ji,jj) = SUM(   v_i(ji,jj,:) -  zviold(ji,jj,:) ) * r1_rdtice
+               diag_trp_vs (ji,jj) = SUM(   v_s(ji,jj,:) -  zvsold(ji,jj,:) ) * r1_rdtice
+               diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice
+            END DO
+         END DO
+
+         ! zap small areas
+         CALL lim_var_zapsmall
+
+         !--- Thickness correction in case too high --------------------------------------------------------
          DO jl = 1, jpl
             DO jj = 1, jpj
@@ -382,4 +373,9 @@
 
                   IF ( v_i(ji,jj,jl) > 0._wp ) THEN
+
+                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
+                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
+                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
+                     
                      zvi  = v_i  (ji,jj,jl)
                      zvs  = v_s  (ji,jj,jl)
@@ -387,48 +383,47 @@
                      zes  = e_s  (ji,jj,1,jl)
                      zei  = SUM( e_i(ji,jj,1:nlay_i,jl) )
-                     zdv  = v_i(ji,jj,jl) - zviold(ji,jj,jl)   
-                     !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl)   
-                     
-                     rswitch = 1._wp
-                     IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. &
-                        & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                          
-                        ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) )
-                        rswitch        = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )
-                        a_i(ji,jj,jl)  = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )
-                     ELSE
-                        ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) )
-                        rswitch        = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )
-                        a_i(ji,jj,jl)  = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )
+
+                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)  
+
+                     IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
+                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
+
+                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
+                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
+
+                        ! small correction due to *rswitch for a_i
+                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl)
+                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl)
+                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl)
+                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
+                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
+
+                        ! Update mass fluxes
+                        wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
+                        wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
+                        sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
+                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
+                        hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0
+
                      ENDIF
 
-                     ! small correction due to *rswitch for a_i
-                     v_i  (ji,jj,jl) = rswitch * v_i  (ji,jj,jl)
-                     v_s  (ji,jj,jl) = rswitch * v_s  (ji,jj,jl)
-                     smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl)
-                     e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl)
-                     e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
-
-                     ! Update mass fluxes
-                     wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
-                     wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
-                     sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
-                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
-                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
                   ENDIF
+
                END DO
             END DO
          END DO
          ! -------------------------------------------------
-
-         ! --- diags ---
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice
-               diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice
-
-               diag_trp_vi(ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice
-               diag_trp_vs(ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice
-            END DO
-         END DO
+         
+         !--------------------------------------
+         ! Impose a_i < amax in mono-category
+         !--------------------------------------
+         !
+         IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), rn_amax )
+               END DO
+            END DO
+         ENDIF
 
          ! --- agglomerate variables -----------------
@@ -436,23 +431,19 @@
          vt_s (:,:) = 0._wp
          at_i (:,:) = 0._wp
-         !
          DO jl = 1, jpl
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  !
-                  vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
-                  vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
-                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
-               END DO
-            END DO
-         END DO
-         ! -------------------------------------------------
-
-         ! open water
+                  vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl)
+                  vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl)
+                  at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
+               END DO
+            END DO
+         END DO
+
+         ! --- open water = 1 if at_i=0 --------------------------------
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ! open water = 1 if at_i=0
                rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
-               ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * zs0ow(ji,jj)
+               ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
             END DO
          END DO      
@@ -463,43 +454,17 @@
       ENDIF
 
-      IF(ln_ctl) THEN   ! Control print
-         CALL prt_ctl_info(' ')
-         CALL prt_ctl_info(' - Cell values : ')
-         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_trp  : cell area :')
-         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_trp  : at_i      :')
-         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_trp  : vt_i      :')
-         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_trp  : vt_s      :')
-         DO jl = 1, jpl
-            CALL prt_ctl_info(' ')
-            CALL prt_ctl_info(' - Category : ', ivar1=jl)
-            CALL prt_ctl_info('   ~~~~~~~~~~')
-            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_trp  : a_i      : ')
-            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_trp  : ht_i     : ')
-            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_trp  : ht_s     : ')
-            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_trp  : v_i      : ')
-            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_trp  : v_s      : ')
-            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_trp  : e_s      : ')
-            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_trp  : t_su     : ')
-            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_trp  : t_snow   : ')
-            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_trp  : sm_i     : ')
-            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_trp  : smv_i    : ')
-            DO jk = 1, nlay_i
-               CALL prt_ctl_info(' ')
-               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
-               CALL prt_ctl_info('   ~~~~~~~')
-               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_trp  : t_i      : ')
-               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_trp  : e_i      : ')
-            END DO
-         END DO
-      ENDIF
+      ! -------------------------------------------------
+      ! control prints
+      ! -------------------------------------------------
+      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
       !
-      CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )
-      CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
-      CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, zs0e )
-
-      CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zaiold, zhimax )   ! clem
+      CALL wrk_dealloc( jpi,jpj,            zsm, zatold, zeiold, zesold )
+      CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
+      CALL wrk_dealloc( jpi,jpj,1,          z0opw )
+      CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
+      CALL wrk_dealloc( jpi,jpj,jpl,        zviold, zvsold, zhimax, zsmvold )
       !
       IF( nn_timing == 1 )  CALL timing_stop('limtrp')
+
    END SUBROUTINE lim_trp
 
@@ -512,5 +477,4 @@
    END SUBROUTINE lim_trp
 #endif
-
    !!======================================================================
 END MODULE limtrp
