Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 9451)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 9452)
@@ -106,11 +106,18 @@
       !!
       !! ** Method  :   Steps :
-      !!       1) Thickness categories boundaries, ice / o.w. concentrations
-      !!          Ridge preparation
-      !!       2) Dynamical inputs (closing rate, divu_adv, opning)
-      !!       3) Ridging iteration
-      !!       4) Ridging diagnostics
-      !!       5) Heat, salt and freshwater fluxes
-      !!       6) Compute increments of tate variables and come back to old values
+      !!       0) Identify grid cells with ice
+      !!       1) Calculate closing rate, divergence and opening
+      !!       2) Identify grid cells with ridging
+      !!       3) Start ridging iterations
+      !!          - prep = ridged and rafted ice + closing_gross
+      !!          - shift = move ice from one category to another
+      !!
+      !! ** Details
+      !!    step1: The net rate of closing is due to convergence and shear, based on Flato and Hibler (1995).
+      !!           The energy dissipation rate is equal to the net closing rate times the ice strength.
+      !!
+      !!    step3: The gross closing rate is equal to the first two terms (open
+      !!           water closing and thin ice ridging) without the third term
+      !!           (thick, newly ridged ice).
       !!
       !! References :   Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
@@ -123,5 +130,5 @@
       !!
       !!     This routine is based on CICE code and authors William H. Lipscomb,
-      !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged
+      !!     and Elizabeth C. Hunke, LANL are gratefully acknowledged
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! number of iteration
@@ -153,7 +160,8 @@
       !--------------------------------
       npti = 0   ;   nptidx(:) = 0
+      ipti = 0   ;   iptidx(:) = 0
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF ( SUM(a_i(ji,jj,:)) > 0._wp ) THEN
+            IF ( at_i(ji,jj) > 0._wp ) THEN
                npti           = npti + 1
                nptidx( npti ) = (jj - 1) * jpi + ji
@@ -161,5 +169,8 @@
          END DO
       END DO
-
+      
+      !--------------------------------------------------------
+      ! 1) Dynamical inputs (closing rate, divergence, opening)
+      !--------------------------------------------------------
       IF( npti > 0 ) THEN
         
@@ -173,118 +184,67 @@
 
          DO ji = 1, npti
-            !-----------------------------------------------------------------------------!
-            ! 2) Dynamical inputs (closing rate, divu_adv, opning)
-            !-----------------------------------------------------------------------------!
+            ! closing_net = rate at which open water area is removed + ice area removed by ridging 
+            !                                                        - ice area added in new ridges
+            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
             !
-            ! 2.1 closing_net
-            !-----------------
-            ! Compute the net rate of closing due to convergence 
-            ! and shear, based on Flato and Hibler (1995).
-            ! 
-            ! The energy dissipation rate is equal to the net closing rate
-            ! times the ice strength.
+            ! divergence given by the advection scheme
+            !   (which may not be equal to divu as computed from the velocity field)
+            zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice
             !
-            ! NOTE: The NET closing rate is equal to the rate that open water 
-            !  area is removed, plus the rate at which ice area is removed by 
-            !  ridging, minus the rate at which area is added in new ridges.
-            !  The GROSS closing rate is equal to the first two terms (open
-            !  water closing and thin ice ridging) without the third term
-            !  (thick, newly ridged ice).
-
-            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
-
-            ! 2.2 divu_adv
-            !--------------
-            ! Compute divu_adv, the divergence rate given by the transport/
-            ! advection scheme, which may not be equal to divu as computed 
-            ! from the velocity field.
-            !
-            ! If divu_adv < 0, make sure the closing rate is large enough
-            ! to give asum = 1.0 after ridging.
-            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.
+            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough
+            !                                                                                        ! to give asum = 1.0 after ridging
+            ! Opening rate (non-negative) that will give asum = 1.0 after ridging.
             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 )
-         ipti = 0 ; iptidx(:) = 0
+         !------------------------------------
+         ! 2) Identify grid cells with ridging
+         !------------------------------------
+         CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
+
          DO ji = 1, npti
-            IF( SUM( ABS(apartf(ji,1:jpl)) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
+            IF( SUM( apartf(ji,1:jpl) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
                ipti = ipti + 1
                iptidx     (ipti)   = nptidx     (ji)
-               a_i_2d     (ipti,:) = a_i_2d     (ji,:)  ! adjust to new indices
-               v_i_2d     (ipti,:) = v_i_2d     (ji,:)  ! adjust to new indices
-               ato_i_1d   (ipti)   = ato_i_1d   (ji)    ! adjust to new indices
-               closing_net(ipti)   = closing_net(ji)    ! adjust to new indices
-               zdivu_adv  (ipti)   = zdivu_adv  (ji)    ! adjust to new indices
-               opning     (ipti)   = opning     (ji)    ! adjust to new indices
+               ! adjust to new indices
+               a_i_2d     (ipti,:) = a_i_2d     (ji,:)
+               v_i_2d     (ipti,:) = v_i_2d     (ji,:)
+               ato_i_1d   (ipti)   = ato_i_1d   (ji)
+               closing_net(ipti)   = closing_net(ji)
+               zdivu_adv  (ipti)   = zdivu_adv  (ji)
+               opning     (ipti)   = opning     (ji)
             ENDIF
          END DO
-         nptidx(:) = iptidx(:)
-         npti      = ipti
 
       ENDIF
 
-      !--------------
-      ! Start ridging
-      !--------------
+      ! grid cells with ridging
+      nptidx(:) = iptidx(:)
+      npti      = ipti
+
+      !-----------------
+      ! 3) Start ridging
+      !-----------------
       IF( npti > 0 ) THEN
          
-         !-----------
-         ! 2D <==> 1D
-         !-----------
-         ! fields used but not modified
-         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
-         ! the following fields are modified in this routine
-         !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
-         !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
-         !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  (:,:,:) )
-         CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
-         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
-         CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
-         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
-         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
-         DO jl = 1, jpl
-            DO jk = 1, nlay_s
-               CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
-            END DO
-            DO jk = 1, nlay_i
-               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
-            END DO
-         END DO
-         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
-         
-         !-----------------------------------------------------------------------------!
-         ! 3) Ridging iteration
-         !-----------------------------------------------------------------------------!
+         CALL ice_dyn_1d2d( 1 )            ! --- Move to 1D arrays --- !
+
          iter            = 1
          iterate_ridging = 1      
-         DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax )
-
-            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d )
-
-            ! 3.2 Redistribute area, volume, and energy.
-            !-----------------------------------------------------------------------------!
+         !                                                        !----------------------!
+         DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax )  !  ridging iterations  !
+            !                                                     !----------------------!
+            ! Calculate participation function (apartf)
+            !       and transfer      function
+            !       and closing_gross (+correction on opening)
+            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
+
+            ! Redistribute area, volume, and energy between categories
             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.
+            ! Do we keep on iterating?
+            !-------------------------
+            ! Check whether a_i + ato_i = 0
+            ! If not, because the closing and opening rates were reduced above, ridge again with new rates
             iterate_ridging = 0
             DO ji = 1, npti
@@ -307,31 +267,7 @@
          END DO
 
-         !-----------
-         ! 1D <==> 2D
-         !-----------
-         CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
-         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
-         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
-         CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
-         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
-         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
-         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
-         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
-         DO jl = 1, jpl
-            DO jk = 1, nlay_s
-               CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
-            END DO
-            DO jk = 1, nlay_i
-               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
-            END DO
-         END DO
-         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
-         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
-         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
-         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
-         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
-         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
-
-      ENDIF ! npti > 0
+         CALL ice_dyn_1d2d( 2 )            ! --- Move to 2D arrays --- !
+
+      ENDIF
    
       CALL ice_var_agg( 1 ) 
@@ -345,5 +281,5 @@
 
 
-   SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i )
+   SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net )
       !!-------------------------------------------------------------------
       !!                ***  ROUTINE rdgrft_prep ***
@@ -352,7 +288,7 @@
       !!
       !! ** Method  :   Compute the thickness distribution of the ice and open water 
-      !!              participating in ridging and of the resulting ridges.
-      !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i 
+      !!                participating in ridging and of the resulting ridges.
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i, pclosing_net 
       REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i, pv_i 
       !!
@@ -372,15 +308,12 @@
       END WHERE
 
+      ! 1) Participation function (apartf): a(h) = b(h).g(h)
       !-----------------------------------------------------------------
-      ! 1) Participation function: a(h) = b(h).g(h) (apartf)
-      !-----------------------------------------------------------------
-      ! Compute the participation function apartf; this is analogous to
-      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
-      ! area lost from category n due to ridging/closing
-      ! apartf(n)   = total area lost due to ridging/closing
-      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). 
-      !
-      ! The expressions for apartf are found by integrating b(h)g(h) between
-      ! the category boundaries.
+      ! Compute the participation function = total area lost due to ridging/closing
+      ! This is analogous to
+      !   a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
+      !   assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar). 
+      !
+      ! apartf = integrating b(h)g(h) between the category boundaries
       ! apartf is always >= 0 and SUM(apartf(0:jpl))=1
       !-----------------------------------------------------------------
@@ -393,8 +326,9 @@
       ELSEWHERE                        ;   z1_asum(1:npti) = 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
+      ! initial value (in h = 0) = open water area
       zGsum(1:npti,-1) = 0._wp
       zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
@@ -465,5 +399,4 @@
       ENDIF
 
-      !-----------------------------------------------------------------
       ! 2) Transfer function
       !-----------------------------------------------------------------
@@ -479,10 +412,8 @@
       !  constrained by hrmin <= (zhmean + hi)/2.
       ! 
-      ! The maximum ridging thickness, hrmax, is determined by
-      !  zhmean and hrmin.
+      ! The maximum ridging thickness, hrmax, is determined by zhmean and hrmin.
       !
       ! These modifications have the effect of reducing the ice strength
-      ! (relative to the Hibler formulation) when very thick ice is
-      ! ridging.
+      ! (relative to the Hibler formulation) when very thick ice is ridging.
       !
       ! zaksum = net area removed/ total area removed
@@ -490,9 +421,8 @@
       !         net area removed = total area removed - area of new ridges
       !-----------------------------------------------------------------
-
       zfac = 1._wp / hi_hrft
       zaksum(1:npti) = apartf(1:npti,0)
-      ! Transfer function
-      DO jl = 1, jpl !all categories have a specific transfer function
+      !
+      DO jl = 1, jpl
          DO ji = 1, npti
             IF ( apartf(ji,jl) > 0._wp ) THEN
@@ -515,39 +445,37 @@
       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.  
+      ! 3) 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:npti) > 0._wp )   ;   closing_gross(1:npti) = closing_net(1:npti) / zaksum(1:npti)
+      WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
       ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp
       END WHERE
       
+      ! correction to closing rate if excessive ice removal
+      !----------------------------------------------------
+      ! Reduce the closing rate if more than 100% of any ice category would be removed
+      ! Reduce the opening rate in proportion
+      DO jl = 1, jpl
+         DO ji = 1, npti
+            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
+            IF( zfac > pa_i(ji,jl) ) THEN
+               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice
+            ENDIF
+         END DO
+      END DO      
+
+      ! 4) correction to opening if excessive open water removal
+      !---------------------------------------------------------
+      ! Reduce the closing rate if more than 100% of the open water would be removed
+      ! Reduce the opening rate in proportion
       DO ji = 1, npti  
-         ! 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
+         zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice
+         IF( zfac < 0._wp ) 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
+         ELSEIF( zfac > zasum(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, npti
-            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
@@ -561,5 +489,5 @@
       !!
       !! ** Method  :   Remove area, volume, and energy from each ridging category
-      !!              and add to thicker ice categories.
+      !!                and add to thicker ice categories.
       !!-------------------------------------------------------------------
       !
@@ -581,16 +509,16 @@
       REAL(wp), DIMENSION(jpij,nlay_s) ::   esrdg     ! enth*volume of new 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.
-      !-------------------------------------------------------------------------------
+      !
+      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation
+      !!-------------------------------------------------------------------
+
+      ! 1) Change in open water area due to closing and opening
+      !--------------------------------------------------------
       DO ji = 1, npti
          ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice )
       END DO
       
-      !-----------------------------------------------------------------
-      ! 2) compute categories participating to ridging/rafting (jl1) 
-      !-----------------------------------------------------------------
+      ! 2) compute categories in which ice is removed (jl1) 
+      !----------------------------------------------------
       DO jl1 = 1, jpl
 
@@ -599,14 +527,9 @@
          DO ji = 1, npti
 
-            !------------------------------------------------
-            ! 2.1) Identify grid cells with nonzero ridging
-            !------------------------------------------------
-            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
+            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
 
                z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
 
-               !--------------------------------------------------------------------
-               ! 2.2) Compute area of ridging ice (airdg1) and of new ridge (airdg2)
-               !--------------------------------------------------------------------
+               ! area of ridging / rafting 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
@@ -615,7 +538,5 @@
                airft2(ji) = airft1 * hi_hrft
 
-               !---------------------------------------------------------------
-               ! 2.3) Compute ridging /rafting fractions, make sure afrdg <=1 
-               !---------------------------------------------------------------
+               ! ridging /rafting fractions
                afrdg = airdg1 * z1_ai(ji)
                afrft = airft1 * z1_ai(ji)
@@ -625,7 +546,5 @@
                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)
-               !---------------------------------------------------------------
+               ! volume etc of ridging / rafting 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 )
@@ -651,7 +570,5 @@
                ENDIF
 
-               !-----------------------------------------------------------------
-               ! 2.5) Ice-ocean exchanges
-               !-----------------------------------------------------------------
+               ! Ice-ocean exchanges associated with ice porosity
                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
@@ -659,5 +576,4 @@
 
                ! 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
@@ -665,5 +581,4 @@
 
                ! Put the melt pond water into the ocean
-               !------------------------------------------            
                ! clem: I think the following lines must be commented since there
                !       is no net mass flux between melt ponds and the ocean (see icethd_pnd.F90 for ex.)
@@ -674,5 +589,4 @@
 
                ! virtual salt flux to keep salinity constant
-               !------------------------------------------            
                IF( nn_icesal /= 2 )  THEN
                   sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i
@@ -681,8 +595,6 @@
                ENDIF
 
-
-               !------------------------------------------            
-               ! 2.6) update jl1 (removing ridged/rafted area)
-               !------------------------------------------            
+               ! Remove area, volume of new ridge to each category jl1
+               !------------------------------------------------------
                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)
@@ -705,5 +617,5 @@
                   afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
                   afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
-                  ! Compute ridging ice and new ridges for es
+                  ! Compute ridging /rafting ice and new ridges for es
                   esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
                   esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
@@ -711,5 +623,7 @@
                   hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2)
                      &                                - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
-                  ! Update jl1
+                  !
+                  ! Remove energy of new ridge to each category jl1
+                  !-------------------------------------------------
                   ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
                ENDIF
@@ -727,5 +641,7 @@
                   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
+                  !
+                  ! Remove energy of new ridge to each category jl1
+                  !-------------------------------------------------
                   ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 
                ENDIF
@@ -733,8 +649,8 @@
          END DO
          
-
-         !-------------------------------------------------------------------------------
-         ! 3) Add area, volume, and energy of new ridge to each category jl2
-         !-------------------------------------------------------------------------------
+         ! 3) compute categories in which ice is added (jl2) 
+         !--------------------------------------------------
+         itest_rdg(1:npti) = 0
+         itest_rft(1:npti) = 0
          DO jl2  = 1, jpl 
             !
@@ -743,14 +659,12 @@
                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.
+                  ! 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) )
+                     !
+                     itest_rdg(ji) = 1   ! test for conservation
                   ELSE
                      farea    = 0._wp 
@@ -759,9 +673,21 @@
 
                   ! 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
+                  IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) >  hi_max(jl2-1) ) THEN
+                     zswitch(ji) = 1._wp
+                     !
+                     itest_rft(ji) = 1   ! test for conservation
+                  ELSE
+                     zswitch(ji) = 0._wp
                   ENDIF
                   !
+                  ! Patch to ensure perfect conservation if ice thickness goes mad
+                  ! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas)
+                  ! Then ice volume is removed from one category but the ridging/rafting scheme
+                  ! does not know where to move it, leading to a conservation issue.  
+                  IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN   ;   farea = 1._wp   ;   fvol(ji) = 1._wp   ;   ENDIF
+                  IF( itest_rft(ji) == 0 .AND. jl2 == jpl )      zswitch(ji) = 1._wp
+                  !
+                  ! Add area, volume of new ridge to category jl2
+                  !----------------------------------------------
                   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) )
@@ -780,5 +706,6 @@
 
             END DO
-            ! for snow enthalpy
+            ! Add snow energy of new ridge to category jl2
+            !---------------------------------------------
             DO jk = 1, nlay_s
                DO ji = 1, npti
@@ -788,5 +715,6 @@
                END DO
             END DO
-            ! for ice enthalpy
+            ! Add ice energy of new ridge to category jl2
+            !--------------------------------------------
             DO jk = 1, nlay_i
                DO ji = 1, npti
@@ -886,4 +814,77 @@
    END SUBROUTINE ice_strength
 
+   
+   SUBROUTINE ice_dyn_1d2d( kn )
+      !!-----------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_dyn_1d2d *** 
+      !!                 
+      !! ** Purpose :   move arrays from 1d to 2d and the reverse
+      !!-----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kn   ! 1= from 2D to 1D   ;   2= from 1D to 2D
+      !
+      INTEGER ::   jl, jk   ! dummy loop indices
+      !!-----------------------------------------------------------------------
+      !
+      SELECT CASE( kn )
+      !                    !---------------------!
+      CASE( 1 )            !==  from 2D to 1D  ==!
+         !                 !---------------------!
+         ! fields used but not modified
+         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
+         ! the following fields are modified in this routine
+         !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
+         !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
+         !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  (:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
+         DO jl = 1, jpl
+            DO jk = 1, nlay_s
+               CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
+            END DO
+            DO jk = 1, nlay_i
+               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
+
+         !                 !---------------------!
+      CASE( 2 )            !==  from 1D to 2D  ==!
+         !                 !---------------------!
+         CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
+         DO jl = 1, jpl
+            DO jk = 1, nlay_s
+               CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
+            END DO
+            DO jk = 1, nlay_i
+               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
+         !
+      END SELECT
+      !
+   END SUBROUTINE ice_dyn_1d2d
+   
 
    SUBROUTINE ice_dyn_rdgrft_init
