Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90	(revision 8553)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90	(revision 8554)
@@ -47,4 +47,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_opw_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_snw_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_dyn_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_rem_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_dif_1d
@@ -61,4 +62,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_sni_1d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_sum_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_dyn_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sub_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_sub_1d 
@@ -74,4 +76,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_res_1d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_spr_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_dyn_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_pnd_1d
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_bri_1d
@@ -84,7 +88,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_sub_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_lam_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_dyn_1d
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_1d        !: <==> the 2D  at_i
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ato_i_1d      !: <==> the 2D  ato_i
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhtur_1d      !: <==> the 2D  fhtur
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhld_1d       !: <==> the 2D  fhld
@@ -100,7 +106,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_i_1d        !: <==> the 2D  a_i
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ib_1d       !: <==> the 2D  a_i_b
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_i_1d       !: <==> the 2D  ht_i
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_ib_1d      !: <==> the 2D  ht_i_b
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_s_1d       !: <==> the 2D  ht_s
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_i_1d       !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_ib_1d      !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_s_1d       !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_su         !: Surface Conduction flux 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_bo_i       !: Bottom  Conduction flux 
@@ -117,4 +123,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_s_1d        !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   smv_i_1d      !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   oa_i_1d       !:
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d   !: corresponding to the 2D var  t_s
@@ -174,5 +181,5 @@
          &      t_bo_1d   (jpij) ,                                         &
          &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,    & 
-         &      hfx_dif_1d(jpij) , hfx_opw_1d(jpij) ,                      &
+         &      hfx_dif_1d(jpij) , hfx_opw_1d(jpij) ,hfx_dyn_1d(jpij) ,    &
          &      rn_amax_1d(jpij) ,                                         &
          &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) ,                      &
@@ -181,14 +188,14 @@
       !
       ii = ii + 1
-      ALLOCATE( sprecip_1d (jpij) , at_i_1d    (jpij) ,                     &
+      ALLOCATE( sprecip_1d (jpij) , at_i_1d    (jpij) , ato_i_1d(jpij) ,                    &
          &      fhtur_1d   (jpij) , wfx_snw_sni_1d (jpij) , wfx_spr_1d (jpij) , wfx_snw_sum_1d(jpij) , &
          &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) ,  &
          &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) ,  &
-         &      wfx_snw_sub_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) ,              &
-         &      wfx_lam_1d(jpij)  , dqns_ice_1d(jpij) , evap_ice_1d (jpij),                     &
+         &      wfx_snw_sub_1d(jpij), wfx_snw_dyn_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) ,   &
+         &      wfx_lam_1d(jpij)  , wfx_dyn_1d(jpij), wfx_pnd_1d(jpij),dqns_ice_1d(jpij) , evap_ice_1d (jpij),   &
          &      qprec_ice_1d(jpij), i0         (jpij) ,                                         &  
          &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  &
          &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  &
-         &      sfx_lam_1d (jpij) , hicol_1d   (jpij)                     , STAT=ierr(ii) )
+         &      sfx_lam_1d (jpij) , sfx_dyn_1d(jpij), hicol_1d   (jpij)    , STAT=ierr(ii) )
       !
       ii = ii + 1
@@ -197,5 +204,6 @@
          &      dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) ,                                  &    
          &      dh_i_bott(jpij) , dh_snowice(jpij) , sm_i_1d (jpij) , s_i_new(jpij) , &
-         &      a_ip_1d  (jpij) , v_ip_1d   (jpij) , v_i_1d  (jpij) , v_s_1d (jpij) , smv_i_1d(jpij) , STAT=ierr(ii) )
+         &      a_ip_1d  (jpij) , v_ip_1d   (jpij) , v_i_1d  (jpij) , v_s_1d (jpij) , &
+         &      smv_i_1d (jpij) , oa_i_1d   (jpij) , STAT=ierr(ii) )
       !
       ii = ii + 1
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 8553)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 8554)
@@ -21,4 +21,5 @@
    USE ice1D          ! sea-ice: thermodynamics
    USE ice            ! sea-ice: variables
+   USE icetab         ! sea-ice: 1D <==> 2D transformation
    USE icevar         ! sea-ice: operations
    USE icectl         ! sea-ice: control prints
@@ -39,17 +40,20 @@
 
    ! Variables shared among ridging subroutines
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apartf   ! participation function; fraction of ridging/closing associated w/ category n
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hi_hrdg  ! thickness of ridging ice / mean ridge thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   closing_net     ! net rate at which area is removed    (1/s)
+      !                                                 ! (ridging ice area - area of new ridges) / dt
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   opning         ! rate of opening due to divergence/shear
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   closing_gross  ! rate at which area removed, not counting area of new ridges
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apartf   ! participation function; fraction of ridging/closing associated w/ category n
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmin    ! minimum ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmax    ! maximum ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hraft    ! thickness of rafted ice
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hi_hrdg  ! thickness of ridging ice / mean ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aridge   ! participating ice ridging
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   araft    ! participating ice rafting
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_i_2d
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_s_2d
    !
    REAL(wp), PARAMETER ::   hrdg_hi_min = 1.1_wp    ! min ridge thickness multiplier: min(hrdg/hi)
    REAL(wp), PARAMETER ::   hi_hrft     = 0.5_wp    ! rafting multipliyer: (hi/hraft)
-   REAL(wp)            ::   zdrho                   ! 
    !
    ! ** namelist (namdyn_rdgrft) **
@@ -83,10 +87,11 @@
 
    INTEGER FUNCTION ice_dyn_rdgrft_alloc()
-      !!---------------------------------------------------------------------!
+      !!-------------------------------------------------------------------
       !!                ***  ROUTINE ice_dyn_rdgrft_alloc ***
-      !!---------------------------------------------------------------------!
-      ALLOCATE( asum (jpi,jpj)     , apartf (jpi,jpj,0:jpl) , aksum (jpi,jpj)     ,     &
-         &      hrmin(jpi,jpj,jpl) , hraft  (jpi,jpj,jpl)   , aridge(jpi,jpj,jpl) ,     &
-         &      hrmax(jpi,jpj,jpl) , hi_hrdg(jpi,jpj,jpl)   , araft (jpi,jpj,jpl) , STAT=ice_dyn_rdgrft_alloc )
+      !!-------------------------------------------------------------------
+      ALLOCATE( closing_net(jpij), opning(jpij)   , closing_gross(jpij),   &
+         &      apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl)    , aridge(jpij,jpl),   &
+         &      hrmax(jpij,jpl), hi_hrdg(jpij,jpl)  , araft (jpij,jpl),  &
+         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
 
       IF( lk_mpp                    )   CALL mpp_sum ( ice_dyn_rdgrft_alloc )
@@ -97,5 +102,5 @@
 
    SUBROUTINE ice_dyn_rdgrft( kt )
-      !!---------------------------------------------------------------------!
+      !!-------------------------------------------------------------------
       !!                ***  ROUTINE ice_dyn_rdgrft ***
       !!
@@ -121,19 +126,17 @@
       !!     This routine is based on CICE code and authors William H. Lipscomb,
       !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged
-      !!--------------------------------------------------------------------!
+      !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! number of iteration
       !!
-      INTEGER  ::   ji, jj, jk, jl     ! dummy loop index
-      INTEGER  ::   niter              ! local integer 
-      INTEGER  ::   iterate_ridging    ! if =1, repeat the ridging
-      REAL(wp) ::   zfac               ! local scalar
-      REAL(wp), DIMENSION(jpi,jpj) ::   closing_net     ! net rate at which area is removed    (1/s)
-      !                                                 ! (ridging ice area - area of new ridges) / dt
-      REAL(wp), DIMENSION(jpi,jpj) ::   divu_adv        ! divu as implied by transport scheme  (1/s)
-      REAL(wp), DIMENSION(jpi,jpj) ::   opning          ! rate of opening due to divergence/shear
-      REAL(wp), DIMENSION(jpi,jpj) ::   closing_gross   ! rate at which area removed, not counting area of new ridges
+      INTEGER  ::   ji, jj, jk, jl             ! dummy loop index
+      INTEGER  ::   niter, iterate_ridging     ! local integer 
+      INTEGER  ::   nidx2                      ! local integer
+      REAL(wp) ::   zfac                       ! local scalar
+      INTEGER , DIMENSION(jpij) ::   idxice2       ! compute ridge/raft or not
+      REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! divu as implied by transport scheme  (1/s)
+      REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
       !
       INTEGER, PARAMETER ::   nitermax = 20    
-      !!-----------------------------------------------------------------------------
+      !!-------------------------------------------------------------------
       ! controls
       IF( nn_timing == 1 )   CALL timing_start('icedyn_rdgrft')                                                             ! timing
@@ -143,18 +146,33 @@
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting'
-         IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~'
-         !
-         zdrho = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE
-         !
+         IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~'
       ENDIF      
 
-      !-----------------------------------------------------------------------------!
-      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
-      !-----------------------------------------------------------------------------!
-      !
-      CALL rdgrft_prep                             ! prepare ridging
-      !
-      DO jj = 1, jpj                                        ! Initialize arrays.
+      CALL ice_var_zapsmall   ! Zero out categories with very small areas
+
+      !--------------------------------
+      ! 0) Identify grid cells with ice
+      !--------------------------------
+      nidx = 0   ;   idxice(:) = 0
+      DO jj = 1, jpj
          DO ji = 1, jpi
+            IF ( SUM(a_i(ji,jj,:)) > 0._wp ) THEN
+               nidx           = nidx + 1
+               idxice( nidx ) = (jj - 1) * jpi + ji
+            ENDIF
+         END DO
+      END DO
+
+      IF( nidx > 0 ) THEN
+        
+         ! just needed here
+         CALL tab_2d_1d( nidx, idxice(1:nidx), zdivu(1:nidx), divu_i(:,:) )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), zdelt(1:nidx), delta_i(:,:) )
+         ! needed here and in the iteration loop
+         CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i(:,:,:) )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i(:,:,:) )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), ato_i_1d(1:nidx)      , ato_i(:,:) )
+
+         DO ji = 1, nidx
             !-----------------------------------------------------------------------------!
             ! 2) Dynamical inputs (closing rate, divu_adv, opning)
@@ -176,5 +194,5 @@
             !  (thick, newly ridged ice).
 
-            closing_net(ji,jj) = rn_csrdg * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
+            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
 
             ! 2.2 divu_adv
@@ -186,106 +204,142 @@
             ! If divu_adv < 0, make sure the closing rate is large enough
             ! to give asum = 1.0 after ridging.
-            
-            divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep
-
-            IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
+            zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice
+
+            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )
 
             ! 2.3 opning
             !------------
             ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
-            opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
-         END DO
-      END DO
-
-      !-----------------------------------------------------------------------------!
-      ! 3) Ridging iteration
-      !-----------------------------------------------------------------------------!
-      niter           = 1                 ! iteration counter
-      iterate_ridging = 1
-
-      DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
-
-         ! 3.1 closing_gross
+            opning(ji) = closing_net(ji) + zdivu_adv(ji)
+         END DO
+
+         !
+         !------------------------------------------------
+         ! 2.1) Identify grid cells with nonzero ridging
+         !------------------------------------------------
+         CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d )
+         nidx2 = 0 ; idxice2(:) = 0
+         DO ji = 1, nidx
+            IF( SUM( ABS(apartf(ji,1:jpl)) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
+               nidx2 = nidx2 + 1
+               idxice2(nidx2) = idxice(ji)
+               a_i_2d(nidx2,:) = a_i_2d(ji,:)  ! adjust to new indices
+               v_i_2d(nidx2,:) = v_i_2d(ji,:)  ! adjust to new indices
+               ato_i_1d(nidx2) = ato_i_1d(ji)  ! adjust to new indices
+            ENDIF
+         END DO
+         idxice(:) = idxice2(:)
+         nidx      = nidx2
+
+      ENDIF
+
+      !--------------
+      ! Start ridging
+      !--------------
+      IF( nidx > 0 ) THEN
+         
+         !-----------
+         ! 2D <==> 1D
+         !-----------
+         ! fields used but not modified
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d(1:nidx), sss_m(:,:) )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sst_1d(1:nidx), sst_m(:,:) )
+         ! the following fields are modified in this routine
+         !!CALL tab_2d_1d( nidx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) )
+         !!CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d(1:nidx,1:jpl), a_i(:,:,:) )
+         !!CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i  (:,:,:) )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s  (:,:,:) )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i (:,:,:) )
+         IF ( nn_pnd_scheme > 0 ) THEN
+            CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) )
+            CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) )
+         ENDIF
+         DO jl = 1, jpl
+            DO jk = 1, nlay_s
+               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) )
+            END DO
+            DO jk = 1, nlay_i
+               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_dyn_1d    (1:nidx), sfx_dyn    (:,:) )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bri_1d    (1:nidx), sfx_bri    (:,:) )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_dyn_1d    (1:nidx), wfx_dyn    (:,:) )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_dyn_1d    (1:nidx), hfx_dyn    (:,:) )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) )
+         IF ( nn_pnd_scheme > 0 ) THEN
+            CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) )
+         ENDIF
+         
          !-----------------------------------------------------------------------------!
-         ! Based on the ITD of ridging and ridged ice, convert the net
-         !  closing rate to a gross closing rate.  
-         ! NOTE: 0 < aksum <= 1
-         closing_gross(:,:) = closing_net(:,:) / aksum(:,:)
-
-         ! correction to closing rate and opening if closing rate is excessive
-         !---------------------------------------------------------------------
-         ! Reduce the closing rate if more than 100% of the open water 
-         ! would be removed.  Reduce the opening rate proportionately.
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zfac   = ( opning(ji,jj) - apartf(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice
-               IF    ( zfac < 0._wp .AND. zfac > - ato_i(ji,jj) ) THEN                  ! would lead to negative ato_i
-                  opning(ji,jj) = apartf(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 
-               ELSEIF( zfac > 0._wp .AND. zfac > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum
-                  opning(ji,jj) = apartf(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 
+         ! 3) Ridging iteration
+         !-----------------------------------------------------------------------------!
+         niter           = 1
+         iterate_ridging = 1      
+         DO WHILE( iterate_ridging > 0 .AND. niter < nitermax )
+
+            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d )
+
+            ! 3.2 Redistribute area, volume, and energy.
+            !-----------------------------------------------------------------------------!
+            CALL rdgrft_shift
+
+            ! 3.4 Do we keep on iterating?
+            !-----------------------------------------------------------------------------!
+            ! Check whether asum = 1.  If not (because the closing and opening
+            ! rates were reduced above), ridge again with new rates.
+            iterate_ridging = 0
+            DO ji = 1, nidx
+               zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) )
+               IF( ABS( zfac ) < epsi10 ) THEN
+                  closing_net(ji) = 0._wp
+                  opning     (ji) = 0._wp
+                  ato_i_1d   (ji) = MAX( 0._wp, 1._wp - SUM( a_i_2d(ji,:) ) )
+               ELSE
+                  iterate_ridging  = 1
+                  zdivu_adv  (ji) = zfac * r1_rdtice
+                  closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) )
+                  opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) )
                ENDIF
             END DO
-         END DO
-
-         ! correction to closing rate / opening if excessive ice removal
-         !---------------------------------------------------------------
-         ! Reduce the closing rate if more than 100% of any ice category 
-         ! would be removed.  Reduce the opening rate proportionately.
+            !
+            niter = niter + 1
+            IF( niter  >  nitermax )    CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' )
+            !
+         END DO
+
+         !-----------
+         ! 1D <==> 2D
+         !-----------
+         CALL tab_1d_2d( nidx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) )
+         CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d(1:nidx,1:jpl), a_i(:,:,:) )
+         CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i  (:,:,:) )
+         CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s  (:,:,:) )
+         CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) )
+         CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i (:,:,:) )
+         IF ( nn_pnd_scheme > 0 ) THEN
+            CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) )
+            CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) )
+         ENDIF
          DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zfac = apartf(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
-                  IF( zfac  >  a_i(ji,jj,jl) ) THEN
-                     closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / zfac
-                  ENDIF
-               END DO
-            END DO
-         END DO
-
-         ! 3.2 Redistribute area, volume, and energy.
-         !-----------------------------------------------------------------------------!
-         CALL rdgrft_shift( opning, closing_gross )
-         
-         ! 3.3 Compute total area of ice plus open water after ridging.
-         !-----------------------------------------------------------------------------!
-         ! This is in general not equal to one because of divergence during transport
-         asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
-
-         ! 3.4 Do we keep on iterating?
-         !-----------------------------------------------------------------------------!
-         ! Check whether asum = 1.  If not (because the closing and opening
-         ! rates were reduced above), ridge again with new rates.
-         iterate_ridging = 0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN
-                  closing_net(ji,jj) = 0._wp
-                  opning     (ji,jj) = 0._wp
-                  ato_i      (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) )
-               ELSE
-                  iterate_ridging    = 1
-                  divu_adv   (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice
-                  closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
-                  opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) )
-               ENDIF
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_max( iterate_ridging )
-
-         ! Repeat if necessary.
-         ! NOTE: If strength smoothing is turned on, the ridging must be
-         ! iterated globally because of the boundary update in the smoothing.
-         niter = niter + 1
-         !
-         IF( iterate_ridging == 1 ) THEN
-            CALL rdgrft_prep
-            IF( niter  >  nitermax ) THEN
-               WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
-               WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
-            ENDIF
+            DO jk = 1, nlay_s
+               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) )
+            END DO
+            DO jk = 1, nlay_i
+               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_dyn_1d    (1:nidx), sfx_dyn    (:,:) )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bri_1d    (1:nidx), sfx_bri    (:,:) )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_dyn_1d    (1:nidx), wfx_dyn    (:,:) )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_dyn_1d    (1:nidx), hfx_dyn    (:,:) )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) )
+         IF ( nn_pnd_scheme > 0 ) THEN
+            CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) )
          ENDIF
-         !
-      END DO !! on the do while over iter
-
+
+      ENDIF ! nidx > 0
+   
       CALL ice_var_agg( 1 ) 
 
@@ -298,26 +352,29 @@
 
 
-   SUBROUTINE rdgrft_prep
-      !!---------------------------------------------------------------------!
+   SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i )
+      !!-------------------------------------------------------------------
       !!                ***  ROUTINE rdgrft_prep ***
       !!
-      !! ** Purpose :   preparation for ridging and strength calculations
+      !! ** Purpose :   preparation for ridging calculations
       !!
       !! ** Method  :   Compute the thickness distribution of the ice and open water 
       !!              participating in ridging and of the resulting ridges.
-      !!---------------------------------------------------------------------!
-      INTEGER  ::   ji, jj, jl    ! dummy loop indices
-      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zdummy   ! local scalar
-      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   zGsum     ! zGsum(n) = sum of areas in categories 0 to n
-      !------------------------------------------------------------------------------!
-
-      z1_gstar      = 1._wp / rn_gstar    
-      z1_astar      = 1._wp / rn_astar    
-
-      CALL ice_var_zapsmall   ! Zero out categories with very small areas
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i 
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i, pv_i 
+      !!
+      INTEGER  ::   ji, jl                     ! dummy loop indices
+      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zfac   ! local scalar
+      REAL(wp), DIMENSION(jpij)        ::   zasum, z1_asum, zaksum   ! sum of a_i+ato_i and reverse 
+      REAL(wp), DIMENSION(jpij,jpl)    ::   zhi                      ! ice thickness
+      REAL(wp), DIMENSION(jpij,-1:jpl) ::   zGsum                    ! zGsum(n) = sum of areas in categories 0 to n
+      !--------------------------------------------------------------------
+
+      z1_gstar = 1._wp / rn_gstar
+      z1_astar = 1._wp / rn_astar
 
       !                       ! Ice thickness needed for rafting
-      WHERE( a_i(:,:,:) >= epsi20 )   ;   ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:)
-      ELSEWHERE                       ;   ht_i(:,:,:) = 0._wp
+      WHERE( pa_i(1:nidx,:) > 0._wp )   ;   zhi(1:nidx,:) = pv_i(1:nidx,:) / pa_i(1:nidx,:)
+      ELSEWHERE                         ;   zhi(1:nidx,:) = 0._wp
       END WHERE
 
@@ -338,31 +395,31 @@
       ! Compute total area of ice plus open water.
       ! This is in general not equal to one because of divergence during transport
-      asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
-      !
+      zasum(1:nidx) = pato_i(1:nidx) + SUM( pa_i(1:nidx,:), dim=2 )
+      !
+      WHERE( zasum(1:nidx) > 0._wp )   ;   z1_asum(1:nidx) = 1._wp / zasum(1:nidx)
+      ELSEWHERE                        ;   z1_asum(1:nidx) = 0._wp
+      END WHERE
       ! Compute cumulative thickness distribution function
       ! Compute the cumulative thickness distribution function zGsum,
       ! where zGsum(n) is the fractional area in categories 0 to n.
       ! initial value (in h = 0) equals open water area
-      zGsum(:,:,-1) = 0._wp
-      zGsum(:,:,0 ) = ato_i(:,:) / asum(:,:)
+      zGsum(1:nidx,-1) = 0._wp
+      zGsum(1:nidx,0 ) = pato_i(1:nidx) * z1_asum(1:nidx)
       DO jl = 1, jpl
-         zGsum(:,:,jl) = ( ato_i(:,:) + SUM( a_i(:,:,1:jl), dim=3 ) ) / asum(:,:) ! sum(1:jl) is ok (and not jpl)
+         zGsum(1:nidx,jl) = ( pato_i(1:nidx) + SUM( pa_i(1:nidx,1:jl), dim=2 ) ) * z1_asum(1:nidx)  ! sum(1:jl) is ok (and not jpl)
       END DO
-
       !
       IF( ln_partf_lin ) THEN          !--- Linear formulation (Thorndike et al., 1975)
          DO jl = 0, jpl    
-            DO jj = 1, jpj 
-               DO ji = 1, jpi
-                  IF    ( zGsum(ji,jj,jl)   < rn_gstar ) THEN
-                     apartf(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * &
-                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar )
-                  ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN
-                     apartf(ji,jj,jl) = z1_gstar * ( rn_gstar        - zGsum(ji,jj,jl-1) ) *  &
-                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar        ) * z1_gstar )
-                  ELSE
-                     apartf(ji,jj,jl) = 0._wp
-                  ENDIF
-               END DO
+            DO ji = 1, nidx
+               IF    ( zGsum(ji,jl)   < rn_gstar ) THEN
+                  apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * &
+                     &                       ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar )
+               ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN
+                  apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  &
+                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar )
+               ELSE
+                  apartf(ji,jl) = 0._wp
+               ENDIF
             END DO
          END DO
@@ -370,10 +427,14 @@
       ELSEIF( ln_partf_exp ) THEN      !--- Exponential, more stable formulation (Lipscomb et al, 2007)
          !                        
-         zdummy = 1._wp / ( 1._wp - EXP(-z1_astar) )        ! precompute exponential terms using zGsum as a work array
+         zfac = 1._wp / ( 1._wp - EXP(-z1_astar) )
          DO jl = -1, jpl
-            zGsum(:,:,jl) = EXP( -zGsum(:,:,jl) * z1_astar ) * zdummy
+            DO ji = 1, nidx
+               zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac
+            END DO
          END DO
          DO jl = 0, jpl
-            apartf(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl)
+            DO ji = 1, nidx
+               apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl)
+            END DO
          END DO
          !
@@ -383,20 +444,30 @@
       IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting
          DO jl = 1, jpl
-            DO jj = 1, jpj 
-               DO ji = 1, jpi
-                  aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jj,jl)
-                  araft (ji,jj,jl) = apartf(ji,jj,jl) - aridge(ji,jj,jl)
-               END DO
+            DO ji = 1, nidx
+               aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl)
+               araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl)
             END DO
          END DO
       ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
-         aridge(:,:,:) = apartf(:,:,1:jpl)
-         araft (:,:,:) = 0._wp
+         DO jl = 1, jpl
+            DO ji = 1, nidx
+               aridge(ji,jl) = apartf(ji,jl)
+               araft (ji,jl) = 0._wp
+            END DO
+         END DO
       ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
-        aridge(:,:,:) = 0._wp
-         araft (:,:,:) = apartf(:,:,1:jpl)
+         DO jl = 1, jpl
+            DO ji = 1, nidx
+               aridge(ji,jl) = 0._wp
+               araft (ji,jl) = apartf(ji,jl)
+            END DO
+         END DO
       ELSE                                               !- no ridging & no rafting
-         aridge(:,:,:) = 0._wp
-         araft (:,:,:) = 0._wp         
+         DO jl = 1, jpl
+            DO ji = 1, nidx
+               aridge(ji,jl) = 0._wp
+               araft (ji,jl) = 0._wp         
+            END DO
+         END DO
       ENDIF
 
@@ -422,40 +493,75 @@
       ! ridging.
       !
-      ! aksum = net area removed/ total area removed
+      ! zaksum = net area removed/ total area removed
       ! where total area removed = area of ice that ridges
       !         net area removed = total area removed - area of new ridges
       !-----------------------------------------------------------------
 
-      aksum(:,:) = apartf(:,:,0)
-      zdummy = 1._wp / hi_hrft
+      zaksum(:) = apartf(:,0)
+      zfac = 1._wp / hi_hrft
       ! Transfer function
       DO jl = 1, jpl !all categories have a specific transfer function
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF ( apartf(ji,jj,jl) > 0._wp ) THEN
-                  zhmean            = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * hrdg_hi_min )
-                  hrmin  (ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) )
-                  hrmax  (ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl)
-                  hraft  (ji,jj,jl) = ht_i(ji,jj,jl) * zdummy
-                  hi_hrdg(ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 )
-                  !
-                  ! Normalization factor : aksum, ensures mass conservation
-                  aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - hi_hrdg(ji,jj,jl) )    &
-                     &                        + araft (ji,jj,jl) * ( 1._wp - hi_hrft )
-               ELSE
-                  hrmin  (ji,jj,jl) = 0._wp 
-                  hrmax  (ji,jj,jl) = 0._wp 
-                  hraft  (ji,jj,jl) = 0._wp 
-                  hi_hrdg(ji,jj,jl) = 1._wp
-               ENDIF
-            END DO
+         DO ji = 1, nidx
+            IF ( apartf(ji,jl) > 0._wp ) THEN
+               zhmean         = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
+               hrmin  (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) )
+               hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
+               hraft  (ji,jl) = zhi(ji,jl) * zfac
+               hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
+               !
+               ! Normalization factor : zaksum, ensures mass conservation
+               zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) )    &
+                  &                    + araft (ji,jl) * ( 1._wp - hi_hrft )
+            ELSE
+               hrmin  (ji,jl) = 0._wp 
+               hrmax  (ji,jl) = 0._wp 
+               hraft  (ji,jl) = 0._wp 
+               hi_hrdg(ji,jl) = 1._wp
+            ENDIF
          END DO
       END DO
       !
+      ! 3.1 closing_gross
+      !-----------------------------------------------------------------------------!
+      ! Based on the ITD of ridging and ridged ice, convert the net
+      !  closing rate to a gross closing rate.  
+      ! NOTE: 0 < aksum <= 1
+      WHERE( zaksum(1:nidx) > 0._wp )   ;   closing_gross(1:nidx) = closing_net(1:nidx) / zaksum(1:nidx)
+      ELSEWHERE                         ;   closing_gross(1:nidx) = 0._wp
+      END WHERE
+      
+      DO ji = 1, nidx
+         
+         ! correction to closing rate and opening if closing rate is excessive
+         !---------------------------------------------------------------------
+         ! Reduce the closing rate if more than 100% of the open water 
+         ! would be removed.  Reduce the opening rate proportionately.
+         zfac = ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice
+         IF    ( zfac < 0._wp .AND. zfac > - pato_i(ji) ) THEN                  ! would lead to negative ato_i
+            opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice 
+         ELSEIF( zfac > 0._wp .AND. zfac > ( zasum(ji) - pato_i(ji) ) ) THEN  ! would lead to ato_i > asum
+            opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice 
+         ENDIF
+      END DO
+      
+      ! correction to closing rate / opening if excessive ice removal
+      !---------------------------------------------------------------
+      ! Reduce the closing rate if more than 100% of any ice category 
+      ! would be removed.  Reduce the opening rate proportionately.
+      DO jl = 1, jpl
+         DO ji = 1, nidx
+            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
+!!            IF( zfac > pa_i(ji,jl) .AND. zfac > 0._wp ) THEN
+            IF( zfac > pa_i(ji,jl) ) THEN
+               closing_gross(ji) = closing_gross(ji) * pa_i(ji,jl) / zfac
+            ENDIF
+         END DO
+      END DO      
+      !
    END SUBROUTINE rdgrft_prep
 
 
-   SUBROUTINE rdgrft_shift( opning, closing_gross )
-      !!----------------------------------------------------------------------
+   SUBROUTINE rdgrft_shift
+      !!-------------------------------------------------------------------
       !!                ***  ROUTINE rdgrft_shift ***
       !!
@@ -464,229 +570,163 @@
       !! ** Method  :   Remove area, volume, and energy from each ridging category
       !!              and add to thicker ice categories.
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   opning         ! rate of opening due to divergence/shear
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   closing_gross  ! rate at which area retreats, excluding area of new ridges
+      !!-------------------------------------------------------------------
       !
       INTEGER  ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
-      INTEGER  ::   ij                         ! horizontal index, combines i and j loops
-      INTEGER  ::   icells                     ! number of cells with a_i > puny
       REAL(wp) ::   hL, hR, farea              ! left and right limits of integration and new area going to jl2
-
-      INTEGER , DIMENSION(jpij) ::   indxi, indxj     ! compressed indices
+      REAL(wp) ::   vsw                        ! vol of water trapped into ridges
+      REAL(wp) ::   afrdg, afrft               ! fraction of category area ridged/rafted 
+      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
+      REAL(wp)                  ::   airft1, oirft1, aprft1
+      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, esrdg  ! area etc of new ridges
+      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, esrft  ! area etc of rafted ice
+      !
+      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
       REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2
-
-      REAL(wp), DIMENSION(jpij) ::   afrac            ! fraction of category area ridged 
-      REAL(wp), DIMENSION(jpij) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
-      REAL(wp), DIMENSION(jpij) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
-      ! MV MP 2016
-      REAL(wp), DIMENSION(jpij) ::   vprdg            ! pond volume of ridging ice
-      REAL(wp), DIMENSION(jpij) ::   aprdg1           ! pond area of ridging ice
-      REAL(wp), DIMENSION(jpij) ::   aprdg2           ! pond area of ridging ice
-      ! END MV MP 2016
-      REAL(wp), DIMENSION(jpij) ::   dhr, dhr2        ! hrmax - hrmin  &  hrmax^2 - hrmin^2
-
-      REAL(wp), DIMENSION(jpij) ::   vrdg1            ! volume of ice ridged
-      REAL(wp), DIMENSION(jpij) ::   vrdg2            ! volume of new ridges
-      REAL(wp), DIMENSION(jpij) ::   vsw              ! volume of seawater trapped into ridges
-      REAL(wp), DIMENSION(jpij) ::   srdg1            ! sal*volume of ice ridged
-      REAL(wp), DIMENSION(jpij) ::   srdg2            ! sal*volume of new ridges
-      REAL(wp), DIMENSION(jpij) ::   smsw             ! sal*volume of water trapped into ridges
-      REAL(wp), DIMENSION(jpij) ::   oirdg1, oirdg2   ! ice age of ice ridged
-
-      REAL(wp), DIMENSION(jpij) ::   afrft            ! fraction of category area rafted
-      REAL(wp), DIMENSION(jpij) ::   arft1 , arft2    ! area of ice rafted and new rafted zone
-      REAL(wp), DIMENSION(jpij) ::   virft , vsrft    ! ice & snow volume of rafting ice
-      ! MV MP 2016
-      REAL(wp), DIMENSION(jpij) ::   vprft            ! pond volume of rafting ice
-      REAL(wp), DIMENSION(jpij) ::   aprft1           ! pond area of rafted ice
-      REAL(wp), DIMENSION(jpij) ::   aprft2           ! pond area of new rafted ice
-      ! END MV MP 2016
-      REAL(wp), DIMENSION(jpij) ::   esrft , smrft    ! snow energy & salinity of rafting ice
-      REAL(wp), DIMENSION(jpij) ::   oirft1, oirft2   ! ice age of ice rafted
-
+      !
+      REAL(wp), DIMENSION(jpij,jpl) ::   z1_ai        ! 1 / a
+      !
       REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft     ! ice energy of rafting ice
-      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg1     ! enth*volume of ice ridged
-      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg2     ! enth*volume of new ridges
-      REAL(wp), DIMENSION(jpij,nlay_i) ::   ersw      ! enth of water trapped into ridges
-      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirdg     ! enth*volume of new ridges
+      !!-------------------------------------------------------------------
 
       !-------------------------------------------------------------------------------
       ! 1) Compute change in open water area due to closing and opening.
       !-------------------------------------------------------------------------------
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) +  &
-               &                     ( opning(ji,jj) - apartf(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
-         END DO
+      DO ji = 1, nidx
+         ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice )
       END DO
-
+      
       !-----------------------------------------------------------------
-      ! 2) Pump everything from ice which is being ridged / rafted
+      ! 2) compute categories participating to ridging/rafting (jl1) 
       !-----------------------------------------------------------------
-      ! Compute the area, volume, and energy of ice ridging in each
-      ! category, along with the area of the resulting ridge.
-
-      DO jl1 = 1, jpl !jl1 describes the ridging category
-
-         !------------------------------------------------
-         ! 2.1) Identify grid cells with nonzero ridging
-         !------------------------------------------------
-         icells = 0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( apartf(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN
-                  icells = icells + 1
-                  indxi(icells) = ji
-                  indxj(icells) = jj
+      DO jl1 = 1, jpl
+
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sm_i_1d(1:nidx), sm_i(:,:,jl1) )
+
+         DO ji = 1, nidx
+
+            !------------------------------------------------
+            ! 2.1) Identify grid cells with nonzero ridging
+            !------------------------------------------------
+            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
+
+               z1_ai(ji,jl1) = 1._wp / a_i_2d(ji,jl1)
+
+               !--------------------------------------------------------------------
+               ! 2.2) Compute area of ridging ice (airdg1) and of new ridge (airdg2)
+               !--------------------------------------------------------------------
+               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
+               airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice
+
+               airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
+               airft2(ji) = airft1 * hi_hrft
+
+               !---------------------------------------------------------------
+               ! 2.3) Compute ridging /rafting fractions, make sure afrdg <=1 
+               !---------------------------------------------------------------
+               afrdg = airdg1 * z1_ai(ji,jl1)
+               afrft = airft1 * z1_ai(ji,jl1)
+
+               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
+               vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg
+               ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
+
+               !---------------------------------------------------------------
+               ! 2.4) Compute ridging ice and new ridges (vi, vs, sm, oi, es, ei)
+               !---------------------------------------------------------------
+               virdg1     = v_i_2d  (ji,jl1)   * afrdg
+               virdg2(ji) = v_i_2d  (ji,jl1)   * afrdg * ( 1. + rn_porordg )
+               vsrdg(ji)  = v_s_2d  (ji,jl1)   * afrdg
+               sirdg1     = smv_i_2d(ji,jl1)   * afrdg
+               sirdg2(ji) = smv_i_2d(ji,jl1)   * afrdg + vsw * sss_1d(ji)
+               oirdg1     = oa_i_2d (ji,jl1)   * afrdg
+               oirdg2(ji) = oa_i_2d (ji,jl1)   * afrdg * hi_hrdg(ji,jl1) 
+               esrdg(ji)  = ze_s_2d (ji,1,jl1) * afrdg
+
+               virft(ji)  = v_i_2d  (ji,jl1)   * afrft
+               vsrft(ji)  = v_s_2d  (ji,jl1)   * afrft
+               sirft(ji)  = smv_i_2d(ji,jl1)   * afrft 
+               oirft1     = oa_i_2d (ji,jl1)   * afrft 
+               oirft2(ji) = oa_i_2d (ji,jl1)   * afrft * hi_hrft 
+               esrft(ji)  = ze_s_2d (ji,1,jl1) * afrft
+
+               !MV MP 2016
+               IF ( nn_pnd_scheme > 0 ) THEN
+                  aprdg1     = a_ip_2d(ji,jl1) * afrdg
+                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
+                  vprdg (ji) = v_ip_2d(ji,jl1) * afrdg
+                  aprft1     = a_ip_2d(ji,jl1) * afrft
+                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
+                  vprft (ji) = v_ip_2d(ji,jl1) * afrft
                ENDIF
-            END DO
-         END DO
-
-         DO ij = 1, icells
-            ji = indxi(ij) ; jj = indxj(ij)
-
-            !--------------------------------------------------------------------
-            ! 2.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
-            !--------------------------------------------------------------------
-            ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
-            arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
-
-            !---------------------------------------------------------------
-            ! 2.3) Compute ridging /rafting fractions, make sure afrac <=1 
-            !---------------------------------------------------------------
-            afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging
-            afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting
-            ardg2(ij) = ardg1(ij) * hi_hrdg(ji,jj,jl1)
-            arft2(ij) = arft1(ij) * hi_hrft
-
-            !--------------------------------------------------------------------------
-            ! 2.4) Substract area, volume, and energy from ridging 
-            !     / rafting category n1.
-            !--------------------------------------------------------------------------
-            vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij)
-            vrdg2(ij) = vrdg1(ij) * ( 1. + rn_porordg )
-            vsw  (ij) = vrdg1(ij) * rn_porordg
-
-            vsrdg (ij) = v_s  (ji,jj,  jl1) * afrac(ij)
-            esrdg (ij) = e_s  (ji,jj,1,jl1) * afrac(ij)
-            !MV MP 2016
-            IF ( nn_pnd_scheme > 0 ) THEN
-               aprdg1(ij) = a_ip(ji,jj, jl1) * afrac(ij)
-               aprdg2(ij) = a_ip(ji,jj, jl1) * afrac(ij) * hi_hrdg(ji,jj,jl1)
-               vprdg(ij)  = v_ip(ji,jj, jl1) * afrac(ij)
+               ! END MV MP 2016
+
+               !-----------------------------------------------------------------
+               ! 2.5) Ice-ocean exchanges
+               !-----------------------------------------------------------------
+               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
+               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_rdtice
+               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_rdtice          ! > 0 [W.m-2] 
+
+               ! Put the snow lost by ridging into the ocean
+               !------------------------------------------------
+               !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow.
+               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean
+                  &                                      + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
+
+               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2)
+                  &                                - esrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
+
+               ! Put the melt pond water into the ocean
+               !------------------------------------------            
+               IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN
+                  wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean
+                     &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice
+               ENDIF
+
+               ! virtual salt flux to keep salinity constant
+               !------------------------------------------            
+               IF( nn_icesal /= 2 )  THEN
+                  sirdg2(ji)      = sirdg2(ji)      - vsw * ( sss_1d(ji) - sm_i_1d(ji) )        ! ridge salinity = sm_i
+                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji)  * vsw * rhoic * r1_rdtice  &  ! put back sss_m into the ocean
+                     &                            - sm_i_1d(ji) * vsw * rhoic * r1_rdtice     ! and get  sm_i  from the ocean 
+               ENDIF
+
+
+               !------------------------------------------            
+               ! 2.6) update jl1 (removing ridged/rafted area)
+               !------------------------------------------            
+               a_i_2d  (ji,jl1) = a_i_2d  (ji,jl1) - airdg1    - airft1
+               v_i_2d  (ji,jl1) = v_i_2d  (ji,jl1) - virdg1    - virft(ji)
+               v_s_2d  (ji,jl1) = v_s_2d  (ji,jl1) - vsrdg(ji) - vsrft(ji)
+               smv_i_2d(ji,jl1) = smv_i_2d(ji,jl1) - sirdg1    - sirft(ji)
+               oa_i_2d (ji,jl1) = oa_i_2d (ji,jl1) - oirdg1    - oirft1
+               ! MV MP 2016
+               IF ( nn_pnd_scheme > 0 ) THEN
+                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
+                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
+               ENDIF
+               ! END MV MP 2016
+               ze_s_2d(ji,1,jl1) = ze_s_2d(ji,1,jl1) - esrdg (ji) - esrft (ji)
             ENDIF
-            ! END MV MP 2016
-            srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
-            oirdg1(ij) = oa_i (ji,jj,  jl1) * afrac(ij)
-            oirdg2(ij) = oa_i (ji,jj,  jl1) * afrac(ij) * hi_hrdg(ji,jj,jl1) 
-
-            ! rafting volumes, heat contents ...
-            virft (ij) = v_i  (ji,jj,  jl1) * afrft(ij)
-            vsrft (ij) = v_s  (ji,jj,  jl1) * afrft(ij)
-            !MV MP 2016
-            IF ( nn_pnd_scheme > 0 ) THEN
-               aprft1(ij) = a_ip (ji,jj,  jl1) * afrft(ij)
-               aprft2(ij) = a_ip (ji,jj,  jl1) * afrft(ij) * hi_hrft
-               vprft(ij)  = v_ip(ji,jj,jl1)    * afrft(ij)
-            ENDIF
-            ! END MV MP 2016
-            srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
-            esrft (ij) = e_s  (ji,jj,1,jl1) * afrft(ij)
-            smrft (ij) = smv_i(ji,jj,  jl1) * afrft(ij) 
-            oirft1(ij) = oa_i (ji,jj,  jl1) * afrft(ij) 
-            oirft2(ij) = oa_i (ji,jj,  jl1) * afrft(ij) * hi_hrft 
-
-            !-----------------------------------------------------------------
-            ! 2.5) Compute properties of new ridges
-            !-----------------------------------------------------------------
-            smsw(ij)  = vsw(ij) * sss_m(ji,jj)                   ! salt content of seawater frozen in voids
-            srdg2(ij) = srdg1(ij) + smsw(ij)                     ! salt content of new ridge
-            
-            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice
-            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
-
-            ! virtual salt flux to keep salinity constant
-            IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
-               srdg2(ij)      = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) )           ! ridge salinity = sm_i
-               sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj)    * vsw(ij) * rhoic * r1_rdtice  &  ! put back sss_m into the ocean
-                  &                            - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice     ! and get  sm_i  from the ocean 
-            ENDIF
-
-            !------------------------------------------            
-            ! 2.6 Put the snow somewhere in the ocean
-            !------------------------------------------            
-            !  Place part of the snow lost by ridging into the ocean. 
-            !  Note that esrdg > 0; the ocean must cool to melt snow.
-            !  If the ocean temp = Tf already, new ice must grow.
-            !  During the next time step, thermo_rates will determine whether
-            !  the ocean cools or new ice grows.
-            wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnwrdg )   & 
-               &                                      + rhosn * vsrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice  ! fresh water source for ocean
-
-            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnwrdg )         & 
-               &                                - esrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
-
-            ! MV MP 2016
-            !------------------------------------------            
-            ! 2.7 Put the melt pond water in the ocean
-            !------------------------------------------            
-            !  Place part of the melt pond volume into the ocean. 
-            IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN
-               wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( rhofw * vprdg(ij) * ( 1._wp - rn_fpndrdg )   & 
-               &                                 + rhofw * vprft(ij) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice  ! fresh water source for ocean
-            ENDIF
-            ! END MV MP 2016
-
-            !-----------------------------------------------------------------
-            ! 2.8 Compute quantities used to apportion ice among categories
-            ! in the n2 loop below
-            !-----------------------------------------------------------------
-            dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1)                    - hrmin(ji,jj,jl1)                    )
-            dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) )
-
-
-            ! update jl1 (removing ridged/rafted area)
-            a_i  (ji,jj,  jl1) = a_i  (ji,jj,  jl1) - ardg1 (ij) - arft1 (ij)
-            v_i  (ji,jj,  jl1) = v_i  (ji,jj,  jl1) - vrdg1 (ij) - virft (ij)
-            v_s  (ji,jj,  jl1) = v_s  (ji,jj,  jl1) - vsrdg (ij) - vsrft (ij)
-            e_s  (ji,jj,1,jl1) = e_s  (ji,jj,1,jl1) - esrdg (ij) - esrft (ij)
-            smv_i(ji,jj,  jl1) = smv_i(ji,jj,  jl1) - srdg1 (ij) - smrft (ij)
-            oa_i (ji,jj,  jl1) = oa_i (ji,jj,  jl1) - oirdg1(ij) - oirft1(ij)
-
-            ! MV MP 2016
-            IF ( nn_pnd_scheme > 0 ) THEN
-               v_ip (ji,jj,jl1) = v_ip (ji,jj,jl1) - vprdg (ij) - vprft (ij)
-               a_ip (ji,jj,jl1) = a_ip (ji,jj,jl1) - aprdg1(ij) - aprft1(ij)
-            ENDIF
-            ! END MV MP 2016
-
-         END DO
-
-         !--------------------------------------------------------------------
-         ! 2.9 Compute ridging ice enthalpy, remove it from ridging ice and
-         !      compute ridged ice enthalpy 
-         !--------------------------------------------------------------------
+
+         END DO ! ji
+
+         ! special loop for e_i because of layers jk
          DO jk = 1, nlay_i
-            DO ij = 1, icells
-               ji = indxi(ij) ; jj = indxj(ij)
-               ! heat content of ridged ice
-               erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij) 
-               eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)               
-               
-               ! enthalpy of the trapped seawater (J/m2, >0)
-               ! clem: if sst>0, then ersw <0 (is that possible?)
-               ersw(ij,jk)  = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i
-
-               ! heat flux to the ocean
-               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux 
-
-               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
-               erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk)
-
-               ! update jl1
-               e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
-
-            END DO
-         END DO
+            DO ji = 1, nidx
+               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
+                  ! Compute ridging /rafting fractions
+                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji,jl1)
+                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji,jl1)
+                  ! Compute ridging ice and new ridges for ei
+                  eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
+                  eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
+                  ! Update jl1
+                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
+               ENDIF
+            END DO
+         END DO
+         
 
          !-------------------------------------------------------------------------------
@@ -694,48 +734,55 @@
          !-------------------------------------------------------------------------------
          DO jl2  = 1, jpl 
-            ! over categories to which ridged/rafted ice is transferred
-            DO ij = 1, icells
-               ji = indxi(ij) ; jj = indxj(ij)
-
-               ! Compute the fraction of ridged ice area and volume going to thickness category jl2.
-               IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN
-                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
-                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
-                  farea    = ( hR      - hL      ) * dhr(ij) 
-                  fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij)
-               ELSE
-                  farea    = 0._wp 
-                  fvol(ij) = 0._wp                  
+            !
+            DO ji = 1, nidx
+
+               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
+
+                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2.
+                  IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
+                     !-----------------------------------------------------------------
+                     ! 2.8 Compute quantities used to apportion ice among categories
+                     ! in the n2 loop below
+                     !-----------------------------------------------------------------
+                     hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
+                     hR = MIN( hrmax(ji,jl1), hi_max(jl2)   )
+                     farea    = ( hR      - hL      ) / ( hrmax(ji,jl1)                 - hrmin(ji,jl1)                 )
+                     fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
+                  ELSE
+                     farea    = 0._wp 
+                     fvol(ji) = 0._wp                  
+                  ENDIF
+
+                  ! Compute the fraction of rafted ice area and volume going to thickness category jl2
+                  !!gm see above               IF( hraft(ji) <= hi_max(jl2) .AND. hraft(ji) >  hi_max(jl2-1) ) THEN
+                  IF( hi_max(jl2-1) < hraft(ji,jl1) .AND. hraft(ji,jl1) <= hi_max(jl2)  ) THEN   ;   zswitch(ji) = 1._wp
+                  ELSE                                                                           ;   zswitch(ji) = 0._wp
+                  ENDIF
+                  !
+                  a_i_2d  (ji,jl2) = a_i_2d  (ji,jl2) + ( airdg2(ji) * farea    + airft2(ji) * zswitch(ji) )
+                  oa_i_2d (ji,jl2) = oa_i_2d (ji,jl2) + ( oirdg2(ji) * farea    + oirft2(ji) * zswitch(ji) )
+                  v_i_2d  (ji,jl2) = v_i_2d  (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
+                  smv_i_2d(ji,jl2) = smv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
+                  v_s_2d  (ji,jl2) = v_s_2d  (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
+                     &                                    vsrft (ji) * rn_fsnwrft * zswitch(ji) )
+                  ! MV MP 2016
+                  IF ( nn_pnd_scheme > 0 ) THEN
+                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
+                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )
+                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         & 
+                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
+                  ENDIF
+                  ! END MV MP 2016
+                  ze_s_2d(ji,1,jl2) = ze_s_2d(ji,1,jl2) + ( esrdg (ji) * rn_fsnwrdg * fvol(ji)  +  &
+                     &                                      esrft (ji) * rn_fsnwrft * zswitch(ji) )
+                  
                ENDIF
 
-               ! Compute the fraction of rafted ice area and volume going to thickness category jl2
-!!gm see above               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
-               IF( hi_max(jl2-1) < hraft(ji,jj,jl1) .AND. hraft(ji,jj,jl1) <= hi_max(jl2)  ) THEN   ;   zswitch(ij) = 1._wp
-               ELSE                                                                                 ;   zswitch(ij) = 0._wp
-               ENDIF
-               !
-               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ( ardg2 (ij) * farea    + arft2 (ij) * zswitch(ij) )
-               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + ( oirdg2(ij) * farea    + oirft2(ij) * zswitch(ij) )
-               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) )
-               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) )
-               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + ( vsrdg (ij) * rn_fsnwrdg * fvol(ij)  +  &
-                  &                                        vsrft (ij) * rn_fsnwrft * zswitch(ij) )
-               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnwrdg * fvol(ij)  +  &
-                  &                                        esrft (ij) * rn_fsnwrft * zswitch(ij) )
-               ! MV MP 2016
-               IF ( nn_pnd_scheme > 0 ) THEN
-                  v_ip (ji,jj,jl2) = v_ip(ji,jj,jl2) + (   vprdg (ij) * rn_fpndrdg * fvol   (ij)   &
-                     &                                   + vprft (ij) * rn_fpndrft * zswitch(ij)   )
-                  a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2) + (   aprdg2(ij) * rn_fpndrdg * farea         & 
-                     &                                   + aprft2(ij) * rn_fpndrft * zswitch(ji)   )
-               ENDIF
-               ! END MV MP 2016
-            END DO
-
-            ! Transfer ice energy to category jl2 by ridging
+            END DO
+
             DO jk = 1, nlay_i
-               DO ij = 1, icells
-                  ji = indxi(ij) ; jj = indxj(ij)
-                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)                  
+               DO ji = 1, nidx
+                  IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
+                     &   ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)                  
                END DO
             END DO
@@ -743,5 +790,9 @@
          END DO ! jl2
          !
-      END DO ! jl1 (deforming categories)
+      END DO ! jl1
+      !
+      ! In case ridging/rafting lead to very small negative values (sometimes it happens)
+      WHERE( a_i_2d(1:nidx,:) < 0._wp )   a_i_2d(1:nidx,:) = 0._wp
+      WHERE( v_i_2d(1:nidx,:) < 0._wp )   v_i_2d(1:nidx,:) = 0._wp
       !
    END SUBROUTINE rdgrft_shift
@@ -760,5 +811,5 @@
       !!              ice strength has to be smoothed
       !!----------------------------------------------------------------------
-      INTEGER             ::   ji,jj, jl   ! dummy loop indices
+      INTEGER             ::   ji, jj, jl  ! dummy loop indices
       INTEGER             ::   ismooth     ! smoothing the resistance to deformation
       INTEGER             ::   itframe     ! number of time steps for the P smoothing
@@ -767,28 +818,7 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps
       !!----------------------------------------------------------------------
-
-      !                              !--------------------------------------------------!
-      CALL rdgrft_prep               ! Thickness distribution of ridging and ridged ice !
-      !                              !--------------------------------------------------!
-
       !                              !--------------------------------------------------!
       IF( ln_str_R75 ) THEN          ! Ice strength => Rothrock (1975) method           !
       !                              !--------------------------------------------------!
-         z1_3 = 1._wp / 3._wp
-         DO jl = 1, jpl
-            WHERE( apartf(:,:,jl) > 0._wp )
-               strength(:,:) =  -         apartf(:,:,jl) * ht_i   (:,:,jl) * ht_i(:,:,jl)   &  ! PE loss from deforming ice
-                  &             + 2._wp * araft (:,:,jl) * ht_i   (:,:,jl) * ht_i(:,:,jl)   &  ! PE gain from rafting ice
-                  &             +         aridge(:,:,jl) * hi_hrdg(:,:,jl) * z1_3 *         &  ! PE gain from ridging ice
-                  &                      ( hrmax(:,:,jl) * hrmax  (:,:,jl) +                &
-                  &                        hrmin(:,:,jl) * hrmin  (:,:,jl) +                &
-                  &                        hrmax(:,:,jl) * hrmin  (:,:,jl) )
-            ELSEWHERE
-               strength(:,:) = 0._wp
-            END WHERE
-         END DO
-         strength(:,:) = rn_perdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1)
-                         ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_perdg accounts for frictional dissipation
-         ismooth = 1
       !                              !--------------------------------------------------!
       ELSEIF( ln_str_H79 ) THEN      ! Ice strength => Hibler (1979) method             !
@@ -805,5 +835,5 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN 
+               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 
                   zworka(ji,jj) = ( 4.0 * strength(ji,jj)              &
                      &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &  
@@ -831,5 +861,5 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN 
+               IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 
                   itframe = 1 ! number of time steps for the running mean
                   IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icetab.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icetab.F90	(revision 8553)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icetab.F90	(revision 8554)
@@ -13,5 +13,4 @@
    !!   tab_1d_2d  : 1-D <==> 2-D
    !!----------------------------------------------------------------------
-   USE par_kind
    USE par_oce
    USE ice, ONLY : jpl
