Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icealb.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icealb.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icealb.F90	(revision 9937)
@@ -148,5 +148,5 @@
                !
                !                       !--- Snow-covered ice albedo (freezing, melting cases)
-               IF( pt_su(ji,jj,jl) < rt0_snow ) THEN
+               IF( pt_su(ji,jj,jl) < rt0 ) THEN
                   zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw(ji,jj,jl) * z1_c3 )
                ELSE
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icecor.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icecor.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icecor.F90	(revision 9937)
@@ -86,5 +86,5 @@
       IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        !
       !                             !-----------------------------------------------------
-         zzc = rhoic * r1_Dt_ice
+         zzc = rhoi * r1_Dt_ice
          DO jl = 1, jpl
             DO jj = 1, jpj 
@@ -139,7 +139,7 @@
                   &                  + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_Dt_ice
                !                 ! salt, volume
-               diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_Dt_ice
-               diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_Dt_ice
-               diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_Dt_ice
+               diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoi * r1_Dt_ice
+               diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoi * r1_Dt_ice
+               diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhos * r1_Dt_ice
             END DO
          END DO
@@ -160,7 +160,7 @@
                   &                                   + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_Dt_ice
                !                 ! salt, volume
-               diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_Dt_ice
-               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_Dt_ice
-               diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_Dt_ice
+               diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoi * r1_Dt_ice
+               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoi * r1_Dt_ice
+               diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhos * r1_Dt_ice
             END DO
          END DO
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icectl.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icectl.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icectl.F90	(revision 9937)
@@ -93,7 +93,7 @@
             &                  ) *  e1e2t(:,:) ) * zconv
 
-         pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv )
-
-         pdiag_s = glob_sum( SUM( sv_i * rhoic             , dim=3 ) * e1e2t * zconv )
+         pdiag_v = glob_sum( SUM( v_i  * rhoi + v_s * rhos, dim=3 ) * e1e2t * zconv )
+
+         pdiag_s = glob_sum( SUM( sv_i * rhoi             , dim=3 ) * e1e2t * zconv )
 
          pdiag_t = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     &
@@ -120,8 +120,8 @@
  
          ! outputs
-         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv  &
+         zv = ( ( glob_sum( SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) * zconv  &
             &     - pdiag_v ) * r1_Dt_ice - zfv ) * rday
 
-         zs = ( ( glob_sum( SUM( sv_i * rhoic             , dim=3 ) * e1e2t ) * zconv  &
+         zs = ( ( glob_sum( SUM( sv_i * rhoi             , dim=3 ) * e1e2t ) * zconv  &
             &     - pdiag_s ) * r1_Dt_ice + zfs ) * rday
 
@@ -131,6 +131,6 @@
 
          ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative
-         zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t  ) * zconv * rday 
-         zetrp = glob_sum( ( diag_trp_ei         + diag_trp_es         ) * e1e2t  ) * zconv
+         zvtrp = glob_sum( ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t  ) * zconv * rday 
+         zetrp = glob_sum( ( diag_trp_ei        + diag_trp_es        ) * e1e2t  ) * zconv
 
          zvmin = glob_min( v_i )
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedia.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedia.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedia.F90	(revision 9937)
@@ -110,7 +110,7 @@
       ! 3 -  Content variations !
       ! ----------------------- !
-      zdiff_vol = r1_rho0 * glob_sum( ( rhoic*vt_i(:,:) + rhosn*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3) 
-      zdiff_sal = r1_rho0 * glob_sum( ( rhoic* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! salt content trend (km3*pss)
-      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:)             - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20 ! heat content trend (1.e20 J)
+      zdiff_vol = r1_rho0 * glob_sum( ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9  ! freshwater trend (km3) 
+      zdiff_sal = r1_rho0 * glob_sum( ( rhoi* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9  ! salt content trend (km3*pss)
+      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:)           - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20 ! heat content trend (1.e20 J)
       !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check)
 
@@ -241,7 +241,7 @@
             frc_sal     = 0._wp                                                 
             ! record initial ice volume, salt and temp
-            vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2)
-            tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J)
-            sal_loc_ini(:,:) = rhoic * SUM( sv_i(:,:,:), dim=3 )      ! ice salt content (pss*kg/m2)
+            vol_loc_ini(:,:) = rhoi * vt_i(:,:) + rhos * vt_s(:,:)   ! ice/snow volume (kg/m2)
+            tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                 ! ice/snow heat content (J)
+            sal_loc_ini(:,:) = rhoi * SUM( sv_i(:,:,:), dim=3 )      ! ice salt content (pss*kg/m2)
          ENDIF
          !
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_adv.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_adv.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_adv.F90	(revision 9937)
@@ -103,7 +103,7 @@
       diag_trp_vi(:,:) = SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice
       diag_trp_vs(:,:) = SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice
-      IF( iom_use('icemtrp') )   CALL iom_put( "icemtrp" , diag_trp_vi * rhoic          )   ! ice mass transport
-      IF( iom_use('snwmtrp') )   CALL iom_put( "snwmtrp" , diag_trp_vs * rhosn          )   ! snw mass transport
-      IF( iom_use('salmtrp') )   CALL iom_put( "salmtrp" , diag_trp_sv * rhoic * 1.e-03 )   ! salt mass transport (kg/m2/s)
+      IF( iom_use('icemtrp') )   CALL iom_put( "icemtrp" , diag_trp_vi * rhoi           )   ! ice mass transport
+      IF( iom_use('snwmtrp') )   CALL iom_put( "snwmtrp" , diag_trp_vs * rhos           )   ! snw mass transport
+      IF( iom_use('salmtrp') )   CALL iom_put( "salmtrp" , diag_trp_sv * rhoi * 1.e-03  )   ! salt mass transport (kg/m2/s)
       IF( iom_use('dihctrp') )   CALL iom_put( "dihctrp" , -diag_trp_ei                 )   ! advected ice heat content (W/m2)
       IF( iom_use('dshctrp') )   CALL iom_put( "dshctrp" , -diag_trp_es                 )   ! advected snw heat content (W/m2)
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_rdgrft.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_rdgrft.F90	(revision 9937)
@@ -543,5 +543,5 @@
                ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
                vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg
-               ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
+               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?)
 
                ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
@@ -570,12 +570,12 @@
 
                ! Ice-ocean exchanges associated with ice porosity
-               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw              * rhoic * r1_Dt_ice   ! increase in ice volume due to seawater frozen in voids
-               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_Dt_ice
-               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji)                 * r1_Dt_ice   ! > 0 [W.m-2] 
+               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw              * rhoi * r1_Dt_ice   ! increase in ice volume due to seawater frozen in voids
+               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_Dt_ice
+               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji)                * r1_Dt_ice   ! > 0 [W.m-2] 
 
                ! Put the snow lost by ridging into the ocean
                !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow.
-               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean
-                  &                                      + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice
+               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhos * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean
+                  &                                      + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice
 
                ! Put the melt pond water into the ocean
@@ -583,6 +583,6 @@
                !       is no net mass flux between melt ponds and the ocean (see icethd_pnd.F90 for ex.)
                !IF ( ln_pnd_fwb ) THEN
-               !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean
-               !      &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_Dt_ice
+               !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhow * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean
+               !      &                              + rhow * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_Dt_ice
                !ENDIF
 
@@ -590,6 +590,6 @@
                IF( nn_icesal /= 2 )  THEN
                   sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i
-                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoic * r1_Dt_ice  &  ! put back sss_m into the ocean
-                     &                            - s_i_1d(ji) * vsw * rhoic * r1_Dt_ice     ! and get  s_i  from the ocean 
+                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_Dt_ice  &  ! put back sss_m into the ocean
+                     &                            - s_i_1d(ji) * vsw * rhoi * r1_Dt_ice     ! and get  s_i  from the ocean 
                ENDIF
 
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_rhg_evp.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_rhg_evp.F90	(revision 9937)
@@ -285,7 +285,7 @@
 
             ! Ice/snow mass at U-V points
-            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
-            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
-            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
+            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
+            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
+            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
             zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
             zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
@@ -799,12 +799,12 @@
                zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
                
-               zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
-               zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
-               
-               zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
-               zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
-               
-               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
-               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
+               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) )   ! ice mass transport, X-component
+               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) )   !        ''           Y-   ''
+               
+               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) )   ! snow mass transport, X-component
+               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) )   !          ''          Y-   ''
+               
+               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )          ! area transport,      X-component
+               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )          !        ''            Y-   ''
                
             END DO
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/iceistate.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/iceistate.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/iceistate.F90	(revision 9937)
@@ -295,8 +295,8 @@
                   ! In case snow load is in excess that would lead to transformation from snow to ice
                   ! Then, transfer the snow excess into the ice (different from icethd_dh)
-                  zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 ) 
+                  zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 ) 
                   ! recompute h_i, h_s avoiding out of bounds values
                   h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh )
-                  h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
+                  h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi * r1_rhos )
                   !
                   ! ice volume, salt content, age content
@@ -321,5 +321,5 @@
                      t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
                      ! Snow energy of melting
-                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
+                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
                      !
                      ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
@@ -340,7 +340,7 @@
                      !
                      ! heat content per unit volume
-                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) )         &
-                        &             + lfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 )  )   &
-                        &             - rcp  * ( ztmelts - rt0 ) )
+                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoi * (   rcpi    * ( ztmelts - t_i(ji,jj,jk,jl) )           &
+                        &             + rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 )  )   &
+                        &             - rcp   * ( ztmelts - rt0 ) )
                      !
                      ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
@@ -410,5 +410,5 @@
       ! 5) Snow-ice mass (case ice is fully embedded)
       !----------------------------------------------
-      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhosn * v_s(:,:,:) + rhoic * v_i(:,:,:), dim=3  )   ! snow+ice mass
+      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass
       snwice_mass_b(:,:) = snwice_mass(:,:)
       !
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd.F90	(revision 9937)
@@ -295,7 +295,7 @@
             ztmelts       = -tmut * sz_i_1d(ji,jk)
             ! Conversion q(S,T) -> T (second order equation)
-            zbbb          = ( rcp - cpic ) * ztmelts + e_i_1d(ji,jk) * r1_rhoic - lfus
-            zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts, 0._wp ) )
-            t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_cpic
+            zbbb          = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus
+            zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) )
+            t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_cpi
             
             ! mask temperature
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_da.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_da.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_da.F90	(revision 9937)
@@ -137,5 +137,5 @@
             
             ! Contribution to salt flux
-            sfx_lam_1d(ji) = sfx_lam_1d(ji) + rhoic *  h_i_1d(ji) * zda * s_i_1d(ji) * r1_Dt_ice
+            sfx_lam_1d(ji) = sfx_lam_1d(ji) + rhoi *  h_i_1d(ji) * zda * s_i_1d(ji) * r1_Dt_ice
             
             ! Contribution to heat flux into the ocean [W.m-2], (<0)  
@@ -144,5 +144,5 @@
             
             ! Contribution to mass flux
-            wfx_lam_1d(ji) =  wfx_lam_1d(ji) + zda * r1_Dt_ice * ( rhoic * h_i_1d(ji) + rhosn * h_s_1d(ji) )
+            wfx_lam_1d(ji) =  wfx_lam_1d(ji) + zda * r1_Dt_ice * ( rhoi * h_i_1d(ji) + rhos * h_s_1d(ji) )
             
             ! new concentration
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_dh.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_dh.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_dh.F90	(revision 9937)
@@ -76,5 +76,5 @@
       REAL(wp) ::   zgrr         ! bottom growth rate
       REAL(wp) ::   zt_i_new     ! bottom formation temperature
-      REAL(wp) ::   z1_rho       ! 1/(rhosn+rho0-rhoic)
+      REAL(wp) ::   z1_rho       ! 1/(rhos+rho0-rhoi)
 
       REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
@@ -182,5 +182,5 @@
             IF( t_s_1d(ji,jk) > rt0 ) THEN
                hfx_res_1d    (ji) = hfx_res_1d    (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! heat flux to the ocean [W.m-2], < 0
-               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn         * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! mass flux
+               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos          * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! mass flux
                ! updates
                dh_s_mlt(ji)    = dh_s_mlt(ji) - zh_s(ji,jk)
@@ -202,9 +202,9 @@
             !
             ! --- precipitation ---
-            zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)   ! thickness change
+            zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji)    ! thickness change
             zqprec    (ji) = - qprec_ice_1d(ji)                                             ! enthalpy of the precip (>0, J.m-3)
             !
             hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji)    * r1_Dt_ice   ! heat flux from snow precip (>0, W.m-2)
-            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn         * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice   ! mass flux, <0
+            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos          * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice   ! mass flux, <0
             
             ! --- melt of falling snow ---
@@ -213,5 +213,5 @@
             zdeltah       (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) )                ! bound melting 
             hfx_snw_1d    (ji)   = hfx_snw_1d    (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji)    * r1_Dt_ice   ! heat used to melt snow (W.m-2, >0)
-            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn         * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip), >0
+            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhos          * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip), >0
             
             ! updates available heat + precipitations after melting
@@ -253,5 +253,5 @@
                
                hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_Dt_ice   ! heat used to melt snow(W.m-2, >0)
-               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn          * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip)
+               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip)
                
                ! updates available heat + thickness
@@ -273,6 +273,6 @@
          IF( evap_ice_1d(ji) > 0._wp ) THEN
             !
-            zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
-            zevap_rema(ji)   = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn   ! remaining evap in kg.m-2 (used for ice melting later on)
+            zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos  * rdt_ice )
+            zevap_rema(ji)   = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhos   ! remaining evap in kg.m-2 (used for ice melting later on)
             zdeltah   (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
             
@@ -280,5 +280,5 @@
                &                 ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) )  &
                &                 * a_i_1d(ji) * r1_Dt_ice
-            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_Dt_ice   ! Mass flux by sublimation
+            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_Dt_ice   ! Mass flux by sublimation
             
             ! new snow thickness
@@ -309,5 +309,5 @@
             e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *            &
               &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  &
-              &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
+              &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) )
          END DO
       END DO
@@ -326,24 +326,24 @@
             IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
 
-               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
+               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]       
                zdE            = 0._wp                                 ! Specific enthalpy difference (J/kg, <0)
                                                                       ! set up at 0 since no energy is needed to melt water...(it is already melted)
                zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
                                                                       ! this should normally not happen, but sometimes, heat diffusion leads to this
-               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0
+               zfmdt          = - zdeltah(ji,jk) * rhoi               ! Mass flux x time step > 0
                          
                dh_i_itm(ji)   = dh_i_itm(ji) + zdeltah(ji,jk)         ! Cumulate internal melting
                
-               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
+               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0]
 
                hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], <0
                !                                                                                                  ice enthalpy zEi is "sent" to the ocean
-               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux
+               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux
                !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok
-               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux
+               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux
 
             ELSE                                        !-- Surface melting
                
-               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
+               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]
                zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0]
                zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
@@ -351,22 +351,22 @@
                zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0]
                
-               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0]
+               zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0]
                
                zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0]
                
-               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
+               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat
                
                dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)         ! Cumulate surface melt
                
-               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0]
+               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0]
                
                zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
                
-               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux >0
+               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux >0
                !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok)
                hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice                           ! Heat flux [W.m-2], < 0
                hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat flux used in this process [W.m-2], > 0  
                ! 
-               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux
+               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux
                
             END IF
@@ -374,9 +374,9 @@
             ! Ice sublimation
             ! ---------------
-            zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic )
+            zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi )
             zdeltah (ji,jk) = zdeltah (ji,jk) + zdum
             dh_i_sub(ji)    = dh_i_sub(ji)    + zdum
             
-            sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0
+            sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0
             !                                                                                          clem: flux is sent to the ocean for simplicity
             !                                                                                                but salt should remain in the ice except
@@ -384,8 +384,8 @@
             hfx_sub_1d(ji)     = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_Dt_ice      ! Heat flux [W.m-2], < 0
 
-            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_Dt_ice          ! Mass flux > 0
+            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_Dt_ice          ! Mass flux > 0
 
             ! update remaining mass flux
-            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic
+            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoi
             
             ! record which layers have disappeared (for bottom melting) 
@@ -450,6 +450,6 @@
                zt_i_new      = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
                
-               zEi           = cpic * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0)
-                  &            - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts
+               zEi           = rcpi  * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0)
+                  &          - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp * ztmelts
 
                zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0)
@@ -457,19 +457,19 @@
                zdE           = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0)
 
-               dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
+               dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )
                
             END DO
             ! Contribution to Energy and Salt Fluxes                                    
-            zfmdt          = - rhoic * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0)
+            zfmdt          = - rhoi * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0)
             
             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], >0
             hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat flux used in this process [W.m-2], <0
             
-            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_Dt_ice    ! Salt flux, <0
-
-            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * r1_Dt_ice                  ! Mass flux, <0
+            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_Dt_ice    ! Salt flux, <0
+
+            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_Dt_ice                  ! Mass flux, <0
 
             ! update heat content (J.m-2) and layer thickness
-            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoic)
+            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi)
             h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji)
 
@@ -489,5 +489,5 @@
                IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting
 
-                  zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
+                  zEi               = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0)
                   zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
                                                                     !    set up at 0 since no energy is needed to melt water...(it is already melted)
@@ -497,11 +497,11 @@
                   dh_i_itm (ji)     = dh_i_itm(ji) + zdeltah(ji,jk)
 
-                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0
+                  zfmdt             = - zdeltah(ji,jk) * rhoi       ! Mass flux x time step > 0
 
                   hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], <0
                   !                                                                                                  ice enthalpy zEi is "sent" to the ocean
-                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux
+                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux
                   !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok
-                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux
+                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux
 
                   ! update heat content (J.m-2) and layer thickness
@@ -511,5 +511,5 @@
                ELSE                                        !-- Basal melting
 
-                  zEi             = - e_i_1d(ji,jk) * r1_rhoic                                ! Specific enthalpy of melting ice (J/kg, <0)
+                  zEi             = - e_i_1d(ji,jk) * r1_rhoi                                 ! Specific enthalpy of melting ice (J/kg, <0)
                   zEw             = rcp * ztmelts                                             ! Specific enthalpy of meltwater (J/kg, <0)
                   zdE             = zEi - zEw                                                 ! Specific enthalpy difference   (J/kg, <0)
@@ -517,13 +517,13 @@
                   zfmdt           = - zq_bo(ji) / zdE                                         ! Mass flux x time step (kg/m2, >0)
 
-                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic                                        ! Gross thickness change
+                  zdeltah(ji,jk)  = - zfmdt * r1_rhoi                                         ! Gross thickness change
 
                   zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change
                   
-                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE )   ! update available heat. MAX is necessary for roundup errors
-
-                  dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                            ! Update basal melt
-
-                  zfmdt           = - zdeltah(ji,jk) * rhoic                                  ! Mass flux x time step > 0
+                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoi * zdE )    ! update available heat. MAX is necessary for roundup errors
+
+                  dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                             ! Update basal melt
+
+                  zfmdt           = - zdeltah(ji,jk) * rhoi                                   ! Mass flux x time step > 0
 
                   zQm             = zfmdt * zEw                                               ! Heat exchanged with ocean
@@ -532,8 +532,8 @@
                   hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat used in this process [W.m-2], >0  
 
-                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice  ! Salt flux
+                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoi *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice  ! Salt flux
                   !                                                                                                   using s_i_1d and not sz_i_1d(jk) is ok
                   
-                  wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux
+                  wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux
 
                   ! update heat content (J.m-2) and layer thickness
@@ -566,5 +566,5 @@
          zq_rema(ji)        = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                               ! update available heat (J.m-2)
          hfx_snw_1d(ji)     = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_Dt_ice   ! Heat used to melt snow, W.m-2 (>0)
-         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice      ! Mass flux
+         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice      ! Mass flux
          dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
          !    
@@ -580,8 +580,8 @@
       ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level, 
       ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
-      z1_rho = 1._wp / ( rhosn+rho0-rhoic )
+      z1_rho = 1._wp / ( rhos + rho0 - rhoi )
       DO ji = 1, npti
          !
-         dh_snowice(ji) = MAX(  0._wp , ( rhosn * h_s_1d(ji) + (rhoic-rho0) * h_i_1d(ji) ) * z1_rho )
+         dh_snowice(ji) = MAX(  0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho )
 
          h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji)
@@ -589,5 +589,5 @@
 
          ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
-         zfmdt          = ( rhosn - rhoic ) * dh_snowice(ji)    ! <0
+         zfmdt          = ( rhos - rhoi ) * dh_snowice(ji)    ! <0
          zEw            = rcp * sst_1d(ji)
          zQm            = zfmdt * zEw 
@@ -599,11 +599,11 @@
          ! Case constant salinity in time: virtual salt flux to keep salinity constant
          IF( nn_icesal /= 2 )  THEN
-            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_Dt_ice  & ! put back sss_m     into the ocean
-               &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_Dt_ice    ! and get  rn_icesal from the ocean 
+            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                 * r1_Dt_ice  & ! put back sss_m     into the ocean
+               &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice    ! and get  rn_icesal from the ocean 
          ENDIF
 
          ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume
-         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_Dt_ice
-         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_Dt_ice
+         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice
+         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_Dt_ice
 
          ! update heat content (J.m-2) and layer thickness
@@ -627,5 +627,5 @@
             e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
             ! recalculate t_s_1d from e_s_1d
-            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhosn * r1_cpic + lfus * r1_cpic )
+            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_cpi + rLfus * r1_cpi )
          END DO
       END DO
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_do.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_do.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_do.F90	(revision 9937)
@@ -140,5 +140,5 @@
          ! Physical constants
          zhicrit = 0.04                                          ! frazil ice thickness
-         ztwogp  = 2. * rho0 / ( grav * 0.3 * ( rho0 - rhoic ) ) ! reduced grav
+         ztwogp  = 2. * rho0 / ( grav * 0.3 * ( rho0 - rhoi ) )  ! reduced grav
          zsqcd   = 1.0 / SQRT( 1.3 * zcai )                      ! 1/SQRT(airdensity*drag)
          zgamafr = 0.03
@@ -264,7 +264,7 @@
          DO ji = 1, npti
             ztmelts       = - tmut * zs_newice(ji)                  ! Melting point (C)
-            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
-               &                       + lfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
-               &                       - rcp  *         ztmelts )
+            ze_newice(ji) =   rhoi * (  rcpi  * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
+               &                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
+               &                      - rcp   *         ztmelts )
          END DO
 
@@ -275,5 +275,5 @@
          DO ji = 1, npti
 
-            zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [J/kg]
+            zEi           = - ze_newice(ji) * r1_rhoi              ! specific enthalpy of forming ice [J/kg]
 
             zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg]
@@ -284,5 +284,5 @@
             zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0) 
                                                                    ! clem: we use qlead instead of zqld (icethd) because we suppose we are at the freezing point   
-            zv_newice(ji) = - zfmdt * r1_rhoic
+            zv_newice(ji) = - zfmdt * r1_rhoi
 
             zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux  
@@ -293,7 +293,7 @@
             hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_Dt_ice
             ! mass flux
-            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_Dt_ice
+            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoi * r1_Dt_ice
             ! salt flux
-            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_Dt_ice
+            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_Dt_ice
          END DO
          
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_pnd.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_pnd.F90	(revision 9937)
@@ -133,5 +133,4 @@
       REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding
       REAL(wp) ::   z1_Tp            ! inverse reference temperature
-      REAL(wp) ::   z1_rhofw         ! inverse freshwater density
       REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
       REAL(wp) ::   zfac, zdum
@@ -139,5 +138,4 @@
       INTEGER  ::   ji   ! loop indices
       !!-------------------------------------------------------------------
-      z1_rhofw       = 1._wp / rhofw 
       z1_zpnd_aspect = 1._wp / zpnd_aspect
       z1_Tp          = 1._wp / zTp 
@@ -157,5 +155,5 @@
             !
             ! available meltwater for melt ponding [m, >0] and fraction
-            zdv_mlt = -( dh_i_sum(ji)*rhoic + dh_s_mlt(ji)*rhosn ) * z1_rhofw * a_i_1d(ji)
+            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * r1_rhow * a_i_1d(ji)
             zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc
             !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper 
@@ -168,5 +166,5 @@
             ! melt pond mass flux (<0)
             IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN
-               zfac = zfr_mlt * zdv_mlt * rhofw * r1_Dt_ice
+               zfac = zfr_mlt * zdv_mlt * rhow * r1_Dt_ice
                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
                !
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_sal.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_sal.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_sal.F90	(revision 9937)
@@ -77,5 +77,5 @@
             !---------------------------------------------------------
             IF( h_i_1d(ji) > 0._wp ) THEN
-               zs_sni   = sss_1d(ji) * ( rhoic - rhosn ) * r1_rhoic                 ! Salinity of snow ice
+               zs_sni   = sss_1d(ji) * ( rhoi - rhos ) * r1_rhoi                    ! Salinity of snow ice
                zs_i_si = ( zs_sni      - s_i_1d(ji) ) * dh_snowice(ji) / h_i_1d(ji) ! snow-ice    
                zs_i_bg = ( s_i_new(ji) - s_i_1d(ji) ) * dh_i_bog  (ji) / h_i_1d(ji) ! bottom growth
@@ -98,5 +98,5 @@
                
                ! Salt flux
-               sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoic * a_i_1d(ji) * h_i_1d(ji) * ( zs_i_fl + zs_i_gd ) * r1_Dt_ice
+               sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoi * a_i_1d(ji) * h_i_1d(ji) * ( zs_i_fl + zs_i_gd ) * r1_Dt_ice
             ENDIF
          END DO
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_zdf_bl99.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_zdf_bl99.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_zdf_bl99.F90	(revision 9937)
@@ -217,11 +217,16 @@
             !
             DO ji = 1, npti
-               ztcond_i(ji,0)      = rcdic + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
-               ztcond_i(ji,nlay_i) = rcdic + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
+               ztcond_i(ji,0)      = rcnd_i + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
+               ztcond_i(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
             END DO
             DO jk = 1, nlay_i-1
                DO ji = 1, npti
-                  ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  &
-                     &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
+!!gm faster coding
+                  ztcond_i(ji,jk) = rcnd_i + zbeta * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) )             &
+                     &                          / MIN( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) - 2._wp * rt0 , -epsi10 )
+!!gm old
+!                  ztcond_i(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) )   &
+!                     &                             / MIN( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 , -epsi10 )
+!!gm 
                END DO
             END DO
@@ -230,14 +235,20 @@
             !
             DO ji = 1, npti
-               ztcond_i(ji,0)      = rcdic + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
-                  &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
-               ztcond_i(ji,nlay_i) = rcdic + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  &
-                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )
+               ztcond_i(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
+                  &                         - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
+               ztcond_i(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  &
+                  &                         - 0.011_wp * ( t_bo_1d(ji) - rt0 )
             END DO
             DO jk = 1, nlay_i-1
                DO ji = 1, npti
-                  ztcond_i(ji,jk) = rcdic + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        &
-                     &                     MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  &
-                     &                    - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )
+!!gm faster coding
+                  zfac   = t_i_1d (ji,jk) + t_i_1d (ji,jk+1) - 2._wp * rt0
+                  ztcond_i(ji,jk) = rcnd_i + 0.09_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / MIN( -epsi10, zfac )  &
+                     &                     - 0.0055_wp * zfac            !NB:  0.0055 = 1/2 * 0.011
+!!gm old
+!                  ztcond_i(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        &
+!                     &                      MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  &
+!                     &                     - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )
+!!gm 
                END DO
             END DO
@@ -299,6 +310,6 @@
          DO jk = 1, nlay_i
             DO ji = 1, npti
-               zcpi = cpic + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
-               zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) 
+               zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
+               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 
             END DO
          END DO
@@ -306,5 +317,5 @@
          DO jk = 1, nlay_s
             DO ji = 1, npti
-               zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)
+               zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_cpi * z1_h_s(ji)
             END DO
          END DO
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/iceupdate.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/iceupdate.F90	(revision 9937)
@@ -172,5 +172,5 @@
             snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step
             !                                               ! new mass per unit area
-            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj)  ) 
+            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj)  ) 
             !                                               ! time evolution of snow+ice mass
             snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_Dt_ice
@@ -432,10 +432,10 @@
             ELSE                                     ! start from rest
                IF(lwp) WRITE(numout,*) '   ==>>   previous run without snow-ice mass output then set it'
-               snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
+               snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) )
                snwice_mass_b(:,:) = snwice_mass(:,:)
             ENDIF
          ELSE                                   !* Start from rest
             IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set the snow-ice mass'
-            snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
+            snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) )
             snwice_mass_b(:,:) = snwice_mass(:,:)
          ENDIF
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icevar.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icevar.F90	(revision 9937)
@@ -228,7 +228,7 @@
                      ztmelts          = - sz_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
                      ! Conversion q(S,T) -> T (second order equation)
-                     zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
-                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
-                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_cpic , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts
+                     zbbb             = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus
+                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) )
+                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_cpi , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts
                      !
                   ELSE                                !--- no ice
@@ -247,5 +247,5 @@
          WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
             t_s(:,:,jk,:) = rt0 + MAX( -100._wp ,  &
-                 &                MIN( r1_cpic * ( -r1_rhosn * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + lfus ) , 0._wp ) )
+                 &                MIN( r1_cpi * ( -r1_rhos * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + rLfus ) , 0._wp ) )
          ELSEWHERE                           !--- no ice
             t_s(:,:,jk,:) = rt0
@@ -498,7 +498,7 @@
             DO ji = 1 , jpi
                ! update exchanges with ocean
-               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoic * r1_Dt_ice
-               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoic * r1_Dt_ice
-               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_Dt_ice
+               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_Dt_ice
+               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_Dt_ice
+               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_Dt_ice
                !
                !-----------------------------------------------------------------
@@ -669,8 +669,8 @@
                ! In case snow load is in excess that would lead to transformation from snow to ice
                ! Then, transfer the snow excess into the ice (different from icethd_dh)
-               zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rho0 ) * zh_i(ji,jl) ) * r1_rho0 ) 
+               zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rho0 ) * zh_i(ji,jl) ) * r1_rho0 ) 
                ! recompute h_i, h_s avoiding out of bounds values
                zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
-               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn )
+               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )
             ENDIF
          END DO
@@ -854,13 +854,13 @@
             ztmelts      = - tmut  * sz_i_1d(ji,jk)
             t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
-                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
-            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
-               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
-               &                    - rcp  *   ztmelts )
+            !                                                   !   (sometimes zdf scheme produces abnormally high temperatures)   
+            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
+               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
+               &                   - rcp   *   ztmelts )
          END DO
       END DO
       DO jk = 1, nlay_s             ! Snow energy of melting
          DO ji = 1, npti
-            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
+            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
          END DO
       END DO
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icewri.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icewri.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icewri.F90	(revision 9937)
@@ -85,5 +85,5 @@
       ! Standard outputs
       !-----------------
-      zrho1 = ( rho0 - rhoic ) * r1_rho0   ;   zrho2 = rhosn * r1_rho0
+      zrho1 = ( rho0 - rhoi ) * r1_rho0   ;   zrho2 = rhos * r1_rho0
       ! masks
       IF( iom_use('icemask'  ) )   CALL iom_put( "icemask"  , zmsk00              )   ! ice mask 0%
@@ -92,6 +92,6 @@
       !
       ! general fields
-      IF( iom_use('icemass'  ) )   CALL iom_put( "icemass", rhoic * vt_i * zmsk00 )   ! Ice mass per cell area 
-      IF( iom_use('snwmass'  ) )   CALL iom_put( "snwmass", rhosn * vt_s * zmsksn )   ! Snow mass per cell area
+      IF( iom_use('icemass'  ) )   CALL iom_put( "icemass", rhoi * vt_i * zmsk00  )   ! Ice mass per cell area 
+      IF( iom_use('snwmass'  ) )   CALL iom_put( "snwmass", rhos * vt_s * zmsksn  )   ! Snow mass per cell area
       IF( iom_use('icepres'  ) )   CALL iom_put( "icepres", zmsk00                )   ! Ice presence (1 or 0) 
       IF( iom_use('iceconc'  ) )   CALL iom_put( "iceconc", at_i  * zmsk00        )   ! ice concentration
@@ -104,4 +104,7 @@
       IF( iom_use('snwvolu'  ) )   CALL iom_put( "snwvolu", vt_s  * zmsksn        )   ! snow volume
       IF( iom_use('icefrb') ) THEN
+!!gm remove the WHERE by using :
+!!         z2d(:,:) = MAX( zrho1 * hm_i(:,:) - zrho2 * hm_s(:,:) , 0._wp )                                         
+!!gm end
          z2d(:,:) = ( zrho1 * hm_i(:,:) - zrho2 * hm_s(:,:) )                                         
          WHERE( z2d < 0._wp )   z2d = 0._wp
@@ -115,5 +118,5 @@
       ! salt
       IF( iom_use('icesalt'  ) )   CALL iom_put( "icesalt", sm_i  * zmsk00        )   ! mean ice salinity
-      IF( iom_use('icesalm'  ) )   CALL iom_put( "icesalm", SUM( sv_i, DIM = 3 ) * rhoic * 1.0e-3 * zmsk00 )   ! Mass of salt in sea ice per cell area
+      IF( iom_use('icesalm'  ) )   CALL iom_put( "icesalm", SUM( sv_i, DIM = 3 ) * rhoi * 1.0e-3 * zmsk00 )   ! Mass of salt in sea ice per cell area
 
       ! heat
@@ -164,5 +167,5 @@
       ! trends
       IF( iom_use('dmithd') )   CALL iom_put( "dmithd", - wfx_bog - wfx_bom - wfx_sum - wfx_sni - wfx_opw - wfx_lam - wfx_res ) ! Sea-ice mass change from thermodynamics
-      IF( iom_use('dmidyn') )   CALL iom_put( "dmidyn", - wfx_dyn + rhoic * diag_trp_vi     )   ! Sea-ice mass change from dynamics(kg/m2/s)
+      IF( iom_use('dmidyn') )   CALL iom_put( "dmidyn", - wfx_dyn + rhoi * diag_trp_vi      )   ! Sea-ice mass change from dynamics(kg/m2/s)
       IF( iom_use('dmiopw') )   CALL iom_put( "dmiopw", - wfx_opw                           )   ! Sea-ice mass change through growth in open water
       IF( iom_use('dmibog') )   CALL iom_put( "dmibog", - wfx_bog                           )   ! Sea-ice mass change through basal growth
@@ -174,7 +177,7 @@
       IF( iom_use('dmisub') )   CALL iom_put( "dmisub", - wfx_ice_sub                       )   ! Sea-ice mass change through sublimation
       IF( iom_use('dmsspr') )   CALL iom_put( "dmsspr", - wfx_spr                           )   ! Snow mass change through snow fall
-      IF( iom_use('dmsssi') )   CALL iom_put( "dmsssi",   wfx_sni*rhosn*r1_rhoic            )   ! Snow mass change through snow-to-ice conversion
+      IF( iom_use('dmsssi') )   CALL iom_put( "dmsssi",   wfx_sni*rhos*r1_rhoi              )   ! Snow mass change through snow-to-ice conversion
       IF( iom_use('dmsmel') )   CALL iom_put( "dmsmel", - wfx_snw_sum                       )   ! Snow mass change through melt
-      IF( iom_use('dmsdyn') )   CALL iom_put( "dmsdyn", - wfx_snw_dyn + rhosn * diag_trp_vs )   ! Snow mass change through dynamics(kg/m2/s)
+      IF( iom_use('dmsdyn') )   CALL iom_put( "dmsdyn", - wfx_snw_dyn + rhos * diag_trp_vs  )   ! Snow mass change through dynamics(kg/m2/s)
 
       ! Global ice diagnostics
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/BDY/bdyice.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/BDY/bdyice.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/BDY/bdyice.F90	(revision 9937)
@@ -124,17 +124,17 @@
 
             ! Then, a) transfer the snow excess into the ice (different from icethd_dh)
-            zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 )
+            zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 )
             ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment)
-            !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhosn / rhoic )
+            !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos / rhoi )
 
             ! recompute h_i, h_s
             h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh )
-            h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic / rhosn ) 
-
-         ENDDO
+            h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos ) 
+
+         END DO
          CALL lbc_bdy_lnk( a_i(:,:,jl), 'T', 1., ib_bdy )
          CALL lbc_bdy_lnk( h_i(:,:,jl), 'T', 1., ib_bdy )
          CALL lbc_bdy_lnk( h_s(:,:,jl), 'T', 1., ib_bdy )
-      ENDDO
+      END DO
       ! retrieve at_i
       at_i(:,:) = 0._wp
@@ -212,5 +212,5 @@
             DO jk = 1, nlay_s
                ! Snow energy of melting
-               e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
+               e_s(ji,jj,jk,jl) = rswitch * rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
                ! Multiply by volume, so that heat content in J/m2
                e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
@@ -219,8 +219,8 @@
                ztmelts          = - tmut * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K                  
                ! heat content per unit volume
-               e_i(ji,jj,jk,jl) = rswitch * rhoic * &
-                  (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
-                  +   lfus    * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
-                  - rcp      * ( ztmelts - rt0 ) )
+               e_i(ji,jj,jk,jl) = rswitch * rhoi * &
+                  (   rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) &
+                  +   rLfus * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
+                  -   rcp   * ( ztmelts - rt0 ) )
                ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
                e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * h_i(ji,jj,jl) * r1_nlay_i
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/phycst.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/phycst.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/phycst.F90	(revision 9937)
@@ -34,61 +34,37 @@
    REAL(wp), PUBLIC ::   rhhmm =  60._wp        !: number of minutes in one hour
    REAL(wp), PUBLIC ::   rmmss =  60._wp        !: number of seconds in one minute
-   REAL(wp), PUBLIC ::   omega                  !: earth rotation parameter           [s-1]
-   REAL(wp), PUBLIC ::   ra    = 6371229._wp    !: earth radius                       [m]
-   REAL(wp), PUBLIC ::   grav  = 9.80665_wp     !: gravity                            [m/s2]
+   REAL(wp), PUBLIC ::   omega                  !: earth rotation parameter               [s-1]
+   REAL(wp), PUBLIC ::   ra    = 6371229._wp    !: earth radius                           [m]
+   REAL(wp), PUBLIC ::   grav  = 9.80665_wp     !: gravity                                [m/s2]
    
-   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_si3
-   REAL(wp), PUBLIC ::   rt0_ice  = 273.15_wp        !: melting point of ice          [Kelvin]
-#else
-   REAL(wp), PUBLIC ::   rt0_ice  = 273.05_wp        !: melting point of ice          [Kelvin]
-#endif
-   REAL(wp), PUBLIC ::   rho0                        !: volumic mass of reference     [kg/m3]
-   REAL(wp), PUBLIC ::   r1_rho0                     !: = 1. / rho0                   [m3/kg]
-   REAL(wp), PUBLIC ::   rcp                         !: ocean specific heat           [J/Kelvin]
-   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [Kelvin/J]
+   REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water     [Kelvin]
+   REAL(wp), PUBLIC ::   rho0                        !: volumic mass of reference         [kg/m3]
+   REAL(wp), PUBLIC ::   r1_rho0                     !: = 1. / rho0                       [m3/kg]
+   REAL(wp), PUBLIC ::   rcp                         !: ocean specific heat               [J/Kelvin]
+   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                        [Kelvin/J]
    REAL(wp), PUBLIC ::   rho0_rcp                    !: = rho0 * rcp 
    REAL(wp), PUBLIC ::   r1_rho0_rcp                 !: = 1. / ( rho0 * rcp )
 
-   REAL(wp), PUBLIC ::   rhosn    =  330._wp         !: volumic mass of snow          [kg/m3]
-   REAL(wp), PUBLIC ::   rhofw    = 1000._wp         !: volumic mass of freshwater in melt ponds [kg/m3]
+   REAL(wp), PUBLIC ::   rhoi     =  917._wp         !: sea ice density                   [kg/m3]
+   REAL(wp), PUBLIC ::   rhos     =  330._wp         !: snow    density                   [kg/m3]
+   REAL(wp), PUBLIC ::   rhow     = 1000._wp         !: water   density (in melt ponds)   [kg/m3]
+   REAL(wp), PUBLIC ::   rcnd_i   =    2.034396_wp   !: thermal conductivity of fresh ice [W/m/K]
+   REAL(wp), PUBLIC ::   rcpi     = 2067.0_wp        !: specific heat of fresh ice        [J/kg/K]
+   REAL(wp), PUBLIC ::   rLsub    =    2.834e+6_wp   !: pure ice latent heat of sublimation   [J/kg]
+   REAL(wp), PUBLIC ::   rLfus    =    0.334e+6_wp   !: latent heat of fusion of fresh ice    [J/kg]
+   REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity
    REAL(wp), PUBLIC ::   emic     =    0.97_wp       !: emissivity of snow or ice
-   REAL(wp), PUBLIC ::   sice     =    6.0_wp        !: salinity of ice               [psu]
-   REAL(wp), PUBLIC ::   soce     =   34.7_wp        !: salinity of sea               [psu]
-   REAL(wp), PUBLIC ::   cevap    =    2.5e+6_wp     !: latent heat of evaporation (water)
-   REAL(wp), PUBLIC ::   srgamma  =    0.9_wp        !: correction factor for solar radiation (Oberhuber, 1974)
+   REAL(wp), PUBLIC ::   sice     =    6.0_wp        !: salinity of ice                   [psu]
+   REAL(wp), PUBLIC ::   soce     =   34.7_wp        !: salinity of sea                   [psu]
    REAL(wp), PUBLIC ::   vkarmn   =    0.4_wp        !: von Karman constant
    REAL(wp), PUBLIC ::   stefan   =    5.67e-8_wp    !: Stefan-Boltzmann constant 
 
-#if defined key_si3 || 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                     [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]
-   REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity
-   REAL(wp), PUBLIC ::   xlsn                        !: = lfus*rhosn (volumetric latent heat fusion of snow)  [J/m3]
-#else
-   REAL(wp), PUBLIC ::   rhoic    =  900._wp         !: volumic mass of sea ice                               [kg/m3]
-   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: conductivity of the ice                               [W/m/K]
-   REAL(wp), PUBLIC ::   rcpic    =    1.8837e+6_wp  !: volumetric specific heat for ice                      [J/m3/K]
-   REAL(wp), PUBLIC ::   cpic                        !: = rcpic / rhoic  (specific heat for ice)              [J/Kg/K]
-   REAL(wp), PUBLIC ::   rcdsn    =    0.22_wp       !: conductivity of the snow                              [W/m/K]
-   REAL(wp), PUBLIC ::   rcpsn    =    6.9069e+5_wp  !: volumetric specific heat for snow                     [J/m3/K]
-   REAL(wp), PUBLIC ::   xlsn     =  110.121e+6_wp   !: volumetric latent heat fusion of snow                 [J/m3]
-   REAL(wp), PUBLIC ::   lfus                        !: = xlsn / rhosn   (latent heat of fusion of fresh ice) [J/Kg]
-   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]
-#endif
-#if defined key_si3
-   REAL(wp), PUBLIC ::   r1_rhoic                    !: 1 / rhoic
-   REAL(wp), PUBLIC ::   r1_rhosn                    !: 1 / rhosn
-   REAL(wp), PUBLIC ::   r1_cpic                     !: 1 / cpic
-#endif
+   REAL(wp), PUBLIC ::   r1_rhoi                     !: 1 / rhoi
+   REAL(wp), PUBLIC ::   r1_rhos                     !: 1 / rhos
+   REAL(wp), PUBLIC ::   r1_rhow                     !: 1 / rhow
+   REAL(wp), PUBLIC ::   r1_cpi                      !: 1 / rcpi
+   REAL(wp), PUBLIC ::   r1_Lsub                     !: 1 / rLsub
+   REAL(wp), PUBLIC ::   r1_Lfus                     !: 1 / rLfus
+
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -105,16 +81,10 @@
       !! ** Purpose :   set and print the constants
       !!----------------------------------------------------------------------
-
+      !
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'phy_cst : initialization of ocean parameters and constants'
       IF(lwp) WRITE(numout,*) '~~~~~~~'
 
-      ! Define & print constants
-      ! ------------------------
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '   Constants'
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '      mathematical constant                 rpi = ', rpi
+      !                 !==  Define derived constant  ==!
 
       rsiyea = 365.25_wp * rday * 2._wp * rpi / 6.283076_wp
@@ -125,69 +95,50 @@
       omega  = 2._wp * rpi / rsiday 
 #endif
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '      day                                rday   = ', rday,   ' s'
-      IF(lwp) WRITE(numout,*) '      sideral year                       rsiyea = ', rsiyea, ' s'
-      IF(lwp) WRITE(numout,*) '      sideral day                        rsiday = ', rsiday, ' s'
-      IF(lwp) WRITE(numout,*) '      omega                              omega  = ', omega,  ' s^-1'
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '      nb of months per year               raamo = ', raamo, ' months'
-      IF(lwp) WRITE(numout,*) '      nb of hours per day                 rjjhh = ', rjjhh, ' hours'
-      IF(lwp) WRITE(numout,*) '      nb of minutes per hour              rhhmm = ', rhhmm, ' mn'
-      IF(lwp) WRITE(numout,*) '      nb of seconds per minute            rmmss = ', rmmss, ' s'
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '      earth radius                         ra   = ', ra, ' m'
-      IF(lwp) WRITE(numout,*) '      gravity                              grav = ', grav , ' m/s^2'
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '      triple point of temperature      rtt      = ', rtt     , ' K'
-      IF(lwp) WRITE(numout,*) '      freezing point of water          rt0      = ', rt0     , ' K'
-      IF(lwp) WRITE(numout,*) '      melting point of snow            rt0_snow = ', rt0_snow, ' K'
-      IF(lwp) WRITE(numout,*) '      melting point of ice             rt0_ice  = ', rt0_ice , ' K'
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '   reference density and heat capacity now defined in eosbn2.f90'
-              
-#if defined key_si3 || defined key_cice
-      xlsn = lfus * rhosn        ! volumetric latent heat fusion of snow [J/m3]
-#else
-      cpic = rcpic / rhoic       ! specific heat for ice   [J/Kg/K]
-      lfus = xlsn / rhosn        ! latent heat of fusion of fresh ice
-#endif
-#if defined key_si3
-      r1_rhoic = 1._wp / rhoic
-      r1_rhosn = 1._wp / rhosn
-      r1_cpic  = 1._wp / cpic
-#endif
-      IF(lwp) THEN
+
+      r1_rhoi = 1._wp / rhoi
+      r1_rhos = 1._wp / rhos
+      r1_cpi  = 1._wp / rcpi
+      r1_Lsub = 1._wp / rLsub
+      r1_Lfus = 1._wp / rLfus
+
+      IF(lwp) THEN      !==  print constants  ==!
          WRITE(numout,*)
-#if defined key_cice
-         WRITE(numout,*) '      thermal conductivity of the snow          = ', rcdsn   , ' 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'
-#if defined key_si3 || defined key_cice
-         WRITE(numout,*) '      latent heat of subl.  of fresh ice / snow = ', lsub    , ' J/kg'
-#else
-         WRITE(numout,*) '      density times specific heat for snow      = ', rcpsn   , ' J/m^3/K' 
-         WRITE(numout,*) '      density times specific heat for ice       = ', rcpic   , ' J/m^3/K'
-         WRITE(numout,*) '      volumetric latent heat fusion of sea ice  = ', xlic    , ' J/m' 
-         WRITE(numout,*) '      latent heat of sublimation of snow        = ', xsn     , ' J/kg' 
-#endif
-         WRITE(numout,*) '      volumetric latent heat fusion of snow     = ', xlsn    , ' J/m^3' 
-         WRITE(numout,*) '      density of sea ice                        = ', rhoic   , ' kg/m^3'
-         WRITE(numout,*) '      density of snow                           = ', rhosn   , ' kg/m^3'
-         WRITE(numout,*) '      density of freshwater (in melt ponds)     = ', rhofw   , ' kg/m^3'
-         WRITE(numout,*) '      emissivity of snow or ice                 = ', emic  
-         WRITE(numout,*) '      salinity of ice                           = ', sice    , ' psu'
-         WRITE(numout,*) '      salinity of sea                           = ', soce    , ' psu'
-         WRITE(numout,*) '      latent heat of evaporation (water)        = ', cevap   , ' J/m^3' 
-         WRITE(numout,*) '      correction factor for solar radiation     = ', srgamma 
-         WRITE(numout,*) '      von Karman constant                       = ', vkarmn 
-         WRITE(numout,*) '      Stefan-Boltzmann constant                 = ', stefan  , ' J/s/m^2/K^4'
+         WRITE(numout,*) '   Constants'
          WRITE(numout,*)
-         WRITE(numout,*) '      conversion: degre ==> radian          rad = ', rad
+         WRITE(numout,*) '      mathematical constant              rpi    = ', rpi
+         WRITE(numout,*) '      conversion: degre ==> radian       rad    = ', rad
          WRITE(numout,*)
-         WRITE(numout,*) '      smallest real computer value       rsmall = ', rsmall
+         WRITE(numout,*) '      day in seconds                     rday   = ', rday  , ' s'
+         WRITE(numout,*) '      sideral year                       rsiyea = ', rsiyea, ' s'
+         WRITE(numout,*) '      sideral day                        rsiday = ', rsiday, ' s'
+         WRITE(numout,*) '      omega = 2 pi / rsiday              omega  = ', omega , ' s^-1'
+         WRITE(numout,*) '      earth radius                       ra     = ', ra    , ' m'
+         WRITE(numout,*) '      gravity                            grav   = ', grav  , ' m/s^2'
+         WRITE(numout,*)
+         WRITE(numout,*) '      nb of months per year              raamo  = ', raamo, ' months'
+         WRITE(numout,*) '      nb of hours per day                rjjhh  = ', rjjhh, ' hours'
+         WRITE(numout,*) '      nb of minutes per hour             rhhmm  = ', rhhmm, ' mn'
+         WRITE(numout,*) '      nb of seconds per minute           rmmss  = ', rmmss, ' s'
+         WRITE(numout,*)
+         WRITE(numout,*) '   reference ocean density and heat capacity now defined in eosbn2.f90'
+         WRITE(numout,*)
+         WRITE(numout,*) '      freezing point of freshwater                rt0    = ', rt0   , ' K'
+         WRITE(numout,*) '      sea ice density                             rhoi   = ', rhoi  , ' kg/m^3'
+         WRITE(numout,*) '      snow    density                             rhos   = ', rhos  , ' kg/m^3'
+         WRITE(numout,*) '      freshwater density (in melt ponds)          rhow   = ', rhow  , ' kg/m^3'
+         WRITE(numout,*) '      thermal conductivity of pure ice            rcnd_i = ', rcnd_i, ' J/s/m/K'
+         WRITE(numout,*) '      fresh ice specific heat                     rcpi   = ', rcpi  , ' J/kg/K'
+         WRITE(numout,*) '      latent heat of fusion of fresh ice / snow   rLfus  = ', rLfus , ' J/kg'
+         WRITE(numout,*) '      latent heat of subl.  of fresh ice / snow   rLsub  = ', rLsub , ' J/kg'
+         WRITE(numout,*) '      emissivity of snow or ice                   emic   = ', emic  
+         WRITE(numout,*) '      salinity of ice                             sice   = ', sice  , ' psu'
+         WRITE(numout,*) '      salinity of sea                             soce   = ', soce  , ' psu'
+         WRITE(numout,*) '      von Karman constant                         vkarmn = ', vkarmn 
+         WRITE(numout,*) '      Stefan-Boltzmann constant                   stefan = ', stefan, ' J/s/m^2/K^4'
+         WRITE(numout,*)
+         WRITE(numout,*)
+         WRITE(numout,*) '      smallest real computer value                rsmall = ', rsmall
       ENDIF
-
+      !
    END SUBROUTINE phy_cst
 
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcblk.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcblk.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcblk.F90	(revision 9937)
@@ -504,5 +504,5 @@
       ENDIF
 
-      zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux
+      zqla(:,:) = L_vap( zst(:,:) ) * zevap(:,:)     ! Latent Heat flux
 
 
@@ -526,10 +526,10 @@
       !
       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
-         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
+         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
          &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
          &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
          &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
          &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
-         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic
+         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
       qns(:,:) = qns(:,:) * tmask(:,:,1)
       !
@@ -643,5 +643,4 @@
       !! ** Purpose : Compute the moist adiabatic lapse-rate.
       !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
-      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
       !!
       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
@@ -652,5 +651,6 @@
       !
       INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) :: zrv, ziRT        ! local scalar
+      REAL(wp) ::   zrv, ziRT      ! local scalar
+      REAL(wp) ::   zLv = 2.5e+6_wp   ! latent heat of vaporisation 
       !!----------------------------------------------------------------------------------
       !
@@ -659,5 +659,5 @@
             zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
             ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
-            gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) )
+            gamma_moist(ji,jj) = grav * ( 1. + zLv*zrv*ziRT ) / ( Cp_dry + zLv*zLv*zrv*reps0*ziRT/ptak(ji,jj) )
          END DO
       END DO
@@ -792,5 +792,5 @@
       REAL(wp) ::   zst3                     ! local variable
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
-      REAL(wp) ::   zztmp, z1_lsub           !   -      -
+      REAL(wp) ::   zztmp                    !   -      -
       REAL(wp) ::   zfr1, zfr2               ! local variables
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
@@ -868,7 +868,6 @@
 
       ! --- evaporation --- !
-      z1_lsub = 1._wp / Lsub
-      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub    ! sublimation
-      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub    ! d(sublimation)/dT
+      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * r1_Lsub    ! sublimation
+      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * r1_Lsub    ! d(sublimation)/dT
       zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )  ! evaporation over ocean
 
@@ -884,7 +883,7 @@
          &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
          &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
-         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
       qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
-         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
 
       ! --- total solar and non solar fluxes --- !
@@ -894,9 +893,9 @@
 
       ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
-      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
 
       ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
       DO jl = 1, jpl
-         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) )
+         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) )
          !                         ! But we do not have Tice => consider it at 0degC => evap=0 
       END DO
@@ -971,5 +970,5 @@
       CASE ( 1 , 2 )
          !
-         zfac  = 1._wp /  ( rn_cnd_s + rcdic )
+         zfac  = 1._wp /  ( rn_cnd_s + rcnd_i )
          zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
          zfac3 = 2._wp / zepsilon
@@ -978,6 +977,6 @@
             DO jj = 1 , jpj
                DO ji = 1, jpi
-                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac                             ! Effective thickness
-                  IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
+                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                             ! Effective thickness
+                  IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) )  ! Enhanced conduction factor
                END DO
             END DO
@@ -990,12 +989,11 @@
       ! -------------------------------------------------------------!
       !
-      zfac = rcdic * rn_cnd_s
+      zfac = rcnd_i * rn_cnd_s
       !
       DO jl = 1, jpl
          DO jj = 1 , jpj
             DO ji = 1, jpi
-               !                    
-               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
-                  &      ( rcdic * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
+               !                                                                       ! Effective conductivity of the snow-ice system divided by thickness
+               zkeff_h = zfac * zgfac(ji,jj,jl) / ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
                ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature
                ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbccpl.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbccpl.F90	(revision 9937)
@@ -1418,9 +1418,9 @@
             zqns(:,:) =  zqns(:,:) - zemp(:,:) * sst_m(:,:) * rcp         ! remove heat content due to mass flux (assumed to be at SST)
             IF( srcv(jpr_snow  )%laction ) THEN
-               zqns(:,:) = zqns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus    ! energy for melting solid precipitation over the free ocean
+               zqns(:,:) = zqns(:,:) - frcv(jpr_snow)%z3(:,:,1) * rLfus   ! energy for melting solid precipitation over the free ocean
             ENDIF
          ENDIF
          !
-         IF( srcv(jpr_icb)%laction )  zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus ! remove heat content associated to iceberg melting
+         IF( srcv(jpr_icb)%laction )  zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * rLfus ! remove heat content associated to iceberg melting
          !
          IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
@@ -1811,8 +1811,8 @@
       !                                     
       ! --- calving (removed from qns_tot) --- !
-      IF( srcv(jpr_cal)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * lfus  ! remove latent heat of calving
-                                                                                                    ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean
+      IF( srcv(jpr_cal)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * rLfus  ! remove latent heat of calving
+                                                                                                     ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean
       ! --- iceberg (removed from qns_tot) --- !
-      IF( srcv(jpr_icb)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus  ! remove latent heat of iceberg melting
+      IF( srcv(jpr_icb)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * rLfus  ! remove latent heat of iceberg melting
 
 #if defined key_si3      
@@ -1823,16 +1823,16 @@
 
       ! 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(:,:) = rcpi * SUM( (tn_ice - rt0) * a_i, dim=3 ) / SUM( a_i, dim=3 )
       ELSEWHERE                             ;   zcptsnw(:,:) = zcptn(:,:)
-      ENDWHERE
+      END WHERE
       ! Heat content per unit mass of rain (J/kg)
       zcptrain(:,:) = rcp * ( SUM( (tn_ice(:,:,:) - rt0) * a_i(:,:,:), dim=3 ) + sst_m(:,:) * ziceld(:,:) ) 
 
       ! --- enthalpy of snow precip over ice in J/m3 (to be used in 1D-thermo) --- !
-      zqprec_ice(:,:) = rhosn * ( zcptsnw(:,:) - lfus )
+      zqprec_ice(:,:) = rhos * ( zcptsnw(:,:) - rLfus )
 
       ! --- 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 atm. does not take it into account
+         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * rcpi ) but atm. does not take it into account
       END DO
 
@@ -1840,8 +1840,8 @@
       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(:,:) =     zsprecip(:,:)                   * zsnw             * ( zcptsnw (:,:) - lfus )   ! solid precip over ice (qevap_ice=0 since atm. does not take it into account)
+         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptsnw (:,:) - rLfus )   ! solid precip over ocean + snow melting
+      zqemp_ice(:,:) =     zsprecip(:,:)                   * zsnw             * ( zcptsnw (:,:) - rLfus )   ! solid precip over ice (qevap_ice=0 since atm. does not take it into account)
 !!    zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * picefr(:,:)      *   zcptsnw (:,:)   &        ! ice evap
-!!       &             +   zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhosn ! solid precip over ice
+!!       &             +   zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhos  ! solid precip over ice
       
       ! --- total non solar flux (including evap/precip) --- !
@@ -1875,5 +1875,5 @@
       ! clem: this formulation is certainly wrong... but better than it was...
       zqns_tot(:,:) = zqns_tot(:,:)                            &          ! zqns_tot update over free ocean with:
-         &          - (  ziceld(:,:) * zsprecip(:,:) * lfus )  &          ! remove the latent heat flux of solid precip. melting
+         &          - (  ziceld(:,:) * zsprecip(:,:) * rLfus )  &         ! remove the latent heat flux of solid precip. melting
          &          - (  zemp_tot(:,:)                         &          ! remove the heat content of mass flux (assumed to be at SST)
          &             - zemp_ice(:,:) ) * zcptn(:,:) 
@@ -1892,13 +1892,13 @@
 #endif
       ! outputs
-      IF ( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * lfus )                       ! latent heat from calving
-      IF ( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * lfus )                       ! latent heat from icebergs melting
+      IF ( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )                      ! latent heat from calving
+      IF ( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )                      ! latent heat from icebergs melting
       IF ( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'   , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average)
       IF ( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) &
            &                                                              * picefr(:,:) ) * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average)
-      IF ( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - Lfus )   )               ! heat flux from snow (cell average)
-      IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - Lfus ) &
+      IF ( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )   )               ! heat flux from snow (cell average)
+      IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &
            &                                                              * ( 1._wp - zsnw(:,:) )                  )               ! heat flux from snow (over ocean)
-      IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - Lfus ) &
+      IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &
            &                                                              *           zsnw(:,:)                    )               ! heat flux from snow (over ice)
       ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcice_cice.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcice_cice.F90	(revision 9937)
@@ -13,5 +13,5 @@
    USE dom_oce         ! ocean space and time domain
    USE domvvl
-   USE phycst   , ONLY : rcp, rho0, r1_rho0, rhosn, rhoic
+   USE phycst   , ONLY : rcp, rho0, r1_rho0, rhos, rhoi
    USE in_out_manager  ! I/O manager
    USE iom, ONLY : iom_put,iom_use              ! I/O manager library !!Joakim edit
@@ -222,5 +222,5 @@
       CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
       CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
-      snwice_mass  (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:)  )
+      snwice_mass  (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:)  )
       snwice_mass_b(:,:) = snwice_mass(:,:)
 
@@ -644,5 +644,5 @@
       CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
       CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
-      snwice_mass  (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:)  )
+      snwice_mass  (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:)  )
       snwice_mass_b(:,:) = snwice_mass(:,:)
       snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcisf.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcisf.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcisf.F90	(revision 9937)
@@ -52,9 +52,8 @@
    LOGICAL, PUBLIC ::   l_isfcpl = .false.       !: isf recieved from oasis
 
-   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     !: specific heat of ice shelf             [J/kg/K]
+   REAL(wp)        , SAVE ::   rcp_isf  = 2000.0_wp     !: specific heat of ice shelf             [J/kg/K]
    REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    !: heat diffusivity through the ice-shelf [m2/s]
    REAL(wp), PUBLIC, SAVE ::   rho_isf  = 920.0_wp      !: volumic mass of ice shelf              [kg/m3]
    REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      !: air temperature on top of ice shelf    [C]
-   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    !: latent heat of fusion of ice shelf     [J/kg]
 
 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
@@ -114,5 +113,5 @@
             ! compute fwf and heat flux
             IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt)
-            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rlfusisf  ! heat        flux
+            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rLfus    ! heat flux
             ENDIF
             !
@@ -127,5 +126,5 @@
                fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
             ENDIF
-            qisf(:,:)   = fwfisf(:,:) * rlfusisf             ! heat flux
+            qisf(:,:)   = fwfisf(:,:) * rLfus                   ! heat flux
             stbl(:,:)   = soce
             !
@@ -137,5 +136,5 @@
                fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1)            ! fwf
             ENDIF
-            qisf(:,:)   = fwfisf(:,:) * rlfusisf               ! heat flux
+            qisf(:,:)   = fwfisf(:,:) * rLfus                     ! heat flux
             stbl(:,:)   = soce
             !
@@ -308,7 +307,7 @@
       qisf    (:,:)    = 0._wp   ;   fwfisf  (:,:) = 0._wp
       risf_tsc(:,:,:)  = 0._wp   ;   fwfisf_b(:,:) = 0._wp
-      !
-      ! define isf tbl tickness, top and bottom indice
-      SELECT CASE ( nn_isf )
+
+      SELECT CASE ( nn_isf )      ! define isf tbl tickness, top and bottom indice
+      !
       CASE ( 1 ) 
          IF(lwp) WRITE(numout,*)
@@ -455,5 +454,5 @@
                            & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk)
              
-               fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                  
+               fwfisf(ji,jj) = qisf(ji,jj) * r1_Lfus                        ! fresh water flux kg/(m2s)                  
                fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
                !add to salinity trend
@@ -527,5 +526,5 @@
                DO ji = 1, jpi
                   zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rho0*(ttbl(ji,jj)-zfrz(ji,jj))
-                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf
+                  zfwflx(ji,jj) = - zhtflx(ji,jj) * r1_Lfus
                END DO
             END DO
@@ -545,6 +544,6 @@
                   ! compute coeficient to solve the 2nd order equation
                   zeps1 = rcp*rho0*zgammat(ji,jj)
-                  zeps2 = rlfusisf*rho0*zgammas(ji,jj)
-                  zeps3 = rho_isf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps)
+                  zeps2 = rLfus*rho0*zgammas(ji,jj)
+                  zeps3 = rho_isf*rcp_isf*rkappa/MAX(risfdep(ji,jj),zeps)
                   zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
                   zeps6 = zeps4-ttbl(ji,jj)
Index: /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcrnf.F90
===================================================================
--- /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcrnf.F90	(revision 9936)
+++ /NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcrnf.F90	(revision 9937)
@@ -128,5 +128,5 @@
             END WHERE
             WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg
-               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rho0 - rnf(:,:) * rlfusisf * r1_rho0_rcp
+               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rho0 - rnf(:,:) * rLfus * r1_rho0_rcp
             END WHERE
          ELSE                                                        ! use SST as runoffs temperature
