Index: /trunk/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- /trunk/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 7812)
+++ /trunk/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 7813)
@@ -125,4 +125,6 @@
                                     !     3: activate G(he) only                 --- temporary option
                                     !     4: activate extra lateral melting only --- temporary option
+   rn_cdsn     = 0.31              !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971)
+                                   !  Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
                   ! -- limthd_dh -- !
    ln_limdH       = .true.          !  activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
Index: /trunk/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /trunk/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7812)
+++ /trunk/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7813)
@@ -438,11 +438,12 @@
 &namsbc_alb    !   albedo parameters
 !-----------------------------------------------------------------------
-   nn_ice_alb  =    1      !  parameterization of ice/snow albedo
-                           !     0: Shine & Henderson-Sellers (JGR 1985)
-                           !     1: "home made" based on Brandt et al. (J. Climate 2005)
-                           !                         and Grenfell & Perovich (JGR 2004)
-   rn_albice   =  0.50     !  albedo of bare puddled ice (values from 0.49 to 0.58)
-                           !     0.53 (default) => if nn_ice_alb=0
-                           !     0.50 (default) => if nn_ice_alb=1
+   nn_ice_alb   =    1   !  parameterization of ice/snow albedo
+                         !     0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo
+                         !     1: "home made" based on Brandt et al. (JClim 2005) and Grenfell & Perovich (JGR 2004),
+                         !        giving cloud-sky albedo
+   rn_alb_sdry  =  0.85  !  dry snow albedo         : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky)
+   rn_alb_smlt  =  0.75  !  melting snow albedo     : 0.65 ( '' )          ; 0.75 ( '' )          ; obs 0.72-0.82 ( '' )
+   rn_alb_idry  =  0.60  !  dry ice albedo          : 0.72 ( '' )          ; 0.60 ( '' )          ; obs 0.54-0.65 ( '' )
+   rn_alb_imlt  =  0.50  !  bare puddled ice albedo : 0.53 ( '' )          ; 0.50 ( '' )          ; obs 0.49-0.58 ( '' )
 /
 !-----------------------------------------------------------------------
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 7812)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 7813)
@@ -232,4 +232,5 @@
    LOGICAL , PUBLIC ::   ln_it_qnsice     !: iterate surface flux with changing surface temperature or not (F)
    INTEGER , PUBLIC ::   nn_monocat       !: virtual ITD mono-category parameterizations (1) or not (0)
+   REAL(wp), PUBLIC ::   rn_cdsn          !: thermal conductivity of the snow [W/m/K]
                                           ! -- limthd_dh -- !
    LOGICAL , PUBLIC ::   ln_limdH         !: activate ice thickness change from growing/melting (T) or not (F)
@@ -458,5 +459,5 @@
    FUNCTION ice_alloc()
       !!-----------------------------------------------------------------
-      !!               *** Routine ice_alloc_2 ***
+      !!               *** Routine ice_alloc ***
       !!-----------------------------------------------------------------
       INTEGER :: ice_alloc
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 7812)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 7813)
@@ -147,15 +147,7 @@
          z2d(:,:) = t_su(:,:,jl)
          CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-      END DO
-
-      DO jl = 1, jpl 
-         WRITE(zchar,'(I2.2)') jl
          znam = 'tempt_sl1'//'_htc'//zchar
          z2d(:,:) = e_s(:,:,1,jl)
          CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-      END DO
-
-      DO jl = 1, jpl 
-         WRITE(zchar,'(I2.2)') jl
          DO jk = 1, nlay_i 
             WRITE(zchar1,'(I2.2)') jk
@@ -376,15 +368,7 @@
          CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
          t_su(:,:,jl) = z2d(:,:)
-      END DO
-
-      DO jl = 1, jpl 
-         WRITE(zchar,'(I2.2)') jl
          znam = 'tempt_sl1'//'_htc'//zchar
          CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
          e_s(:,:,1,jl) = z2d(:,:)
-      END DO
-
-      DO jl = 1, jpl 
-         WRITE(zchar,'(I2.2)') jl
          DO jk = 1, nlay_i 
             WRITE(zchar1,'(I2.2)') jk
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 7812)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 7813)
@@ -571,7 +571,7 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicethd/ rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon,ln_it_qnsice,nn_monocat,  &
-         &                ln_limdH, rn_betas,                                                          &
-         &                ln_limdA, rn_beta, rn_dmin,                                                  &
+      NAMELIST/namicethd/ rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, ln_it_qnsice, nn_monocat, rn_cdsn,  &
+         &                ln_limdH, rn_betas,                                                                     &
+         &                ln_limdA, rn_beta, rn_dmin,                                                             &
          &                ln_limdO, rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, rn_himin
       !!-------------------------------------------------------------------
@@ -585,4 +585,10 @@
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp )
       IF(lwm) WRITE ( numoni, namicethd )
+      !
+      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN 
+         nn_monocat = 0
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
+      ENDIF
       !
       IF(lwp) THEN                          ! control print
@@ -596,4 +602,5 @@
          WRITE(numout,*)'      iterate the surface non-solar flux (T) or not (F)       ln_it_qnsice = ', ln_it_qnsice
          WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat
+         WRITE(numout,*)'      thermal conductivity of the snow                        rn_cdsn      = ', rn_cdsn
          WRITE(numout,*)'   -- limthd_dh --'
          WRITE(numout,*)'      activate ice thick change from top/bot (T) or not (F)   ln_limdH     = ', ln_limdH
@@ -614,9 +621,4 @@
          WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i
       ENDIF
-      IF( jpl > 1 .AND. nn_monocat == 1 ) THEN 
-         nn_monocat = 0
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
-      ENDIF
       !
    END SUBROUTINE lim_thd_init
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 7812)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 7813)
@@ -376,6 +376,6 @@
 
                ! Effective thickness he (zhe)
-               zfac     = 1._wp / ( rcdsn + zkimean )
-               zratio_s = rcdsn   * zfac
+               zfac     = 1._wp / ( rn_cdsn + zkimean )
+               zratio_s = rn_cdsn   * zfac
                zratio_i = zkimean * zfac
                zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji)
@@ -400,11 +400,11 @@
          DO ji = kideb, kiut
             zfac                  =  1. / MAX( epsi10 , zh_s(ji) )
-            zkappa_s(ji,0)        = zghe(ji) * rcdsn * zfac
-            zkappa_s(ji,nlay_s)   = zghe(ji) * rcdsn * zfac
+            zkappa_s(ji,0)        = zghe(ji) * rn_cdsn * zfac
+            zkappa_s(ji,nlay_s)   = zghe(ji) * rn_cdsn * zfac
          END DO
 
          DO jk = 1, nlay_s-1
             DO ji = kideb , kiut
-               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0 * zh_s(ji) )
+               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rn_cdsn / MAX( epsi10, 2.0 * zh_s(ji) )
             END DO
          END DO
@@ -422,6 +422,6 @@
             zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac
             zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac
-            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / & 
-           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) )
+            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rn_cdsn * ztcond_i(ji,0) / & 
+           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cdsn * zh_i(ji) ) )
             zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
          END DO
Index: /trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 7812)
+++ /trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 7813)
@@ -677,5 +677,5 @@
                zht_i(ji,1:jpl) = 0._wp
                za_i (ji,1:jpl) = 0._wp
-               itest(:)        = 0           
+               itest(:)        = 0      
                
                ! *** case very thin ice: fill only category 1
@@ -722,5 +722,5 @@
                
                ENDIF ! case ice is thick or thin
-               
+            
                !---------------------
                ! Compatibility tests
@@ -747,5 +747,5 @@
          ENDIF ! if zhti > 0
       END DO ! i loop
-      
+
       ! ------------------------------------------------
       ! Adding Snow in each category where za_i is not 0
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7812)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7813)
@@ -64,7 +64,6 @@
 #if defined key_lim3 || defined key_cice
    REAL(wp), PUBLIC ::   rhoic    =  917._wp         !: volumic mass of sea ice                               [kg/m3]
-   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice
-   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow
-   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice 
+   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice                     [W/m/K]
+   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat of fresh ice                            [J/kg/K]
    REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg]
    REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg]
@@ -82,4 +81,7 @@
    REAL(wp), PUBLIC ::   xlic     =  300.33e+6_wp    !: volumetric latent heat fusion of ice                  [J/m3]
    REAL(wp), PUBLIC ::   xsn      =    2.8e+6_wp     !: volumetric latent heat of sublimation of snow         [J/m3]
+#endif
+#if defined key_cice
+   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow                          [W/m/K], now namelist parameter for LIM3
 #endif
 #if defined key_lim3
@@ -157,6 +159,8 @@
       IF(lwp) THEN
          WRITE(numout,*)
+#if defined key_cice
          WRITE(numout,*) '          thermal conductivity of the snow          = ', rcdsn   , ' J/s/m/K'
-         WRITE(numout,*) '          thermal conductivity of the ice           = ', rcdic   , ' J/s/m/K'
+#endif
+         WRITE(numout,*) '          thermal conductivity of pure ice          = ', rcdic   , ' J/s/m/K'
          WRITE(numout,*) '          fresh ice specific heat                   = ', cpic    , ' J/kg/K'
          WRITE(numout,*) '          latent heat of fusion of fresh ice / snow = ', lfus    , ' J/kg'
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 7812)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 7813)
@@ -39,5 +39,5 @@
    !                             !!* namelist namsbc_alb
    INTEGER  ::   nn_ice_alb
-   REAL(wp) ::   rn_albice
+   REAL(wp) ::   rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt
 
    !!----------------------------------------------------------------------
@@ -89,5 +89,4 @@
       INTEGER  ::   ji, jj, jl         ! dummy loop indices
       INTEGER  ::   ijpl               ! number of ice categories (3rd dim of ice input arrays)
-      REAL(wp)            ::   ralb_im, ralb_sf, ralb_sm, ralb_if
       REAL(wp)            ::   zswitch, z1_c1, z1_c2
       REAL(wp)                            ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
@@ -100,5 +99,4 @@
 
       IF( albd_init == 0 )   CALL albedo_init      ! initialization 
-
       
       SELECT CASE ( nn_ice_alb )
@@ -109,12 +107,7 @@
       CASE( 0 )
        
-         ralb_sf = 0.80       ! dry snow
-         ralb_sm = 0.65       ! melting snow
-         ralb_if = 0.72       ! bare frozen ice
-         ralb_im = rn_albice  ! bare puddled ice 
-         
          !  Computation of ice albedo (free of snow)
-         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
-         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = rn_alb_imlt
+         ELSE WHERE                                              ;   zalb(:,:,:) = rn_alb_idry
          END  WHERE
       
@@ -132,15 +125,14 @@
                   ! freezing snow
                   ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
-                  !                                        !  freezing snow        
                   zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
                   zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
-                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   &
-                     &        +         zswitch   * ralb_sf  
+                     &                           + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) / c1  )   &
+                     &        +         zswitch   * rn_alb_sdry  
 
                   ! melting snow
                   ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2
                   zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
-                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   &
-                      &     +         zswitch   *   ralb_sm 
+                  zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) / c2 )   &
+                      &     +         zswitch   *   rn_alb_smlt 
                   !
                   ! snow albedo
@@ -163,23 +155,22 @@
       CASE( 1 ) 
 
-         ralb_im = rn_albice  ! bare puddled ice
 ! compilation of values from literature
-         ralb_sf = 0.85      ! dry snow
-         ralb_sm = 0.75      ! melting snow
-         ralb_if = 0.60      ! bare frozen ice
+!        rn_alb_sdry = 0.85      ! dry snow
+!        rn_alb_smlt = 0.75      ! melting snow
+!        rn_alb_idry = 0.60      ! bare frozen ice
 ! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
-!         ralb_sf = 0.85       ! dry snow
-!         ralb_sm = 0.72       ! melting snow
-!         ralb_if = 0.65       ! bare frozen ice
+!        rn_alb_sdry = 0.85      ! dry snow
+!        rn_alb_smlt = 0.72      ! melting snow
+!        rn_alb_idry = 0.65      ! bare frozen ice
 ! Brandt et al 2005 (East Antarctica)
-!         ralb_sf = 0.87      ! dry snow
-!         ralb_sm = 0.82      ! melting snow
-!         ralb_if = 0.54      ! bare frozen ice
+!        rn_alb_sdry = 0.87      ! dry snow
+!        rn_alb_smlt = 0.82      ! melting snow
+!        rn_alb_idry = 0.54      ! bare frozen ice
 ! 
          !  Computation of ice albedo (free of snow)
          z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
          z1_c2 = 1. / 0.05
-         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im
-         ELSE WHERE                                              ;   zalb = ralb_if
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = rn_alb_imlt
+         ELSE WHERE                                              ;   zalb = rn_alb_idry
          END  WHERE
          
@@ -196,8 +187,8 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zalb_sf = ralb_sf - ( ralb_sf - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
-                  zalb_sm = ralb_sm - ( ralb_sm - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
-
-                   ! snow albedo
+                  zalb_sf = rn_alb_sdry - ( rn_alb_sdry - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
+                  zalb_sm = rn_alb_smlt - ( rn_alb_smlt - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
+
+                  ! snow albedo
                   zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
                   zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
@@ -248,5 +239,5 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
+      NAMELIST/namsbc_alb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry , rn_alb_imlt
       !!----------------------------------------------------------------------
       !
@@ -267,6 +258,9 @@
          WRITE(numout,*) '~~~~~~~'
          WRITE(numout,*) '   Namelist namsbc_alb : albedo '
-         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb
-         WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice
+         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb  = ', nn_ice_alb
+         WRITE(numout,*) '      albedo of dry snow                                  rn_alb_sdry = ', rn_alb_sdry
+         WRITE(numout,*) '      albedo of melting snow                              rn_alb_smlt = ', rn_alb_smlt
+         WRITE(numout,*) '      albedo of dry ice                                   rn_alb_idry = ', rn_alb_idry
+         WRITE(numout,*) '      albedo of bare puddled ice                          rn_alb_imlt = ', rn_alb_imlt
       ENDIF
       !
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7812)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7813)
@@ -1803,21 +1803,27 @@
       WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:)
 
-      ! --- heat flux associated with emp (W/m2) --- !
-      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn(:,:)   &       ! evap
-         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &       ! liquid precip
-         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus )  ! solid precip over ocean + snow melting
-!      zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
-!         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
-      zqemp_ice(:,:) =      zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice (only)
-                                                                                                       ! qevap_ice=0 since we consider Tice=0degC
-      
+      ! Heat content per unit mass of snow (J/kg)
+      WHERE( SUM( a_i, dim=3 ) > 1.e-10 )   ;   zcptsnw(:,:) = cpic * SUM( (tn_ice - rt0) * a_i, dim=3 ) / SUM( a_i, dim=3 )
+      ELSEWHERE                             ;   zcptsnw(:,:) = zcptn(:,:)
+      ENDWHERE
+      ! Heat content per unit mass of rain (J/kg)
+      zcptrain(:,:) = rcp * ( SUM( (tn_ice(:,:,:) - rt0) * a_i(:,:,:), dim=3 ) + sst_m(:,:) * p_frld(:,:) ) 
+
       ! --- enthalpy of snow precip over ice in J/m3 (to be used in 1D-thermo) --- !
-      zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus )
+      zqprec_ice(:,:) = rhosn * ( zcptsnw(:,:) - lfus )
 
       ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
       DO jl = 1, jpl
-         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but we do not have Tice, so we consider Tice=0degC
+         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but atm. does not take it into account
       END DO
 
+      ! --- heat flux associated with emp (W/m2) --- !
+      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn   (:,:)   &        ! evap
+         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptrain(:,:)   &        ! liquid precip
+         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * zqprec_ice(:,:) * r1_rhosn ! solid precip over ocean + snow melting
+      zqemp_ice(:,:) =     zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhosn ! solid precip over ice (qevap_ice=0 since atm. does not take it into account)
+!!    zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptsnw (:,:)   &        ! ice evap
+!!       &             +   zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhosn ! solid precip over ice
+      
       ! --- total non solar flux (including evap/precip) --- !
       zqns_tot(:,:) = zqns_tot(:,:) + zqemp_ice(:,:) + zqemp_oce(:,:)
@@ -1974,5 +1980,5 @@
       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
 
-      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
+      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zcptrain, zcptsnw, zicefr, zmsk, zsnw )
       CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
       CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
Index: /trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7812)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7813)
@@ -323,5 +323,5 @@
                   zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) )
                   !                                           ! TKE Langmuir circulation source term
-                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * (1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   &
+                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   &
                      &                              / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
@@ -454,5 +454,5 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   &
-                     &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                     &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -463,5 +463,5 @@
                jk = nmln(ji,jj)
                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   &
-                  &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                  &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
             END DO
          END DO
@@ -476,5 +476,5 @@
                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   &
-                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                     &                        * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
