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 9923)
+++ 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 9923)
+++ 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 9923)
+++ 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 9923)
+++ 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 9923)
+++ 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 9923)
+++ 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 9923)
+++ 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
