Index: branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 7276)
+++ branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 7278)
@@ -74,5 +74,5 @@
 
       REAL(wp) ::   ztmelts             ! local scalar
-      REAL(wp) ::   zfdum       
+      REAL(wp) ::   zdum       
       REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
       REAL(wp) ::   zs_snic      ! snow-ice salinity
@@ -95,4 +95,5 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
       REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
+      REAL(wp), POINTER, DIMENSION(:) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2)
 
       REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt 
@@ -105,6 +106,4 @@
 
       REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2)
-      REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2)
-      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3)
       REAL(wp), POINTER, DIMENSION(:) ::   zsnw        ! distribution of snow after wind blowing
 
@@ -118,20 +117,19 @@
       ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values)
       SELECT CASE( nn_icesal )                       ! varying salinity or not
-         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile
-         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile
+         CASE( 1, 3 ) ;   zswitch_sal = 0       ! prescribed salinity profile
+         CASE( 2 )    ;   zswitch_sal = 1       ! varying salinity profile
       END SELECT
 
-      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw )
-      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
+      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
+      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i )
       CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i )
       CALL wrk_alloc( jpij, nlay_i, icount )
        
-      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp
+      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp ; dh_i_sub(:) = 0._wp
       dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp   
 
       zqprec   (:) = 0._wp ; zq_su    (:) = 0._wp ; zq_bo    (:) = 0._wp ; zf_tt(:) = 0._wp
-      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp
+      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp ; zevap_rema(:) = 0._wp ;
       zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp
-      zqh_s    (:) = 0._wp ; zq_s     (:) = 0._wp     
 
       zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp       
@@ -159,8 +157,8 @@
       !
       DO ji = kideb, kiut
-         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
+         zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
          zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
 
-         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
+         zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
          zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
       END DO
@@ -187,10 +185,4 @@
       !  2) Computing layer thicknesses and enthalpies.            !
       !------------------------------------------------------------!
-      !
-      DO jk = 1, nlay_s
-         DO ji = kideb, kiut
-            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s
-         END DO
-      END DO
       !
       DO jk = 1, nlay_i
@@ -275,17 +267,15 @@
       END DO
 
-      !----------------------
-      ! 3.2 Snow sublimation 
-      !----------------------
+      !------------------------------
+      ! 3.2 Sublimation (part1: snow) 
+      !------------------------------
       ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
       ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean)
-      ! clem comment: ice should also sublimate
       zdeltah(:,:) = 0._wp
-      ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice
-      ! forced  mode: snow thickness change due to sublimation
-      DO ji = kideb, kiut
-         zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
-         ! Heat flux by sublimation [W.m-2], < 0
-         !      sublimate first snow that had fallen, then pre-existing snow
+      DO ji = kideb, kiut
+         zdh_s_sub(ji)  = MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
+         ! remaining evap in kg.m-2 (used for ice melting later on)
+         zevap_rema(ji)  = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn
+         ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow)
          zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
          hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  &
@@ -309,12 +299,10 @@
       !-------------------------------------------
       ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
-      zq_s(:) = 0._wp 
       DO jk = 1, nlay_s
          DO ji = kideb,kiut
-            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  )
-            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          &
-              &            ( (   zdh_s_pre(ji)             ) * zqprec(ji) +  &
-              &              (   ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
-            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk)
+            rswitch       = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) )
+            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *           &
+              &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  &
+              &              ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
          END DO
       END DO
@@ -370,5 +358,5 @@
                zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
                
-               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
+               ! Contribution to salt flux >0 (clem: using sm_i_1d and not s_i_1d(jk) is ok)
                sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
                
@@ -383,4 +371,20 @@
                
             END IF
+            ! ----------------------
+            ! Sublimation part2: ice
+            ! ----------------------
+            zdum      = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic )
+            zdeltah(ji,jk) = zdeltah(ji,jk) + zdum
+            dh_i_sub(ji)  = dh_i_sub(ji) + zdum
+            ! Salt flux > 0 (clem2016: flux is sent to the ocean for simplicity but salt should remain in the ice except if all ice is melted.
+            !                          It must be corrected at some point)
+            sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * sm_i_1d(ji) * r1_rdtice
+            ! Heat flux [W.m-2], < 0
+            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * q_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
+            ! Mass flux > 0
+            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
+            ! update remaining mass flux
+            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic
+            
             ! record which layers have disappeared (for bottom melting) 
             !    => icount=0 : no layer has vanished
@@ -389,5 +393,5 @@
             icount(ji,jk) = NINT( rswitch )
             zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
-
+                        
             ! update heat content (J.m-2) and layer thickness
             qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
@@ -397,5 +401,11 @@
       ! update ice thickness
       DO ji = kideb, kiut
-         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) )
+         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) )
+      END DO
+
+      ! remaining "potential" evap is sent to ocean
+      DO ji = kideb, kiut
+         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
+         wfx_err_sub(ii,ij) = wfx_err_sub(ii,ij) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice  ! <=0 (net evap for the ocean in kg.m-2.s-1)
       END DO
 
@@ -653,4 +663,10 @@
          sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 
           
+         ! virtual salt flux to keep salinity constant
+         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
+            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_m(ii,ij) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m into the ocean
+               &                            - sm_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  sm_i  from the ocean 
+         ENDIF
+
          ! Contribution to mass flux
          ! All snow is thrown in the ocean, and seawater is taken to replace the volume
@@ -686,6 +702,6 @@
       WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp
       
-      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw )
-      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
+      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
+      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i )
       CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i )
       CALL wrk_dealloc( jpij, nlay_i, icount )
