Index: NEMO/trunk/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 11558)
+++ NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 11560)
@@ -278,4 +278,5 @@
       ! controls
       IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints
+      IF( ln_icectl    )   CALL ice_prt     (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ')                             ! prints
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
       IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
Index: NEMO/trunk/src/ICE/icevar.F90
===================================================================
--- NEMO/trunk/src/ICE/icevar.F90	(revision 11558)
+++ NEMO/trunk/src/ICE/icevar.F90	(revision 11560)
@@ -868,19 +868,8 @@
       !! ** Purpose :  converting 1-cat ice to jpl ice categories
       !!
-      !!                  ice thickness distribution follows a gaussian law
-      !!               around the concentration of the most likely ice thickness
-      !!                           (similar as iceistate.F90)
-      !!
-      !! ** Method:   Iterative procedure
-      !!                
-      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
-      !!
-      !!               2) Check whether the distribution conserves area and volume, positivity and
-      !!                  category boundaries
+      !!
+      !! ** Method:   ice thickness distribution follows a gamma function from Abraham et al. (2015)
+      !!              it has the property of conserving total concentration and volume
       !!              
-      !!               3) If not (input ice is too thin), the last category is empty and
-      !!                  the number of categories is reduced (jpl-1)
-      !!
-      !!               4) Iterate until ok (SUM(itest(:) = 4)
       !!
       !! ** Arguments : phti: 1-cat ice thickness
@@ -890,5 +879,7 @@
       !! ** Output    : jpl-cat 
       !!
-      !!  (Example of application: BDY forcings when input are cell averaged)  
+      !!  Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015.
+      !!               Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice.
+      !!               Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614 
       !!-------------------------------------------------------------------
       REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
@@ -897,16 +888,10 @@
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
       !
-      INTEGER , DIMENSION(4) ::   itest
-      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra
+      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra, z1_hti
       INTEGER  ::   ji, jk, jl
-      INTEGER  ::   idim, i_fill, jl0  
-      REAL(wp) ::   zarg, zV, zconv, zdh, zdv
-      !!-------------------------------------------------------------------
-      !
-      ! == thickness and concentration == !
-      ! distribution over the jpl ice categories
-      !    a gaussian distribution for ice concentration is used
-      !    then we check whether the distribution fullfills
-      !    volume and area conservation, positivity and ice categories bounds
+      INTEGER  ::   idim
+      REAL(wp) ::   zv, zdh
+      !!-------------------------------------------------------------------
+      !
       idim = SIZE( phti , 1 )
       !
@@ -915,89 +900,62 @@
       pa_i(1:idim,1:jpl) = 0._wp
       !
+      ALLOCATE( z1_hti(idim) )
+      WHERE( phti(:) /= 0._wp )   ;   z1_hti(:) = 1._wp / phti(:)
+      ELSEWHERE                   ;   z1_hti(:) = 0._wp
+      END WHERE
+      !
+      ! == thickness and concentration == !
+      ! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl)
+      DO jl = 1, jpl-1
+         DO ji = 1, idim
+            !
+            IF( phti(ji) > 0._wp ) THEN
+               ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
+               pa_i(ji,jl) = pati(ji) * z1_hti(ji) * (  ( phti(ji) + 2.*hi_max(jl-1) ) * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) &
+                  &                                   - ( phti(ji) + 2.*hi_max(jl  ) ) * EXP( -2.*hi_max(jl  )*z1_hti(ji) ) )
+               !
+               ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl)
+               zv = pati(ji) * z1_hti(ji) * (  ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl-1) + 2.*hi_max(jl-1)*hi_max(jl-1) ) &
+                  &                            * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) &
+                  &                          - ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl) + 2.*hi_max(jl)*hi_max(jl) ) &
+                  &                            * EXP(-2.*hi_max(jl)*z1_hti(ji)) )
+               ! thickness
+               IF( pa_i(ji,jl) > epsi06 ) THEN
+                  ph_i(ji,jl) = zv / pa_i(ji,jl)
+               ELSE
+                  ph_i(ji,jl) = 0.
+                  pa_i(ji,jl) = 0.
+               ENDIF
+            ENDIF
+            !
+         ENDDO
+      ENDDO
+      !
+      ! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity
       DO ji = 1, idim
          !
          IF( phti(ji) > 0._wp ) THEN
-            !
-            ! find which category (jl0) the input ice thickness falls into
-            jl0 = jpl
-            DO jl = 1, jpl
-               IF ( ( phti(ji) >= hi_max(jl-1) ) .AND. ( phti(ji) < hi_max(jl) ) ) THEN
-                  jl0 = jl
-                  CYCLE
-               ENDIF
-            END DO
-            !
-            itest(:) = 0
-            i_fill   = jpl + 1                                            !------------------------------------
-            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories
-               !                                                          !------------------------------------
-               i_fill = i_fill - 1
-               !
-               ph_i(ji,1:jpl) = 0._wp
-               pa_i(ji,1:jpl) = 0._wp
-               itest(:)       = 0      
-               !
-               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1
-                  ph_i(ji,1) = phti(ji)
-                  pa_i(ji,1) = pati (ji)
-               ELSE                         !-- case ice is thicker: fill categories >1
-                  ! thickness
-                  DO jl = 1, i_fill - 1
-                     ph_i(ji,jl) = hi_mean(jl)
-                  END DO
-                  !
-                  ! concentration
-                  pa_i(ji,jl0) = pati(ji) / SQRT(REAL(jpl))
-                  DO jl = 1, i_fill - 1
-                     IF ( jl /= jl0 ) THEN
-                        zarg        = ( ph_i(ji,jl) - phti(ji) ) / ( phti(ji) * 0.5_wp )
-                        pa_i(ji,jl) =   pa_i (ji,jl0) * EXP(-zarg**2)
-                     ENDIF
-                  END DO
-                  !
-                  ! last category
-                  pa_i(ji,i_fill) = pati(ji) - SUM( pa_i(ji,1:i_fill-1) )
-                  zV = SUM( pa_i(ji,1:i_fill-1) * ph_i(ji,1:i_fill-1) )
-                  ph_i(ji,i_fill) = ( phti(ji) * pati(ji) - zV ) / MAX( pa_i(ji,i_fill), epsi10 ) 
-                  !
-                  ! correction if concentration of upper cat is greater than lower cat
-                  !    (it should be a gaussian around jl0 but sometimes it is not)
-                  IF ( jl0 /= jpl ) THEN
-                     DO jl = jpl, jl0+1, -1
-                        IF ( pa_i(ji,jl) > pa_i(ji,jl-1) ) THEN
-                           zdv = ph_i(ji,jl) * pa_i(ji,jl)
-                           ph_i(ji,jl    ) = 0._wp
-                           pa_i (ji,jl    ) = 0._wp
-                           pa_i (ji,1:jl-1) = pa_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * phti(ji), epsi10 )
-                        END IF
-                     END DO
-                  ENDIF
-                  !
-               ENDIF
-               !
-               ! Compatibility tests
-               zconv = ABS( pati(ji) - SUM( pa_i(ji,1:jpl) ) ) 
-               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation
-               !
-               zconv = ABS( phti(ji)*pati(ji) - SUM( pa_i(ji,1:jpl)*ph_i(ji,1:jpl) ) )
-               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation
-               !
-               IF ( ph_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                   ! Test 3: thickness of the last category is in-bounds ?
-               !
-               itest(4) = 1
-               DO jl = 1, i_fill
-                  IF ( pa_i(ji,jl) < 0._wp ) itest(4) = 0                                  ! Test 4: positivity of ice concentrations
-               END DO
-               !                                         !----------------------------
-            END DO                                       ! end iteration on categories
-            !                                            !----------------------------
+            ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity
+            pa_i(ji,jpl) = pati(ji) * z1_hti(ji) * ( phti(ji) + 2.*hi_max(jpl-1) ) * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) )
+
+            ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity
+            zv = pati(ji) * z1_hti(ji) * ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jpl-1) + 2.*hi_max(jpl-1)*hi_max(jpl-1) ) &
+               &                         * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) )
+            ! thickness
+            IF( pa_i(ji,jpl) > epsi06 ) THEN
+               ph_i(ji,jpl) = zv / pa_i(ji,jpl)
+            else
+               ph_i(ji,jpl) = 0.
+               pa_i(ji,jpl) = 0.
+            ENDIF
          ENDIF
-      END DO
-
+         !
+      ENDDO
+      !
       ! Add Snow in each category where pa_i is not 0
       DO jl = 1, jpl
          DO ji = 1, idim
             IF( pa_i(ji,jl) > 0._wp ) THEN
-               ph_s(ji,jl) = ph_i(ji,jl) * ( phts(ji) / phti(ji) )
+               ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji)
                ! In case snow load is in excess that would lead to transformation from snow to ice
                ! Then, transfer the snow excess into the ice (different from icethd_dh)
@@ -1009,4 +967,6 @@
          END DO
       END DO
+      !
+      DEALLOCATE( z1_hti )
       !
       ! == temperature and salinity == !
