Index: NEMO/branches/2020/dev_12905_xios_ancil/src/ICE/icevar.F90
===================================================================
--- NEMO/branches/2020/dev_12905_xios_ancil/src/ICE/icevar.F90	(revision 13003)
+++ NEMO/branches/2020/dev_12905_xios_ancil/src/ICE/icevar.F90	(revision 13766)
@@ -51,4 +51,6 @@
    !!   ice_var_sshdyn    : compute equivalent ssh in lead
    !!   ice_var_itd       : convert N-cat to M-cat
+   !!   ice_var_snwfra    : fraction of ice covered by snow
+   !!   ice_var_snwblow   : distribute snow fall between ice and ocean
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
@@ -77,4 +79,6 @@
    PUBLIC   ice_var_sshdyn
    PUBLIC   ice_var_itd
+   PUBLIC   ice_var_snwfra
+   PUBLIC   ice_var_snwblow
 
    INTERFACE ice_var_itd
@@ -84,4 +88,13 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+
+   INTERFACE ice_var_snwfra
+      MODULE PROCEDURE ice_var_snwfra_1d, ice_var_snwfra_2d, ice_var_snwfra_3d
+   END INTERFACE
+
+   INTERFACE ice_var_snwblow
+      MODULE PROCEDURE ice_var_snwblow_1d, ice_var_snwblow_2d
+   END INTERFACE
+
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -115,4 +128,5 @@
       at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
       vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
+      vt_il(:,:) = SUM( v_il(:,:,:), dim=3 )
       !
       ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction  
@@ -166,6 +180,6 @@
          !
          !                           ! mean melt pond depth
-         WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:)
-         ELSEWHERE                      ;   hm_ip(:,:) = 0._wp
+         WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:)   ;   hm_il(:,:) = vt_il(:,:) / at_ip(:,:)
+         ELSEWHERE                      ;   hm_ip(:,:) = 0._wp                     ;   hm_il(:,:) = 0._wp
          END WHERE         
          !
@@ -191,5 +205,7 @@
       REAL(wp) ::   zhmax, z1_zhmax                 !   -      -
       REAL(wp) ::   zlay_i, zlay_s                  !   -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
+      REAL(wp), PARAMETER ::   zhl_max =  0.015_wp  ! pond lid thickness above which the ponds disappear from the albedo calculation
+      REAL(wp), PARAMETER ::   zhl_min =  0.005_wp  ! pond lid thickness below which the full pond area is used in the albedo calculation
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i, z1_a_ip, za_s_fra
       !!-------------------------------------------------------------------
 
@@ -210,4 +226,8 @@
       ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
       END WHERE
+      !
+      WHERE( a_ip(:,:,:) > epsi20 )  ;   z1_a_ip(:,:,:) = 1._wp / a_ip(:,:,:)
+      ELSEWHERE                      ;   z1_a_ip(:,:,:) = 0._wp
+      END WHERE
       !                                           !--- ice thickness
       h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
@@ -224,9 +244,17 @@
       !                                           !--- ice age      
       o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
-      !                                           !--- pond fraction and thickness      
+      !                                           !--- pond and lid thickness      
+      h_ip(:,:,:) = v_ip(:,:,:) * z1_a_ip(:,:,:)
+      h_il(:,:,:) = v_il(:,:,:) * z1_a_ip(:,:,:)
+      !                                           !--- melt pond effective area (used for albedo)
       a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
-      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
-      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp
-      END WHERE
+      WHERE    ( h_il(:,:,:) <= zhl_min )  ;   a_ip_eff(:,:,:) = a_ip_frac(:,:,:)       ! lid is very thin.  Expose all the pond
+      ELSEWHERE( h_il(:,:,:) >= zhl_max )  ;   a_ip_eff(:,:,:) = 0._wp                  ! lid is very thick. Cover all the pond up with ice and snow
+      ELSEWHERE                            ;   a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * &   ! lid is in between. Expose part of the pond
+         &                                                       ( h_il(:,:,:) - zhl_min ) / ( zhl_max - zhl_min )
+      END WHERE
+      !
+      CALL ice_var_snwfra( h_s, za_s_fra )           ! calculate ice fraction covered by snow
+      a_ip_eff = MIN( a_ip_eff, 1._wp - za_s_fra )   ! make sure (a_ip_eff + a_s_fra) <= 1
       !
       !                                           !---  salinity (with a minimum value imposed everywhere)     
@@ -243,5 +271,5 @@
       zlay_i   = REAL( nlay_i , wp )    ! number of layers
       DO jl = 1, jpl
-         DO_3D_11_11( 1, nlay_i )
+         DO_3D( 1, 1, 1, 1, 1, nlay_i )
             IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area 
                !
@@ -292,4 +320,5 @@
       sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
       v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
+      v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:)
       !
    END SUBROUTINE ice_var_eqv2glo
@@ -347,5 +376,5 @@
          z1_dS = 1._wp / ( zsi1 - zsi0 )
          DO jl = 1, jpl
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
                !                             ! force a constant profile when SSS too low (Baltic Sea)
@@ -356,5 +385,5 @@
          ! Computation of the profile
          DO jl = 1, jpl
-            DO_3D_11_11( 1, nlay_i )
+            DO_3D( 1, 1, 1, 1, 1, nlay_i )
                !                          ! linear profile with 0 surface value
                zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
@@ -486,5 +515,5 @@
          ! Zap ice energy and use ocean heat to melt ice
          !-----------------------------------------------------------------
-         DO_3D_11_11( 1, nlay_i )
+         DO_3D( 1, 1, 1, 1, 1, nlay_i )
             ! update exchanges with ocean
             hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0
@@ -493,5 +522,5 @@
          END_3D
          !
-         DO_3D_11_11( 1, nlay_s )
+         DO_3D( 1, 1, 1, 1, 1, nlay_s )
             ! update exchanges with ocean
             hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0
@@ -503,5 +532,5 @@
          ! zap ice and snow volume, add water and salt to ocean
          !-----------------------------------------------------------------
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ! update exchanges with ocean
             sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_Dt_ice
@@ -521,4 +550,5 @@
             a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
             v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
+            v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj)
             !
          END_2D
@@ -542,5 +572,5 @@
 
 
-   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_var_zapneg ***
@@ -557,4 +587,5 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -574,5 +605,5 @@
          ! zap ice energy and send it to the ocean
          !----------------------------------------
-         DO_3D_11_11( 1, nlay_i )
+         DO_3D( 1, 1, 1, 1, 1, nlay_i )
             IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
                hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
@@ -581,5 +612,5 @@
          END_3D
          !
-         DO_3D_11_11( 1, nlay_s )
+         DO_3D( 1, 1, 1, 1, 1, nlay_s )
             IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
                hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
@@ -591,5 +622,5 @@
          ! zap ice and snow volume, add water and salt to ocean
          !-----------------------------------------------------
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
                wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
@@ -613,10 +644,10 @@
       WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp
       WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
-      !                                                        but it does not change conservation, so keep it this way is ok
+      WHERE( pv_il (:,:,:) < 0._wp )   pv_il (:,:,:) = 0._wp !    but it does not change conservation, so keep it this way is ok
       !
    END SUBROUTINE ice_var_zapneg
 
 
-   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
+   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_var_roundoff ***
@@ -631,18 +662,23 @@
       REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
       !!-------------------------------------------------------------------
       !
-      WHERE( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
-      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
-      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
-      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
-      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
-      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
-      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
-      IF( ln_pnd_H12 ) THEN
-         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
-         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
+
+      WHERE( pa_i (1:npti,:)   < 0._wp )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
+      WHERE( pv_i (1:npti,:)   < 0._wp )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
+      WHERE( pv_s (1:npti,:)   < 0._wp )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
+      WHERE( psv_i(1:npti,:)   < 0._wp )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
+      WHERE( poa_i(1:npti,:)   < 0._wp )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
+      WHERE( pe_i (1:npti,:,:) < 0._wp )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
+      WHERE( pe_s (1:npti,:,:) < 0._wp )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
+      IF( ln_pnd_LEV ) THEN
+         WHERE( pa_ip(1:npti,:) < 0._wp )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
+         WHERE( pv_ip(1:npti,:) < 0._wp )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
+         IF( ln_pnd_lids ) THEN
+            WHERE( pv_il(1:npti,:) < 0._wp .AND. pv_il(1:npti,:) > -epsi10 ) pv_il(1:npti,:)   = 0._wp   ! v_il must be >= 0
+         ENDIF
       ENDIF
       !
@@ -763,6 +799,6 @@
    !! ** Purpose :  converting N-cat ice to jpl ice categories
    !!-------------------------------------------------------------------
-   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
+   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !! ** Purpose :  converting 1-cat ice to 1 ice category
@@ -770,6 +806,6 @@
       REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
-      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !!-------------------------------------------------------------------
       ! == thickness and concentration == !
@@ -785,9 +821,10 @@
       pa_ip(:) = patip(:)
       ph_ip(:) = phtip(:)
+      ph_il(:) = phtil(:)
       
    END SUBROUTINE ice_var_itd_1c1c
 
-   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
+   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !! ** Purpose :  converting N-cat ice to 1 ice category
@@ -795,6 +832,6 @@
       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:)  , INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
-      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs
@@ -831,6 +868,10 @@
       ! == ponds == !
       pa_ip(:) = SUM( patip(:,:), dim=2 )
-      WHERE( pa_ip(:) /= 0._wp )   ;   ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
-      ELSEWHERE                    ;   ph_ip(:) = 0._wp
+      WHERE( pa_ip(:) /= 0._wp )
+         ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
+         ph_il(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / pa_ip(:)
+      ELSEWHERE
+         ph_ip(:) = 0._wp
+         ph_il(:) = 0._wp
       END WHERE
       !
@@ -839,6 +880,6 @@
    END SUBROUTINE ice_var_itd_Nc1c
    
-   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
+   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !!
@@ -862,6 +903,6 @@
       REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra, z1_hti
@@ -953,5 +994,4 @@
          pt_su(:,jl) = ptmsu(:)
          ps_i (:,jl) = psmi (:)
-         ps_i (:,jl) = psmi (:)         
       END DO
       !
@@ -974,10 +1014,19 @@
          END WHERE
       END DO
+      ! keep the same v_il/v_i ratio for each category
+      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtil(:) * patip(:) ) / ( phti(:) * pati(:) )
+      ELSEWHERE                                 ;   zfra(:) = 0._wp
+      END WHERE
+      DO jl = 1, jpl
+         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
+         ELSEWHERE                       ;   ph_il(:,jl) = 0._wp
+         END WHERE
+      END DO
       DEALLOCATE( zfra )
       !
    END SUBROUTINE ice_var_itd_1cMc
 
-   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, &
-      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip )
+   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
+      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
       !!-------------------------------------------------------------------
       !!
@@ -994,5 +1043,5 @@
       !!
       !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
-       !!                   by removing 25% ice area from jlmin and jlmax (resp.) 
+      !!                   by removing 25% ice area from jlmin and jlmax (resp.) 
       !!              
       !!               3) Expand the filling to the empty cat between jlmin and jlmax 
@@ -1010,6 +1059,6 @@
       REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
-      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
       !
       INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2
@@ -1040,4 +1089,5 @@
          pa_ip(:,:) = patip(:,:)
          ph_ip(:,:) = phtip(:,:)
+         ph_il(:,:) = phtil(:,:)
          !                              ! ---------------------- !
       ELSEIF( icat == 1 ) THEN          ! input cat = 1          !
@@ -1045,6 +1095,6 @@
          CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), &
             &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), &
-            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), &
-            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:)  )
+            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), phtil(:,1), &
+            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ph_il(:,:)  )
          !                              ! ---------------------- !
       ELSEIF( jpl == 1 ) THEN           ! output cat = 1         !
@@ -1052,6 +1102,6 @@
          CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), &
             &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), &
-            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), &
-            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1)  )
+            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), phtil(:,:), &
+            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ph_il(:,1)  )
          !                              ! ----------------------- !
       ELSE                              ! input cat /= output cat !
@@ -1195,4 +1245,15 @@
             END WHERE
          END DO
+         ! keep the same v_il/v_i ratio for each category
+         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
+            zfra(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
+         ELSEWHERE
+            zfra(:) = 0._wp
+         END WHERE
+         DO jl = 1, jpl
+            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
+            ELSEWHERE                       ;   ph_il(:,jl) = 0._wp
+            END WHERE
+         END DO
          DEALLOCATE( zfra )
          !
@@ -1200,4 +1261,88 @@
       !
    END SUBROUTINE ice_var_itd_NcMc
+
+   !!-------------------------------------------------------------------
+   !! INTERFACE ice_var_snwfra
+   !!
+   !! ** Purpose :  fraction of ice covered by snow
+   !!
+   !! ** Method  :  In absence of proper snow model on top of sea ice,
+   !!               we argue that snow does not cover the whole ice because
+   !!               of wind blowing...
+   !!                
+   !! ** Arguments : ph_s: snow thickness
+   !!                
+   !! ** Output    : pa_s_fra: fraction of ice covered by snow
+   !!
+   !!-------------------------------------------------------------------
+   SUBROUTINE ice_var_snwfra_3d( ph_s, pa_s_fra )
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ph_s        ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
+      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover
+         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
+         ELSEWHERE             ; pa_s_fra = 0._wp
+         END WHERE
+      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style)
+         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
+      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style)
+         pa_s_fra = ph_s / ( ph_s + 0.02_wp )
+      ENDIF
+   END SUBROUTINE ice_var_snwfra_3d
+
+   SUBROUTINE ice_var_snwfra_2d( ph_s, pa_s_fra )
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ph_s        ! snow thickness
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
+      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover
+         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
+         ELSEWHERE             ; pa_s_fra = 0._wp
+         END WHERE
+      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style)
+         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
+      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style)
+         pa_s_fra = ph_s / ( ph_s + 0.02_wp )
+      ENDIF
+   END SUBROUTINE ice_var_snwfra_2d
+
+   SUBROUTINE ice_var_snwfra_1d( ph_s, pa_s_fra )
+      REAL(wp), DIMENSION(:), INTENT(in   ) ::   ph_s        ! snow thickness
+      REAL(wp), DIMENSION(:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
+      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover
+         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
+         ELSEWHERE             ; pa_s_fra = 0._wp
+         END WHERE
+      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style)
+         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
+      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style)
+         pa_s_fra = ph_s / ( ph_s + 0.02_wp )
+      ENDIF
+   END SUBROUTINE ice_var_snwfra_1d
+   
+   !!--------------------------------------------------------------------------
+   !! INTERFACE ice_var_snwblow
+   !!
+   !! ** Purpose :   Compute distribution of precip over the ice
+   !!
+   !!                Snow accumulation in one thermodynamic time step
+   !!                snowfall is partitionned between leads and ice.
+   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads
+   !!                but because of the winds, more snow falls on leads than on sea ice
+   !!                and a greater fraction (1-at_i)^beta of the total mass of snow 
+   !!                (beta < 1) falls in leads.
+   !!                In reality, beta depends on wind speed, 
+   !!                and should decrease with increasing wind speed but here, it is 
+   !!                considered as a constant. an average value is 0.66
+   !!--------------------------------------------------------------------------
+!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
+   SUBROUTINE ice_var_snwblow_2d( pin, pout )
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
+      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
+      pout = ( 1._wp - ( pin )**rn_snwblow )
+   END SUBROUTINE ice_var_snwblow_2d
+
+   SUBROUTINE ice_var_snwblow_1d( pin, pout )
+      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
+      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
+      pout = ( 1._wp - ( pin )**rn_snwblow )
+   END SUBROUTINE ice_var_snwblow_1d
 
 #else
