Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 2382)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 2383)
@@ -232,12 +232,10 @@
 
             !  computing freshwater exchanges at the ice/ocean interface
-            zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj) )  &   !  evaporation over oceanic fraction
-               &   + tprecip(ji,jj) *         at_i(ji,jj)    &   !  total precipitation
-               ! old fashioned way               
-               !              &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )  &   !  remov. snow precip over ice
-               &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   !  remov. snow precip over ice
-               &   - rdmsnif(ji,jj) * r1_rdtice                       &   !  freshwaterflux due to snow melting 
-               ! new contribution from snow falling when ridging
-               &   + fmmec(ji,jj)
+            zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction
+               &   + tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean
+               &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice
+               &   - rdmsnif(ji,jj) * r1_rdtice                       &   ! freshwaterflux due to snow melting 
+               &   + fmmec(ji,jj)                                         ! snow falling when ridging
+
 
             !  computing salt exchanges at the ice/ocean interface
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 2382)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 2383)
@@ -7,6 +7,7 @@
    !!            2.0  ! 2002-07 (C. Ethe, G. Madec)  LIM-2 (F90 rewriting)
    !!            3.0  ! 2005-11 (M. Vancoppenolle)  LIM-3 : Multi-layer thermodynamics + salinity variations
-   !!             -   ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec and lim_thd_con_dif
+   !!             -   ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec, lim_thd_con_dh and lim_thd_con_dif
    !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif
+   !!            3.3  ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas)
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -14,27 +15,27 @@
    !!   'key_lim3'                                      LIM3 sea-ice model
    !!----------------------------------------------------------------------
-   !!   lim_thd      : thermodynamic of sea ice
-   !!   lim_thd_init : initialisation of sea-ice thermodynamic
+   !!   lim_thd        : thermodynamic of sea ice
+   !!   lim_thd_init   : initialisation of sea-ice thermodynamic
    !!----------------------------------------------------------------------
    USE phycst          ! physical constants
    USE dom_oce         ! ocean space and time domain variables
-   USE ice             ! LIM sea-ice variables
-   USE par_ice
+   USE ice             ! LIM: sea-ice variables
+   USE par_ice         ! LIM: sea-ice parameters
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
    USE thd_ice         ! LIM thermodynamic sea-ice variables
    USE dom_ice         ! LIM sea-ice domain
-   USE domvvl
-   USE iceini
-   USE limthd_dif
-   USE limthd_dh
-   USE limthd_sal
-   USE limthd_ent
-   USE limtab
-   USE limvar
+   USE domvvl          ! domain: variable volume level
+   USE iceini          ! LIM: sea-ice initialization
+   USE limthd_dif      ! LIM: thermodynamics, vertical diffusion
+   USE limthd_dh       ! LIM: thermodynamics, ice and snow thickness variation
+   USE limthd_sal      ! LIM: thermodynamics, ice salinity
+   USE limthd_ent      ! LIM: thermodynamics, ice enthalpy redistribution
+   USE limtab          ! LIM: 1D <==> 2D transformation
+   USE limvar          ! LIM: ???
+   USE lbclnk          ! lateral boundary condition - MPP links
+   USE lib_mpp         ! MPP library
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
-   USE lbclnk
-   USE lib_mpp
 
    IMPLICIT NONE
@@ -43,10 +44,10 @@
    PUBLIC   lim_thd    ! called by lim_step
 
-   REAL(wp) ::   epsi20 = 1e-20   ! constant values
-   REAL(wp) ::   epsi16 = 1e-16   !
-   REAL(wp) ::   epsi06 = 1e-06   !
-   REAL(wp) ::   epsi04 = 1e-04   !
-   REAL(wp) ::   zzero  = 0.e0    !
-   REAL(wp) ::   zone   = 1.e0    !
+   REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values
+   REAL(wp) ::   epsi16 = 1e-16_wp   !
+   REAL(wp) ::   epsi06 = 1e-06_wp   !
+   REAL(wp) ::   epsi04 = 1e-04_wp   !
+   REAL(wp) ::   zzero  = 0._wp      !
+   REAL(wp) ::   zone   = 1._wp      !
 
    !! * Substitutions
@@ -56,7 +57,6 @@
    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -86,5 +86,5 @@
       REAL(wp) ::   zfric_umax = 2e-02    ! upper bound for the friction velocity
       REAL(wp) ::   zinda, zindb, zthsnice, zfric_u    ! temporary scalar
-      REAL(wp) ::   zfnsol, zfontn, zfntlat, zpareff   !    -         -
+      REAL(wp) ::   zfntlat, zpareff   !    -         -
       REAL(wp) ::   zeps, zareamin, zcoef
       REAL(wp), DIMENSION(jpi,jpj) ::   zqlbsbq   ! link with lead energy budget qldif
@@ -104,5 +104,4 @@
       !--------------------
       ! Change the units of heat content; from global units to J.m3
-
       DO jl = 1, jpl
          DO jk = 1, nlay_i
@@ -171,5 +170,5 @@
             pfrld(ji,jj)   = 1.0 - at_i(ji,jj)
             zinda          = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )
-
+            !
             !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
             !           !  practically no "direct lateral ablation"
@@ -186,15 +185,13 @@
             qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice
             !                       
-
-            ! still need to be updated : fdtcn !!!!
-            !           !-- Lead heat budget (part 1, next one is in limthd_dh
-            !           !-- qldif -- (or qldif_1d in 1d routines)
-            zfontn         = sprecip(ji,jj) * lfus              ! energy of melting
-            zfnsol         = qns(ji,jj)                         ! total non solar flux
-            qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj)                               &
-               &                               + zfnsol + fdtcn(ji,jj) - zfontn      &
-               &                               + ( 1.0 - zindb ) * fsbbq(ji,jj)  )   &
-               &                               * ( 1.0 - at_i(ji,jj) ) * rdt_ice    
-
+            !           !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) 
+            !           !   caution: exponent betas used as more snow can fallinto leads
+            qldif(ji,jj) =  tms(ji,jj) * rdt_ice  * (                             &
+               &   pfrld(ji,jj)        * (  qsr(ji,jj)                            &   ! solar heat
+               &                            + qns(ji,jj)                          &   ! non solar heat
+               &                            + fdtcn(ji,jj)                        &   ! turbulent ice-ocean heat
+               &                            + fsbbq(ji,jj) * ( 1.0 - zindb )  )   &   ! residual heat from previous step
+               & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus                    )   ! latent heat of sprecip melting
+            !
             ! Positive heat budget is used for bottom ablation
             zfntlat        = 1.0 - MAX( zzero , SIGN( zone ,  - qldif(ji,jj) ) )
@@ -203,12 +200,12 @@
             != 0 if ice and positive heat budget and 1 if one of those two is false
             zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( at_i(ji,jj) * rdt_ice , epsi16 )
-
+            !
             ! Heat budget of the lead, energy transferred from ice to ocean
             qldif  (ji,jj) = zpareff * qldif(ji,jj)
             qdtcn  (ji,jj) = zpareff * qdtcn(ji,jj)
-
+            !
             ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx)
             qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda )
-
+            !
             ! oceanic heat flux (limthd_dh)
             fbif   (ji,jj) = zindb * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) )
@@ -228,5 +225,5 @@
          ENDIF
 
-         zareamin = 1.0e-10
+         zareamin = 1.e-10
          nbpb = 0
          DO jj = 1, jpj
@@ -364,10 +361,10 @@
             CALL tab_1d_2d( nbpb, rdvonif, npb, dvnbq_1d  (1:nbpb), jpi, jpj ) 
             CALL tab_1d_2d( nbpb, fseqv  , npb, fseqv_1d  (1:nbpb), jpi, jpj )
-
+            !
             IF( num_sal == 2 ) THEN
                CALL tab_1d_2d( nbpb, fsbri, npb, fsbri_1d(1:nbpb), jpi, jpj )
                CALL tab_1d_2d( nbpb, fhbri, npb, fhbri_1d(1:nbpb), jpi, jpj )
             ENDIF
-
+            !
             !+++++
             !temporary stuff for a dummy version
@@ -380,5 +377,5 @@
             CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qnsr_ice_1d(1:nbpb), jpi, jpj)
             !+++++
-
+            !
             IF( lk_mpp )   CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
          ENDIF
@@ -393,5 +390,4 @@
       ! 5.1) Ice heat content              
       !------------------------
-
       ! Enthalpies are global variables we have to readjust the units
       zcoef = 1.e0 / ( unit_fac * REAL(nlay_i) )
@@ -406,5 +402,4 @@
       ! 5.2) Snow heat content              
       !------------------------
-
       ! Enthalpies are global variables we have to readjust the units
       zcoef = 1.e0 / ( unit_fac * REAL(nlay_s) )
@@ -429,5 +424,5 @@
       IF( con_i )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)
 
-      IF(ln_ctl) THEN   ! Control print
+      IF(ln_ctl) THEN            ! Control print
          CALL prt_ctl_info(' ')
          CALL prt_ctl_info(' - Cell values : ')
@@ -459,7 +454,6 @@
             END DO
          END DO
-
       ENDIF
-
+      !
    END SUBROUTINE lim_thd
 
@@ -529,17 +523,14 @@
       END DO
       ! layer by layer
-      dq_i_layer(:,:)    = q_i_layer_fin(:,:) - q_i_layer_in(:,:)
+      dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:)
 
       !----------------------------------------
       ! Atmospheric heat flux, ice heat budget
       !----------------------------------------
-
-      DO ji = kideb, kiut
-         zji = MOD( npb(ji) - 1, jpi ) + 1
-         zjj = ( npb(ji) - 1 ) / jpi + 1
-
-         fatm(ji,jl) = qnsr_ice_1d(ji) + (1.0-i0(ji))*qsr_ice_1d(ji)
-
-         sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji)*i0(ji) - fstroc(zji,zjj,jl)
+      DO ji = kideb, kiut
+         zji = MOD( npb(ji) - 1 , jpi ) + 1
+         zjj =    ( npb(ji) - 1 ) / jpi + 1
+         fatm     (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji)
+         sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(zji,zjj,jl)
       END DO
 
@@ -547,10 +538,9 @@
       ! Conservation error
       !--------------------
-
       DO ji = kideb, kiut
          cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
       END DO
 
-      numce = 0
+      numce  = 0
       meance = 0.0
       DO ji = kideb, kiut
@@ -559,6 +549,6 @@
             meance = meance + cons_error(ji,jl)
          ENDIF
-      ENDDO
-      IF (numce .GT. 0 ) meance = meance / numce
+      END DO
+      IF( numce .GT. 0 ) meance = meance / numce
 
       WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err
@@ -570,16 +560,15 @@
       ! Surface error due to imbalance between Fatm and Fcsu
       !-------------------------------------------------------
-      numce  = 0.0
+      numce  = 0
       meance = 0.0
 
       DO ji = kideb, kiut
          surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) )
-         IF ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. &
-            max_surf_err ) ) THEN
+         IF( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) THEN
             numce = numce + 1 
             meance = meance + surf_error(ji,jl)
          ENDIF
       ENDDO
-      IF (numce .GT. 0 ) meance = meance / numce
+      IF( numce .GT. 0 ) meance = meance / numce
 
       WRITE(numout,*) ' Maximum tolerated surface error : ', max_surf_err
@@ -595,5 +584,4 @@
       ! Write ice state in case of big errors
       !---------------------------------------
-
       DO ji = kideb, kiut
          IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. &
@@ -601,5 +589,5 @@
             zji                 = MOD( npb(ji) - 1, jpi ) + 1
             zjj                 = ( npb(ji) - 1 ) / jpi + 1
-
+            !
             WRITE(numout,*) ' alerte 1     '
             WRITE(numout,*) ' Untolerated conservation / surface error after '
@@ -617,6 +605,5 @@
             !        WRITE(numout,*) ' qt_i_fin  : ', qt_i_fin(ji,jl)
             !        WRITE(numout,*) ' qt_s_fin  : ', qt_s_fin(ji,jl)
-            !        WRITE(numout,*) ' qt        : ', qt_i_fin(ji,jl) + &
-            !                                         qt_s_fin(ji,jl)
+            !        WRITE(numout,*) ' qt        : ', qt_i_fin(ji,jl) + qt_s_fin(ji,jl)
             WRITE(numout,*) ' ht_i       : ', ht_i_b(ji)
             WRITE(numout,*) ' ht_s       : ', ht_s_b(ji)
@@ -645,17 +632,12 @@
             WRITE(numout,*)
             WRITE(numout,*) ' Layer by layer ... '
-            WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - &
-               qt_s_in(ji,jl) )  & 
-               / rdt_ice
-            WRITE(numout,*) ' dfc_snow  : ', fc_s(ji,1) -      &
-               fc_s(ji,0)
+            WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice
+            WRITE(numout,*) ' dfc_snow  : ', fc_s(ji,1) - fc_s(ji,0)
             DO jk = 1, nlay_i
                WRITE(numout,*) ' layer  : ', jk
                WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice  
                WRITE(numout,*) ' radab  : ', radab(ji,jk)
-               WRITE(numout,*) ' dfc_i  : ', fc_i(ji,jk) -               &
-                  fc_i(ji,jk-1)
-               WRITE(numout,*) ' tot f  : ', fc_i(ji,jk) -               &
-                  fc_i(ji,jk-1) - radab(ji,jk)
+               WRITE(numout,*) ' dfc_i  : ', fc_i(ji,jk) - fc_i(ji,jk-1)
+               WRITE(numout,*) ' tot f  : ', fc_i(ji,jk) - fc_i(ji,jk-1) - radab(ji,jk)
             END DO
 
@@ -667,12 +649,9 @@
 
 
-   SUBROUTINE lim_thd_con_dh(kideb,kiut,jl)
+   SUBROUTINE lim_thd_con_dh( kideb, kiut, jl )
       !!-----------------------------------------------------------------------
       !!                   ***  ROUTINE lim_thd_con_dh  *** 
       !!                 
       !! ** Purpose :   Test energy conservation after enthalpy redistr.
-      !!
-      !! history :
-      !!  9.9  ! 07-04 (M.Vancoppenolle) original code
       !!-----------------------------------------------------------------------
       INTEGER, INTENT(in) ::        &
@@ -750,10 +729,9 @@
       ! Write ice state in case of big errors
       !---------------------------------------
-
       DO ji = kideb, kiut
          IF ( cons_error(ji,jl) .GT. max_cons_err  ) THEN
             zji                 = MOD( npb(ji) - 1, jpi ) + 1
             zjj                 = ( npb(ji) - 1 ) / jpi + 1
-
+            !
             WRITE(numout,*) ' alerte 1 - category : ', jl
             WRITE(numout,*) ' Untolerated conservation error after limthd_ent '
@@ -776,10 +754,8 @@
             WRITE(numout,*) ' qt_s_in    : ', qt_s_in(ji,jl) / rdt_ice
             WRITE(numout,*) ' qt_i_in    : ', qt_i_in(ji,jl) / rdt_ice
-            WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + &
-               qt_s_in(ji,jl) ) / rdt_ice
+            WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) / rdt_ice
             WRITE(numout,*) ' qt_s_fin   : ', qt_s_fin(ji,jl) / rdt_ice
             WRITE(numout,*) ' qt_i_fin   : ', qt_i_fin(ji,jl) / rdt_ice
-            WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + &
-               qt_s_fin(ji,jl) ) / rdt_ice
+            WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) / rdt_ice
             WRITE(numout,*) ' * '
             WRITE(numout,*) ' Ice variables --- : '
@@ -790,5 +766,4 @@
             WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji)
             WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
-
          ENDIF
          !
@@ -805,5 +780,4 @@
       !!
       !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
-      !!
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
@@ -832,5 +806,4 @@
 
    SUBROUTINE lim_thd_init
-
       !!-----------------------------------------------------------------------
       !!                   ***  ROUTINE lim_thd_init *** 
@@ -843,5 +816,5 @@
       !!
       !! ** input   :   Namelist namicether
-       !!-------------------------------------------------------------------
+      !!-------------------------------------------------------------------
       NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb,   &
          &                hicmin, hiclim, amax  ,                                &
