Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 7813)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 7814)
@@ -419,5 +419,5 @@
    FUNCTION ice_alloc()
       !!-----------------------------------------------------------------
-      !!               *** Routine ice_alloc_2 ***
+      !!               *** Routine ice_alloc ***
       !!-----------------------------------------------------------------
       INTEGER :: ice_alloc
@@ -506,5 +506,5 @@
 
       ice_alloc = MAXVAL( ierr(:) )
-      IF( ice_alloc /= 0 )   CALL ctl_warn('ice_alloc_2: failed to allocate arrays.')
+      IF( ice_alloc /= 0 )   CALL ctl_warn('ice_alloc: failed to allocate arrays.')
       !
    END FUNCTION ice_alloc
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 7813)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 7814)
@@ -259,4 +259,5 @@
                   closing_net(ji,jj) = 0._wp
                   opning     (ji,jj) = 0._wp
+                  ato_i      (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) )
                ELSE
                   iterate_ridging    = 1
@@ -844,5 +845,5 @@
          END DO
    
-         strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:)
+         strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1)
                          ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
          ksmooth = 1
@@ -853,5 +854,5 @@
       ELSE                      ! kstrngth ne 1:  Hibler (1979) form
          !
-         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  )
+         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1)
          !
          ksmooth = 1
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 7813)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 7814)
@@ -348,4 +348,5 @@
          &                   '   control of time parameter  nrstdt' )
 
+      ! Prognostic variables 
       DO jl = 1, jpl 
          WRITE(zchar,'(I2)') jl
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 7813)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 7814)
@@ -62,6 +62,7 @@
 
       IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*) ' lim_update1 ' 
-         WRITE(numout,*) ' ~~~~~~~~~~~ '
+         WRITE(numout,*)'' 
+         WRITE(numout,*)' lim_update1 ' 
+         WRITE(numout,*)' ~~~~~~~~~~~ '
       ENDIF
 
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 7813)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 7814)
@@ -62,6 +62,7 @@
 
       IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*) ' lim_update2 '
-         WRITE(numout,*) ' ~~~~~~~~~~~ '
+         WRITE(numout,*)''
+         WRITE(numout,*)' lim_update2 '
+         WRITE(numout,*)' ~~~~~~~~~~~ '
       ENDIF
 
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 7813)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 7814)
@@ -27,5 +27,5 @@
    !!                        - et_i(jpi,jpj)  !total ice thermal content 
    !!                        - smt_i(jpi,jpj) !mean ice salinity
-   !!                        - ot_i(jpi,jpj)  !average ice age
+   !!                        - tm_i (jpi,jpj) !mean ice temperature
    !!======================================================================
    !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
@@ -400,4 +400,5 @@
       !
    END SUBROUTINE lim_var_salprof
+
 
    SUBROUTINE lim_var_bv
@@ -654,5 +655,5 @@
       INTEGER  :: ji, jk, jl             ! dummy loop indices
       INTEGER  :: ijpij, i_fill, jl0  
-      REAL(wp) :: zarg, zV, zconv, zdh
+      REAL(wp) :: zarg, zV, zconv, zdh, zdv
       REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
@@ -675,78 +676,92 @@
          IF( zhti(ji) > 0._wp ) THEN
 
-         ! initialisation of tests
-         itest(:)  = 0
+            ! find which category (jl0) the input ice thickness falls into
+            jl0 = jpl
+            DO jl = 1, jpl
+               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
+                  jl0 = jl
+                  CYCLE
+               ENDIF
+            END DO
+
+            ! initialisation of tests
+            itest(:)  = 0
          
-         i_fill = jpl + 1                                             !====================================
-         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories  
-            ! iteration                                               !====================================
-            i_fill = i_fill - 1
+            i_fill = jpl + 1                                             !====================================
+            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
+               ! iteration                                               !====================================
+               i_fill = i_fill - 1
+               
+               ! initialisation of ice variables for each try
+               zht_i(ji,1:jpl) = 0._wp
+               za_i (ji,1:jpl) = 0._wp
+               itest(:)        = 0      
+               
+               ! *** case very thin ice: fill only category 1
+               IF ( i_fill == 1 ) THEN
+                  zht_i(ji,1) = zhti(ji)
+                  za_i (ji,1) = zai (ji)
+                  
+               ! *** case ice is thicker: fill categories >1
+               ELSE
+
+                  ! Fill ice thicknesses in the (i_fill-1) cat by hmean 
+                  DO jl = 1, i_fill - 1
+                     zht_i(ji,jl) = hi_mean(jl)
+                  END DO
+                  
+                  ! Concentrations in the (i_fill-1) categories 
+                  za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
+                  DO jl = 1, i_fill - 1
+                     IF ( jl /= jl0 ) THEN
+                        zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
+                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
+                     ENDIF
+                  END DO
+                  
+                  ! Concentration in the last (i_fill) category
+                  za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
+                  
+                  ! Ice thickness in the last (i_fill) category
+                  zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
+                  zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
+                  
+                  ! clem: 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 ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
+                           zdv = zht_i(ji,jl) * za_i(ji,jl)
+                           zht_i(ji,jl    ) = 0._wp
+                           za_i (ji,jl    ) = 0._wp
+                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
+                        END IF
+                     ENDDO
+                  ENDIF
+               
+               ENDIF ! case ice is thick or thin
             
-            ! initialisation of ice variables for each try
-            zht_i(ji,1:jpl) = 0._wp
-            za_i (ji,1:jpl) = 0._wp
-            itest(:)        = 0      
-      
-            ! *** case very thin ice: fill only category 1
-            IF ( i_fill == 1 ) THEN
-               zht_i(ji,1) = zhti(ji)
-               za_i (ji,1) = zai (ji)
-
-            ! *** case ice is thicker: fill categories >1
-            ELSE
-
-               ! Fill ice thicknesses except the last one (i_fill) by hmean 
-               DO jl = 1, i_fill - 1
-                  zht_i(ji,jl) = hi_mean(jl)
-               END DO
+               !---------------------
+               ! Compatibility tests
+               !--------------------- 
+               ! Test 1: area conservation
+               zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
+               IF ( zconv < epsi06 ) itest(1) = 1
+            
+               ! Test 2: volume conservation
+               zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
+               IF ( zconv < epsi06 ) itest(2) = 1
                
-               ! find which category (jl0) the input ice thickness falls into
-               jl0 = i_fill
+               ! Test 3: thickness of the last category is in-bounds ?
+               IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
+               
+               ! Test 4: positivity of ice concentrations
+               itest(4) = 1
                DO jl = 1, i_fill
-                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
-                     jl0 = jl
-		     CYCLE
-                  ENDIF
-               END DO
-               
-               ! Concentrations in the (i_fill-1) categories 
-               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
-               DO jl = 1, i_fill - 1
-                  IF ( jl == jl0 ) CYCLE
-                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
-                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
-               END DO
-               
-               ! Concentration in the last (i_fill) category
-               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
-               
-               ! Ice thickness in the last (i_fill) category
-               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
-               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 
-               
-            ENDIF ! case ice is thick or thin
-            
-            !---------------------
-            ! Compatibility tests
-            !--------------------- 
-            ! Test 1: area conservation
-            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
-            IF ( zconv < epsi06 ) itest(1) = 1
-            
-            ! Test 2: volume conservation
-            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
-            IF ( zconv < epsi06 ) itest(2) = 1
-            
-            ! Test 3: thickness of the last category is in-bounds ?
-            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
-            
-            ! Test 4: positivity of ice concentrations
-            itest(4) = 1
-            DO jl = 1, i_fill
-               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
-            END DO            
-                                                           !============================
-         END DO                                            ! end iteration on categories
-                                                           !============================
+                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
+               END DO
+               !                                         !============================
+            END DO                                       ! end iteration on categories
+               !                                         !============================
          ENDIF ! if zhti > 0
       END DO ! i loop
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7813)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7814)
@@ -40,9 +40,8 @@
    REAL(wp), PUBLIC ::   rtt      = 273.16_wp        !: triple point of temperature   [Kelvin]
    REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin]
+   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin]
 #if defined key_lim3
-   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin]
    REAL(wp), PUBLIC ::   rt0_ice  = 273.15_wp        !: melting point of ice          [Kelvin]
 #else
-   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin]
    REAL(wp), PUBLIC ::   rt0_ice  = 273.05_wp        !: melting point of ice          [Kelvin]
 #endif
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 7813)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 7814)
@@ -113,9 +113,4 @@
       CASE( 0 )
        
-         !ralb_sf = 0.80       ! dry snow
-         !ralb_sm = 0.65       ! melting snow
-         !ralb_if = 0.72       ! bare frozen ice
-         !ralb_im = ...        ! bare puddled ice 
-
          !  Computation of ice albedo (free of snow)
          WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
@@ -167,5 +162,4 @@
       CASE( 1 ) 
 
-!        ralb_im = ...        ! bare puddled ice
 ! compilation of values from literature
 !        ralb_sf = 0.85      ! dry snow
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7813)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7814)
@@ -1622,6 +1622,5 @@
       !
       IF( srcv(jpr_icb)%laction )  zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus ! remove heat content associated to iceberg melting
-      !
-!!      !
+
 
 #if defined key_lim3      
@@ -1632,29 +1631,26 @@
 
       ! 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 )
+      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(:,:) ) 
-
-      ! --- heat flux associated with emp (W/m2) --- !
-      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn(:,:)   &       ! evap
-         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptrain(:,:)   &       ! liquid precip
-         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptsnw(:,:) - 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             * ( zcptsnw(:,:) - lfus ) ! solid precip over ice (only)
-                                                                                                       ! qevap_ice=0 since we consider Tice=0degC
-      
+      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 * ( zcptsnw(:,:) - lfus )
-      !zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - 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(:,:)
