Index: LMDZ6/branches/Ocean_skin/libf/phylmd/calcul_fluxs_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/calcul_fluxs_mod.F90	(revision 3676)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/calcul_fluxs_mod.F90	(revision 3687)
@@ -13,9 +13,10 @@
        fqsat, petAcoef, peqAcoef, petBcoef, peqBcoef, &
        tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
-       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
+       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
  
     
     USE dimphy, ONLY : klon
     USE indice_sol_mod
+    use sens_heat_rain_m, only: sens_heat_rain
 
     INCLUDE "clesphys.h"
@@ -71,9 +72,12 @@
     REAL, DIMENSION(klon), INTENT(IN)    :: ps, q1lay
     REAL, DIMENSION(klon), INTENT(IN)    :: tsurf, p1lay, cal, beta, cdragh,cdragq
-    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow ! pas utiles
+    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
     REAL, DIMENSION(klon), INTENT(IN)    :: radsol, dif_grnd
     REAL, DIMENSION(klon), INTENT(IN)    :: t1lay, u1lay, v1lay,gustiness
     REAL,                  INTENT(IN)    :: fqsat ! correction factor on qsat (generally 0.98 over salty water, 1 everywhere else)
 
+    real, intent(in), optional:: rhoa(:) ! (knon)
+    ! density of moist air  (kg / m3)
+
 ! Parametres entree-sorties
 !****************************************************************************************
@@ -85,5 +89,5 @@
     REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new, evap, fluxsens, fluxlat
     REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l
-    REAL, DIMENSION(klon), OPTIONAL      :: sens_prec_liq, sens_prec_sol
+    REAL, intent(out), OPTIONAL:: sens_prec_liq(:), sens_prec_sol(:) ! (knon)
     REAL, DIMENSION(klon), OPTIONAL      :: lat_prec_liq, lat_prec_sol
 
@@ -140,6 +144,4 @@
     dflux_s = 0.
     dflux_l = 0.
-    if (PRESENT(sens_prec_liq)) sens_prec_liq = 0.
-    if (PRESENT(sens_prec_sol)) sens_prec_sol = 0.
     if (PRESENT(lat_prec_liq)) lat_prec_liq = 0.
     if (PRESENT(lat_prec_sol)) lat_prec_sol = 0.
@@ -267,10 +269,8 @@
 !                          precip_snow(i)
 ! On calcule par rapport a T=0
-       if (PRESENT(sens_prec_liq))                   &
-         sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * &
-                          precip_rain(i)
-       if (PRESENT(sens_prec_sol))                   &
-         sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * &
-                          precip_snow(i)
+       if (PRESENT(sens_prec_liq)) sens_prec_liq(i) &
+            = - sens_heat_rain(precip_rain(i) + precip_snow(i), t1lay(i), &
+            q1lay(i), rhoa(i), rlvtt, tsurf_new(i), ps(i))
+       if (PRESENT(sens_prec_sol)) sens_prec_sol(i) = 0.
        if (PRESENT(lat_prec_liq))                    &
          lat_prec_liq(i) =  precip_rain(i) * (RLVTT - RLVTT) 
@@ -280,13 +280,5 @@
 
     
-!       if (PRESENT(sens_prec_liq))                  &
-!          WRITE(*,*)' calculs_fluxs sens_prec_liq (min, max)', &
-!          MINVAL(sens_prec_liq(1:knon)), MAXVAL(sens_prec_liq(1:knon))
-!       if (PRESENT(sens_prec_sol))                  &
-!          WRITE(*,*)' calculs_fluxs sens_prec_sol (min, max)', &
-!          MINVAL(sens_prec_sol(1:knon)), MAXVAL(sens_prec_sol(1:knon))
- 
-!
-!****************************************************************************************
+!**************************************************************************
 !
   END SUBROUTINE calcul_fluxs
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/cpl_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/cpl_mod.F90	(revision 3676)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/cpl_mod.F90	(revision 3687)
@@ -669,5 +669,5 @@
     REAL, DIMENSION(klon), INTENT(IN)       :: evap, tsurf, fder, albsol
     REAL, DIMENSION(klon), INTENT(IN)       :: taux, tauy, windsp
-    REAL, DIMENSION(klon), INTENT(IN)       :: sens_prec_liq, sens_prec_sol
+    REAL, INTENT(IN):: sens_prec_liq(:), sens_prec_sol(:) ! (knon)
     REAL, DIMENSION(klon), INTENT(IN)       :: lat_prec_liq, lat_prec_sol
     real, intent(in):: t_int(:) ! (klon) ocean-air interface temperature, in K
@@ -905,5 +905,5 @@
     REAL, DIMENSION(klon), INTENT(IN)       :: albsol, taux, tauy
     REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
-    REAL, DIMENSION(klon), INTENT(IN)       :: sens_prec_liq, sens_prec_sol
+    REAL, INTENT(IN):: sens_prec_liq(:), sens_prec_sol(:) ! (knon)
     REAL, DIMENSION(klon), INTENT(IN)       :: lat_prec_liq, lat_prec_sol
     LOGICAL, INTENT(IN)                     :: lafin
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_cpl_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_cpl_mod.F90	(revision 3676)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_cpl_mod.F90	(revision 3687)
@@ -55,5 +55,5 @@
        radsol, snow, agesno, &
        qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l, sss, s_int)
+       tsurf_new, dflux_s, dflux_l, sens_prec_liq, sss, s_int, rhoa)
 
 !
@@ -93,4 +93,5 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
     real, intent(in):: s_int(:) ! (knon) ocean-air interface salinity, in ppt
+    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
 
 ! In/Output arguments
@@ -107,4 +108,5 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
+    REAL, intent(out):: sens_prec_liq(:) ! (knon)
 
     REAL, INTENT(OUT):: sss(:) ! (klon)
@@ -122,5 +124,5 @@
     REAL, DIMENSION(klon) :: u1_lay, v1_lay
     LOGICAL               :: check=.FALSE.
-    REAL, DIMENSION(klon) :: sens_prec_liq, sens_prec_sol    
+    REAL sens_prec_sol(knon)  
     REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol    
 
@@ -145,5 +147,5 @@
     dif_grnd = 0.
     agesno(:) = 0.
-    sens_prec_liq = 0.; sens_prec_sol = 0.; lat_prec_liq = 0.; lat_prec_sol = 0.
+    lat_prec_liq = 0.; lat_prec_sol = 0.
     
 
@@ -160,5 +162,5 @@
          f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
-         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
+         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
 
     if (activate_ocean_skin == 2) then
@@ -226,5 +228,5 @@
        radsol, snow, qsurf, &
        alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l)
+       tsurf_new, dflux_s, dflux_l, rhoa)
 !
 ! This subroutine treats the ocean where there is ice. The subroutine first receives 
@@ -261,4 +263,5 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
     REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
+    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
 
 ! In/output arguments
@@ -288,5 +291,5 @@
     REAL, DIMENSION(klon)   :: u0, v0
     REAL, DIMENSION(klon)   :: u1_lay, v1_lay
-    REAL, DIMENSION(klon) :: sens_prec_liq, sens_prec_sol    
+    REAL sens_prec_liq(knon), sens_prec_sol(knon)   
     REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol    
 
@@ -296,5 +299,5 @@
     IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon 
 
-    sens_prec_liq = 0.; sens_prec_sol = 0.; lat_prec_liq = 0.; lat_prec_sol = 0.
+    lat_prec_liq = 0.; lat_prec_sol = 0.
 
 !****************************************************************************************
@@ -329,5 +332,5 @@
          f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
-         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
+         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
     do j = 1, knon
       i = knindex(j)
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_forced_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_forced_mod.F90	(revision 3676)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_forced_mod.F90	(revision 3687)
@@ -22,5 +22,5 @@
        radsol, snow, agesno, & 
        qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l)
+       tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa)
 !
 ! This subroutine treats the "open ocean", all grid points that are not entierly covered
@@ -55,4 +55,5 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
     REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
+    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
 
 ! In/Output arguments
@@ -69,4 +70,5 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
+    REAL, intent(out):: sens_prec_liq(:) ! (knon)
 
 ! Local variables
@@ -78,5 +80,5 @@
     REAL, DIMENSION(klon)       :: u1_lay, v1_lay
     LOGICAL                     :: check=.FALSE.
-    REAL, DIMENSION(klon) :: sens_prec_liq, sens_prec_sol    
+    REAL sens_prec_sol(knon) 
     REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol    
 
@@ -111,5 +113,5 @@
     alb_neig(:) = 0.
     agesno(:) = 0.
-    sens_prec_liq = 0.; sens_prec_sol = 0.; lat_prec_liq = 0.; lat_prec_sol = 0.
+    lat_prec_liq = 0.; lat_prec_sol = 0.
 
 ! Suppose zero surface speed
@@ -127,5 +129,5 @@
          f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
-         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
+         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
     if (activate_ocean_skin == 2) tsurf_new = tsurf_lim
 
@@ -158,5 +160,5 @@
        radsol, snow, qsol, agesno, tsoil, &
        qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l)
+       tsurf_new, dflux_s, dflux_l, rhoa)
 !
 ! This subroutine treats the ocean where there is ice. 
@@ -191,4 +193,5 @@
     REAL, DIMENSION(klon), INTENT(IN)    :: ps
     REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
+    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
 
 ! In/Output arguments
@@ -220,5 +223,5 @@
     REAL, DIMENSION(klon)       :: u0, v0
     REAL, DIMENSION(klon)       :: u1_lay, v1_lay
-    REAL, DIMENSION(klon)       :: sens_prec_liq, sens_prec_sol    
+    REAL sens_prec_liq(knon), sens_prec_sol (knon)   
     REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol    
 
@@ -254,5 +257,5 @@
 
     beta = 1.0
-    sens_prec_liq = 0.; sens_prec_sol = 0.; lat_prec_liq = 0.; lat_prec_sol = 0.
+    lat_prec_liq = 0.; lat_prec_sol = 0.
 
 ! Suppose zero surface speed
@@ -267,5 +270,5 @@
          f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
-         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
+         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
     do j = 1, knon
       i = knindex(j)
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/pbl_surface_mod.F90	(revision 3676)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/pbl_surface_mod.F90	(revision 3687)
@@ -293,5 +293,5 @@
     USE ioipsl_getin_p_mod, ONLY : getin_p
     use phys_state_var_mod, only: s_int, ds_ns, dt_ns
-    use phys_output_var_mod, only: t_int, dter, dser, tkt, tks, rf, taur
+    use phys_output_var_mod, only: t_int, dter, dser, tkt, tks, taur
 #ifdef CPP_XIOS
     USE wxios, ONLY: missing_val
@@ -831,6 +831,6 @@
 
     real, DIMENSION(klon):: yt_int, ys_int, yds_ns, ydt_ns, ydter, ydser, &
-         ytkt, ytks, yrf, ytaur
-    ! compression of t_int, s_int, ds_ns, dt_ns, dter, dser, tkt, tks, rf,
+         ytkt, ytks, ytaur
+    ! compression of t_int, s_int, ds_ns, dt_ns, dter, dser, tkt, tks,
     ! taur on ocean points
 
@@ -1956,5 +1956,5 @@
                y_flux_u1, y_flux_v1, yt_int(:knon), ys_int(:knon), &
                yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
-               ytkt(:knon), ytks(:knon), yrf(:knon), ytaur(:knon))
+               ytkt(:knon), ytks(:knon), ytaur(:knon))
       IF (prt_level >=10) THEN
           print *,'arg de surf_ocean: ycdragh ',ycdragh
@@ -2569,5 +2569,4 @@
           tkt = missing_val
           tks = missing_val
-          rf = missing_val
           taur = missing_val
           
@@ -2580,5 +2579,4 @@
           tkt(ni(:knon)) = ytkt(:knon)
           tks(ni(:knon)) = ytks(:knon)
-          rf(ni(:knon)) = yrf(:knon)
           taur(ni(:knon)) = ytaur(:knon)
        end if
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3676)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3687)
@@ -1979,9 +1979,4 @@
         "thickness of salinity microlayer", "m", '')
 
-   type(ctrl_out), save:: o_rf &
-        = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'rf', &
-        "sensible heat flux at the surface due to rainfall, positive upward", &
-        "W m-2", '')
-
    type(ctrl_out), save:: o_taur &
         = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'taur', &
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_var_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_var_mod.F90	(revision 3676)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_var_mod.F90	(revision 3687)
@@ -152,12 +152,8 @@
   ! épaisseur (m) de la couche de diffusion de masse (microlayer)
   
-  REAL, SAVE, ALLOCATABLE:: rf(:)
-  ! sensible heat flux at the surface due to rainfall, in  W m-2
-  ! positive upward
-  
   REAL, SAVE, ALLOCATABLE:: taur(:)
   ! momentum flux due to rain, in Pa
   
-  !$OMP THREADPRIVATE(t_int, dter, dser, tkt, tks, rf, taur)
+  !$OMP THREADPRIVATE(t_int, dter, dser, tkt, tks, taur)
 
 CONTAINS
@@ -221,5 +217,5 @@
 
     if (activate_ocean_skin >= 1) allocate(t_int(klon), dter(klon), &
-         dser(klon), tkt(klon), tks(klon), rf(klon), taur(klon))
+         dser(klon), tkt(klon), tks(klon), taur(klon))
 
   END SUBROUTINE phys_output_var_init
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_write_mod.F90	(revision 3676)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_write_mod.F90	(revision 3687)
@@ -203,5 +203,5 @@
          o_flx_co2_ocean, o_flx_co2_land, o_flx_co2_ff, o_flx_co2_bb, &
          o_t_int, o_s_int, o_ds_ns, o_dt_ns, o_dter, o_dser, o_tkt, o_tks, &
-         o_rf, o_taur
+         o_taur
 
 #ifdef CPP_StratAer
@@ -360,5 +360,5 @@
          alt_tropo, &
 !Ionela
-         ok_4xCO2atm, t_int, dter, dser, tkt, tks, rf, taur
+         ok_4xCO2atm, t_int, dter, dser, tkt, tks, taur
 
     USE ocean_slab_mod, ONLY: nslay, tslab, slab_bilg, tice, seaice, &
@@ -2315,5 +2315,4 @@
           CALL histwrite_phy(o_tkt, tkt)
           CALL histwrite_phy(o_tks, tks)
-          CALL histwrite_phy(o_rf, rf)
           CALL histwrite_phy(o_taur, taur)
        end if
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/surf_ocean_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/surf_ocean_mod.F90	(revision 3676)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/surf_ocean_mod.F90	(revision 3687)
@@ -20,6 +20,5 @@
        z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, & 
        tsurf_new, dflux_s, dflux_l, lmt_bils, &
-       flux_u1, flux_v1, t_int, s_int, ds_ns, dt_ns, dter, dser, tkt, tks, rf, &
-       taur)
+       flux_u1, flux_v1, t_int, s_int, ds_ns, dt_ns, dter, dser, tkt, tks, taur)
 
     use albedo, only: alboc, alboc_cd
@@ -32,5 +31,4 @@
     USE indice_sol_mod, ONLY : nbsrf, is_oce
     USE limit_read_mod
-    use sens_heat_rain_m, only: sens_heat_rain
     use config_ocean_skin_m, only: activate_ocean_skin
     !
@@ -122,8 +120,4 @@
     ! épaisseur (m) de la couche de diffusion de masse (microlayer)
 
-    REAL, intent(out):: rf(:) ! (knon)
-    ! sensible heat flux at the surface due to rainfall, in  W m-2
-    ! positive upward
-
     REAL, intent(out):: taur(:) ! (knon)
     ! momentum flux due to rain, in Pa
@@ -139,6 +133,5 @@
     CHARACTER(len=20),PARAMETER :: modname="surf_ocean"
     real rhoa(knon) ! density of moist air  (kg / m3)
-    real xlv(knon) ! chaleur latente d'évaporation (J / kg)
-    real precip_tot(knon) ! rain + snow
+    REAL sens_prec_liq(knon)
 
     real sss(klon)
@@ -175,5 +168,5 @@
     ENDIF
 
-
+    rhoa = PS(:KNON) / (Rd * temp_air(:knon) * (1. + retv * spechum(:knon)))
     !******************************************************************************
     ! Switch according to type of ocean (couple, slab or forced)
@@ -191,5 +184,5 @@
             radsol, snow, agesno, &
             qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-            tsurf_new, dflux_s, dflux_l, sss, s_int)
+            tsurf_new, dflux_s, dflux_l, sens_prec_liq, sss, s_int, rhoa)
 
     CASE('slab')
@@ -214,5 +207,5 @@
             radsol, snow, agesno, &
             qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-            tsurf_new, dflux_s, dflux_l)
+            tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa)
     END SELECT
 
@@ -319,16 +312,13 @@
 
     if (activate_ocean_skin >= 1) then
-       rhoa = PS(:KNON) / (Rd * temp_air(:knon) * (1. + retv * spechum(:knon)))
-       xlv = rlvtt
-       precip_tot = precip_rain(:knon) + precip_snow(:knon)
-       rf =  sens_heat_rain(precip_tot, temp_air(:knon), spechum(:knon), rhoa, &
-            xlv, tsurf_in(:knon), ps(:knon))
        if (type_ocean /= 'couple') sss(:knon) = 35.
        call bulk_flux(tkt, tks, taur, dter, dser, t_int, s_int, ds_ns, dt_ns, &
             u = windsp(:knon), t_ocean_1 = tsurf_new(:knon), s1 = sss(:knon), &
-            rain = precip_tot, hf = - fluxsens(:knon), hlb = - fluxlat(:knon), &
+            rain = precip_rain(:knon) + precip_snow(:knon), &
+            hf = - fluxsens(:knon), hlb = - fluxlat(:knon), &
             rnl = - lwnet(:knon), &
             tau = sqrt(flux_u1(:knon)**2 + flux_v1(:knon)**2), rhoa = rhoa, &
-            xlv = xlv, rf = rf, dtime = dtime, rns = swnet(:knon))
+            xlv = [(rlvtt, i = 1, knon)], rf = - sens_prec_liq, dtime = dtime, &
+            rns = swnet(:knon))
        if (activate_ocean_skin == 2) tsurf_new(:knon) = t_int
     end if
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/surf_seaice_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/surf_seaice_mod.F90	(revision 3676)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/surf_seaice_mod.F90	(revision 3687)
@@ -37,4 +37,7 @@
     INCLUDE "dimsoil.h"
     INCLUDE "clesphys.h"
+
+    INCLUDE "YOMCST.h"
+    ! for rd and retv
 
 ! Input arguments
@@ -87,5 +90,7 @@
     REAL, DIMENSION(klon) :: alb1_new,alb2_new
 !albedo SB <<<
-!
+
+    real rhoa(knon) ! density of moist air  (kg / m3)
+
 ! End definitions
 !****************************************************************************************
@@ -98,4 +103,6 @@
     radsol(:) = 0.0
     radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
+
+    rhoa = PS(:KNON) / (Rd * temp_air(:knon) * (1. + retv * spechum(:knon)))
 
 !****************************************************************************************
@@ -116,5 +123,5 @@
             radsol, snow, qsurf, &
             alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-            tsurf_new, dflux_s, dflux_l)
+            tsurf_new, dflux_s, dflux_l, rhoa)
        
     ELSE IF (type_ocean == 'slab'.AND.version_ocean=='sicINT') THEN
@@ -138,5 +145,5 @@
             radsol, snow, qsol, agesno, tsoil, &
             qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-            tsurf_new, dflux_s, dflux_l)
+            tsurf_new, dflux_s, dflux_l, rhoa)
 
     END IF
