Index: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5146)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5167)
@@ -94,5 +94,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns
       !!-------------------------------------------------------------------
-      CALL wrk_alloc( jpi, jpj, zqsr, zqns )
+      CALL wrk_alloc( jpi,jpj, zqsr, zqns )
 
       IF( nn_timing == 1 )  CALL timing_start('limthd')
@@ -101,4 +101,5 @@
       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
+      CALL lim_var_glo2eqv
       !------------------------------------------------------------------------!
       ! 1) Initialization of some variables                                    !
@@ -209,6 +210,6 @@
             ! Net heat flux on top of ice-ocean [W.m-2]
             ! -----------------------------------------
-            !     First  step here      : heat flux at the ocean surface + precip
-            !     Second step below     : heat flux at the ice   surface (after limthd_dif) 
+            !     heat flux at the ocean surface + precip
+            !   + heat flux at the ice   surface 
             hfx_in(ji,jj) = hfx_in(ji,jj)                                                                                         & 
                ! heat flux above the ocean
@@ -216,5 +217,7 @@
                ! latent heat of precip (note that precip is included in qns but not in qns_ice)
                &    +   ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  &
-               &    +   ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )
+               &    +   ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )          &
+               ! heat flux above the ice
+               &    +   SUM(    a_i_b(ji,jj,:)   * ( qns_ice(ji,jj,:) + qsr_ice(ji,jj,:) ) )
 
             ! -----------------------------------------------------------------------------
@@ -226,5 +229,5 @@
             hfx_out(ji,jj) = hfx_out(ji,jj)                                                                                       & 
                ! Non solar heat flux received by the ocean
-               &    +        pfrld(ji,jj) * qns(ji,jj)                                                                            &
+               &    +        pfrld(ji,jj) * zqns(ji,jj)                                                                            &
                ! latent heat of precip (note that precip is included in qns but not in qns_ice)
                &    +      ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)       &
@@ -311,5 +314,5 @@
             ! --- lateral melting if monocat --- !
             !------------------------------------!
-            IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN
+            IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN
                CALL lim_thd_lam( 1, nbpb )
             END IF
@@ -324,5 +327,5 @@
          ENDIF
          !
-      END DO
+      END DO !jl
 
       !------------------------------------------------------------------------------!
@@ -358,5 +361,7 @@
       ! Change thickness to volume
       !----------------------------------
-      CALL lim_var_eqv2glo
+      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
+      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
+      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
 
       CALL lim_var_zapsmall
@@ -399,30 +404,31 @@
       !
       !
-      CALL wrk_dealloc( jpi, jpj, zqsr, zqns )
-
       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
+      CALL wrk_dealloc( jpi,jpj, zqsr, zqns )
+
       !------------------------------------------------------------------------------|
       !  6) Transport of ice between thickness categories.                           |
       !------------------------------------------------------------------------------|
+      ! Given thermodynamic growth rates, transport ice between thickness categories.
       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-      ! Given thermodynamic growth rates, transport ice between thickness categories.
-      IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt )
-      !
-      CALL lim_var_glo2eqv    ! only for info
-      CALL lim_var_agg(1)
+      IF( jpl > 1 )      CALL lim_itd_th_rem( 1, jpl, kt )
 
       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
       !------------------------------------------------------------------------------|
       !  7) Add frazil ice growing in leads.
       !------------------------------------------------------------------------------|
       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
       CALL lim_thd_lac
-      CALL lim_var_glo2eqv    ! only for info
       
-      ! conservation test
       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-      IF(ln_ctl) THEN   ! Control print
+      ! Control print
+      IF(ln_ctl) THEN
+         CALL lim_var_glo2eqv
+
          CALL prt_ctl_info(' ')
          CALL prt_ctl_info(' - Cell values : ')
@@ -503,17 +509,19 @@
       REAL(wp)            ::   zhi_bef            ! ice thickness before thermo
       REAL(wp)            ::   zdh_mel, zda_mel   ! net melting
-      REAL(wp)            ::   zv                 ! ice volume 
+      REAL(wp)            ::   zvi, zvs           ! ice/snow volumes 
 
       DO ji = kideb, kiut
          zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) )
-         IF( zdh_mel < 0._wp )  THEN
-            zv         = a_i_1d(ji) * ht_i_1d(ji)
+         IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN
+            zvi          = a_i_1d(ji) * ht_i_1d(ji)
+            zvs          = a_i_1d(ji) * ht_s_1d(ji)
             ! lateral melting = concentration change
             zhi_bef     = ht_i_1d(ji) - zdh_mel
-            zda_mel     =  a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) )
-            a_i_1d(ji)  = MAX( 0._wp, a_i_1d(ji) + zda_mel ) 
-            ! adjust thickness
-            rswitch     = MAX( 0._wp , SIGN( 1._wp , a_i_1d(ji) - epsi20 ) )
-            ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 )
+            rswitch     = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) )
+            zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
+            a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel ) 
+             ! adjust thickness
+            ht_i_1d(ji) = zvi / a_i_1d(ji)            
+            ht_s_1d(ji) = zvs / a_i_1d(ji)            
             ! retrieve total concentration
             at_i_1d(ji) = a_i_1d(ji)
@@ -676,5 +684,5 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       &
-         &                rn_himin, parsub, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &
+         &                rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &
          &                nn_monocat, ln_it_qnsice
       !!-------------------------------------------------------------------
@@ -700,5 +708,4 @@
       ENDIF
 
-      IF( lk_cpl .AND. parsub /= 0.0 )   CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' )
       !
       IF(lwp) THEN                          ! control print
@@ -712,5 +719,4 @@
          WRITE(numout,*)'      minimum ice thickness                                   rn_himin     = ', rn_himin 
          WRITE(numout,*)'      numerical carac. of the scheme for diffusion in ice '
-         WRITE(numout,*)'      switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub  
          WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          rn_betas     = ', rn_betas
          WRITE(numout,*)'      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
