Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4622)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 5038)
@@ -5,5 +5,5 @@
    !!======================================================================
    !! History :  LIM  ! 2006-02  (M. Vancoppenolle) Original code 
-   !!            3.2  ! 2009-07  (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_mec
+   !!            3.2  ! 2009-07  (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn
    !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
    !!----------------------------------------------------------------------
@@ -22,5 +22,4 @@
    USE limthd_lac       ! LIM
    USE limvar           ! LIM
-   USE limcons          ! LIM
    USE in_out_manager   ! I/O manager
    USE lbclnk           ! lateral boundary condition - MPP exchanges
@@ -30,7 +29,8 @@
   ! Check budget (Rousset)
    USE iom              ! I/O manager
-   USE lib_fortran     ! glob_sum
+   USE lib_fortran      ! glob_sum
    USE limdiahsb
-   USE timing          ! Timing
+   USE timing           ! Timing
+   USE limcons          ! conservation tests
 
    IMPLICIT NONE
@@ -42,8 +42,4 @@
    PUBLIC   lim_itd_me_zapsmall
    PUBLIC   lim_itd_me_alloc        ! called by iceini.F90
-
-   REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values
-   REAL(wp) ::   epsi10 = 1.e-10_wp   ! constant values
-   REAL(wp) ::   epsi06 = 1.e-06_wp   ! constant values
 
    !-----------------------------------------------------------------------
@@ -143,16 +139,10 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2)
       REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories
-      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
-      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
-      ! mass and salt flux (clem)
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume...
+      !
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!-----------------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limitd_me')
 
       CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
-
-      CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem
-
-      IF( numit == nstart  )   CALL lim_itd_me_init   ! Initialization (first time-step only)
 
       IF(ln_ctl) THEN
@@ -162,19 +152,7 @@
 
       IF( ln_limdyn ) THEN          !   Start ridging and rafting   !
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) )
-         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) )
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
-
-      ! mass and salt flux init (clem)
-      zviold(:,:,:) = v_i(:,:,:)
-      zvsold(:,:,:) = v_s(:,:,:)
-      zsmvold(:,:,:) = smv_i(:,:,:)
+
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       !-----------------------------------------------------------------------------!
@@ -362,6 +340,6 @@
             ! 5) Heat, salt and freshwater fluxes
             !-----------------------------------------------------------------------------!
-            fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean
-            fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean
+            wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean
+            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * unit_fac / area(ji,jj) * r1_rdtice  ! heat sink for ocean (<0, W.m-2)
 
          END DO
@@ -399,17 +377,4 @@
       CALL lim_itd_me_zapsmall
 
-      !--------------------------------
-      ! Update mass/salt fluxes (clem)
-      !--------------------------------
-      DO jl = 1, jpl
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-               diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice
-               rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 
-               rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 
-               sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice 
-            END DO
-         END DO
-      END DO
 
       IF(ln_ctl) THEN     ! Control print
@@ -445,33 +410,10 @@
       ENDIF
 
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
-         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
- 
-         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice
-         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
-
-         zchk_vmin = glob_min(v_i)
-         zchk_amax = glob_max(SUM(a_i,dim=3))
-         zchk_amin = glob_min(a_i)
-       
-         IF(lwp) THEN
-            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * rday)
-            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday)
-            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3)
-            IF ( zchk_amax >  kamax+epsi10  ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax
-            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',zchk_amin
-         ENDIF
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       ENDIF  ! ln_limdyn=.true.
       !
       CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
-      !
-      CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem
       !
       IF( nn_timing == 1 )  CALL timing_stop('limitd_me')
@@ -670,5 +612,4 @@
       !!---------------------------------------------------------------------!
       INTEGER ::   ji,jj, jl    ! dummy loop indices
-      INTEGER ::   krdg_index   ! 
       REAL(wp) ::   Gstari, astari, hi, hrmean, zdummy   ! local scalar
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka    ! temporary array used here
@@ -746,8 +687,6 @@
       !-----------------------------------------------------------------
 
-      krdg_index = 1
-
-      IF( krdg_index == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
-         DO jl = 0, ice_cat_bounds(1,2)       ! only undeformed ice participates
+      IF( partfun_swi == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
+         DO jl = 0, jpl    
             DO jj = 1, jpj 
                DO ji = 1, jpi
@@ -772,11 +711,11 @@
             Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
          END DO !jl
-         DO jl = 0, ice_cat_bounds(1,2)
+         DO jl = 0, jpl
              athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
          END DO
          !
-      ENDIF ! krdg_index
-
-      IF( raftswi == 1 ) THEN      ! Ridging and rafting ice participation functions
+      ENDIF ! partfun_swi
+
+      IF( raft_swi == 1 ) THEN      ! Ridging and rafting ice participation functions
          !
          DO jl = 1, jpl
@@ -794,5 +733,5 @@
          END DO ! jl
 
-      ELSE  ! raftswi = 0
+      ELSE  ! raft_swi = 0
          !
          DO jl = 1, jpl
@@ -802,5 +741,5 @@
       ENDIF
 
-      IF ( raftswi == 1 ) THEN
+      IF ( raft_swi == 1 ) THEN
 
          IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN
@@ -908,6 +847,6 @@
       INTEGER ::   ij                ! horizontal index, combines i and j loops
       INTEGER ::   icells            ! number of cells with aicen > puny
-      REAL(wp) ::   zindb, zsrdg2   ! local scalar
       REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration
+      REAL(wp) ::   zsstK            ! SST in Kelvin
 
       INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
@@ -917,5 +856,5 @@
 
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   aicen_init, vicen_init   ! ice  area    & volume before ridging
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnon_init, esnon_init   ! snow volume  & energy before ridging
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnwn_init, esnwn_init   ! snow volume  & energy before ridging
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   smv_i_init, oa_i_init    ! ice salinity & age    before ridging
 
@@ -952,7 +891,7 @@
       CALL wrk_alloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw )
       CALL wrk_alloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
-      CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init )
-      CALL wrk_alloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw )
-      CALL wrk_alloc( jpi, jpj, jkmax, jpl, eicen_init )
+      CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
+      CALL wrk_alloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw )
+      CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
 
       ! Conservation check
@@ -1008,5 +947,5 @@
          aicen_init(:,:,jl) = a_i(:,:,jl)
          vicen_init(:,:,jl) = v_i(:,:,jl)
-         vsnon_init(:,:,jl) = v_s(:,:,jl)
+         vsnwn_init(:,:,jl) = v_s(:,:,jl)
          !
          smv_i_init(:,:,jl) = smv_i(:,:,jl)
@@ -1014,5 +953,5 @@
       END DO !jl
 
-      esnon_init(:,:,:) = e_s(:,:,1,:)
+      esnwn_init(:,:,:) = e_s(:,:,1,:)
 
       DO jl = 1, jpl  
@@ -1091,17 +1030,17 @@
             !     / rafting category n1.
             !--------------------------------------------------------------------------
-            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )
+            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)
             vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por )
             vsw  (ji,jj) = vrdg1(ji,jj) * ridge_por
 
-            vsrdg(ji,jj) = vsnon_init(ji,jj,jl1) * afrac(ji,jj)
-            esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj)
-            srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )
-            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
+            vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj)
+            esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj)
+            srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
+            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless
 
             ! rafting volumes, heat contents ...
             virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
-            vsrft(ji,jj) = vsnon_init(ji,jj,jl1) * afrft(ji,jj)
-            esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj)
+            vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj)
+            esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj)
             smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 
 
@@ -1120,14 +1059,11 @@
             ! Salinity
             !-------------
-            smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of seawater frozen in voids
-
-            zsrdg2       = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge
-
-            srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity
+            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014
+            srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge
+
+            !srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity
             
-            !                                                             ! excess of salt is flushed into the ocean
-            !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice
-
-            !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic    ! gurvan: increase in ice volume du to seawater frozen in voids             
+            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice
+            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice   ! gurvan: increase in ice volume du to seawater frozen in voids             
 
             !------------------------------------            
@@ -1158,6 +1094,7 @@
                &                                + rhosn*vsrft(ji,jj)*(1.0-fsnowrft)
 
-            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included
-               &                                + esrft(ji,jj)*(1.0-fsnowrft)          
+            ! in 1e-9 Joules (same as e_s)
+            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included
+               &                                - esrft(ji,jj)*(1.0-fsnowrft)          
 
             !-----------------------------------------------------------------
@@ -1184,29 +1121,27 @@
                jj = indxj(ij)
                ! heat content of ridged ice
-               erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 
+               erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) 
                eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
                e_i  (ji,jj,jk,jl1)  = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)
-               ! sea water heat content
-               ztmelts          = - tmut * sss_m(ji,jj) + rtt
-               ! heat content per unit volume
-               zdummy0          = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj)
-
-               ! corrected sea water salinity
-               zindb  = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - epsi20 ) )
-               zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), epsi20 )
-
-               ztmelts          = - tmut * zdummy + rtt
-               ersw(ji,jj,jk)   = - rcp * ( ztmelts - rtt ) * vsw(ji,jj)
-
-               ! heat flux
-               fheat_mec(ji,jj) = fheat_mec(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) * r1_rdtice
+               
+               
+               ! enthalpy of the trapped seawater (J/m2, >0)
+               ! clem: if sst>0, then ersw <0 (is that possible?)
+               zsstK  = sst_m(ji,jj) + rt0
+               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i )
+
+               ! heat flux to the ocean
+               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux 
 
                ! Correct dimensions to avoid big values
-               ersw(ji,jj,jk)   = ersw(ji,jj,jk) * 1.e-09
-
-               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
-               ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i )
+               ersw(ji,jj,jk)   = ersw(ji,jj,jk) / unit_fac
+
+               ! Mutliply by ice volume, and divide by number of layers to get heat content in 1.e9 J
+               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 
+               !! MV HC 2014
+               ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj)
 
                erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
+
             END DO ! ij
          END DO !jk
@@ -1253,5 +1188,5 @@
          !-------------------------------------------------------------------------------
          !        jl1 looping 1-jpl
-         DO jl2  = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
+         DO jl2  = 1, jpl 
             ! over categories to which ridged ice is transferred
 !CDIR NODEP
@@ -1298,5 +1233,5 @@
          END DO                 ! jl2 (new ridges)            
 
-         DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
+         DO jl2 = 1, jpl 
 
 !CDIR NODEP
@@ -1361,7 +1296,7 @@
       CALL wrk_dealloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw )
       CALL wrk_dealloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
-      CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init )
-      CALL wrk_dealloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw )
-      CALL wrk_dealloc( jpi, jpj, jkmax, jpl, eicen_init )
+      CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
+      CALL wrk_dealloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw )
+      CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
       !
    END SUBROUTINE lim_itd_me_ridgeshift
@@ -1404,9 +1339,7 @@
       !!-------------------------------------------------------------------
       INTEGER :: ios                 ! Local integer output status for namelist read
-      NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,& 
-         Gstar, astar,                                &
-         Hstar, raftswi, hparmeter, Craft, ridge_por, &
-         sal_max_ridge,  partfun_swi, transfun_swi,   &
-         brinstren_swi
+      NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,              & 
+        &                   Gstar, astar, Hstar, raft_swi, hparmeter, Craft, ridge_por, &
+        &                   partfun_swi, brinstren_swi
       !!-------------------------------------------------------------------
       !
@@ -1418,5 +1351,5 @@
       READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
-      WRITE ( numoni, namiceitdme )
+      IF(lwm) WRITE ( numoni, namiceitdme )
       !
       IF (lwp) THEN                          ! control print
@@ -1432,11 +1365,9 @@
          WRITE(numout,*)'   Equivalent to G* for an exponential part function       astar           ', astar
          WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     Hstar           ', Hstar
-         WRITE(numout,*)'   Rafting of ice sheets or not                            raftswi         ', raftswi
+         WRITE(numout,*)'   Rafting of ice sheets or not                            raft_swi        ', raft_swi
          WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       hparmeter       ', hparmeter
          WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  Craft           ', Craft  
          WRITE(numout,*)'   Initial porosity of ridges                              ridge_por       ', ridge_por
-         WRITE(numout,*)'   Maximum salinity of ridging ice                         sal_max_ridge   ', sal_max_ridge
          WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    partfun_swi     ', partfun_swi
-         WRITE(numout,*)'   Switch for tran. function (0) linear (1) exponential    transfun_swi    ', transfun_swi
          WRITE(numout,*)'   Switch for including brine volume in ice strength comp. brinstren_swi   ', brinstren_swi
       ENDIF
@@ -1462,5 +1393,5 @@
 
       REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace
-      REAL(wp)                          ::   zmask_glo
+      REAL(wp)                          ::   zmask_glo, zsal, zvi, zvs, zei, zes
 !!gm      REAL(wp) ::   xtmp      ! temporary variable
       !!-------------------------------------------------------------------
@@ -1468,9 +1399,10 @@
       CALL wrk_alloc( jpi, jpj, zmask )
 
+      ! to be sure that at_i is the sum of a_i(jl)
+      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
+
       DO jl = 1, jpl
-
          !-----------------------------------------------------------------
          ! Count categories to be zapped.
-         ! Abort model in case of negative area.
          !-----------------------------------------------------------------
          icells = 0
@@ -1478,8 +1410,7 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   &
-                  & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   &
-                  & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   &
-                  & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp
+               IF( a_i(ji,jj,jl) <= epsi10 .OR. v_i(ji,jj,jl) <= epsi10 .OR. at_i(ji,jj) <= epsi10 ) THEN
+                  zmask(ji,jj) = 1._wp
+               ENDIF
             END DO
          END DO
@@ -1494,8 +1425,9 @@
             DO jj = 1 , jpj
                DO ji = 1 , jpi
-!!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) * r1_rdtice
-!!gm                  xtmp = xtmp * unit_fac
-                  ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
+                  zei  = e_i(ji,jj,jk,jl)
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) )
+                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) + rtt * zmask(ji,jj)
+                  ! update exchanges with ocean
+                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
                END DO
             END DO
@@ -1504,5 +1436,9 @@
          DO jj = 1 , jpj
             DO ji = 1 , jpi
-
+               
+               zsal = smv_i(ji,jj,jl)
+               zvi  = v_i(ji,jj,jl)
+               zvs  = v_s(ji,jj,jl)
+               zes  = e_s(ji,jj,1,jl)
                !-----------------------------------------------------------------
                ! Zap snow energy and use ocean heat to melt snow
@@ -1514,9 +1450,4 @@
                ! fluxes are positive to the ocean
                ! here the flux has to be negative for the ocean
-!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice
-               !           fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
-
-!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ???????
-
                t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) )
 
@@ -1524,13 +1455,5 @@
                ! zap ice and snow volume, add water and salt to ocean
                !-----------------------------------------------------------------
-
-               !           xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt
-               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj)                  )   &
-               !                                            * rhosn * v_s(ji,jj,jl) * r1_rdtice
-               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) )   & 
-               !                                            * rhoic * v_i(ji,jj,jl) * r1_rdtice
-               !           sfx (i,j)      = sfx (i,j)      + xtmp
-
-               ato_i(ji,jj)    = a_i  (ji,jj,jl) *       zmask(ji,jj)   + ato_i(ji,jj)
+               ato_i(ji,jj)    = a_i  (ji,jj,jl) *           zmask(ji,jj)   + ato_i(ji,jj)
                a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
                v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
@@ -1539,9 +1462,21 @@
                oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
                smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
-               !
+               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) )
+               ! additional condition
+               IF( v_s(ji,jj,jl) <= epsi10 ) THEN
+                  v_s(ji,jj,jl)   = 0._wp
+                  e_s(ji,jj,1,jl) = 0._wp
+               ENDIF
+               ! update exchanges with ocean
+               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
+               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
+               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
+               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
             END DO
          END DO
-         !
-      END DO                 ! jl 
+      END DO ! jl 
+
+      ! to be sure that at_i is the sum of a_i(jl)
+      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
       !
       CALL wrk_dealloc( jpi, jpj, zmask )
