Index: LMDZ6/branches/Ocean_skin/libf/phylmd/cpl_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/cpl_mod.F90	(revision 3626)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/cpl_mod.F90	(revision 3627)
@@ -48,4 +48,8 @@
   REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_snow, cpl_evap, cpl_tsol
   !$OMP THREADPRIVATE(cpl_snow,cpl_evap,cpl_tsol)
+
+  REAL, ALLOCATABLE, SAVE:: cpl_t_int(:)
+  !$OMP THREADPRIVATE(cpl_t_int)
+  
   REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_fder, cpl_albe, cpl_taux, cpl_tauy
   !$OMP THREADPRIVATE(cpl_fder,cpl_albe,cpl_taux,cpl_tauy)
@@ -66,4 +70,9 @@
   REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_sit     ! sea ice temperature
   !$OMP THREADPRIVATE(read_sit)
+
+  REAL, ALLOCATABLE, SAVE:: read_sss(:, :)
+  ! bulk salinity of the surface layer of the ocean, in ppt
+  !$OMP THREADPRIVATE(read_sss)
+
   REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: read_sic     ! sea ice fraction
   !$OMP THREADPRIVATE(read_sic)
@@ -84,4 +93,8 @@
   REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_snow2D, cpl_evap2D, cpl_tsol2D
   !$OMP THREADPRIVATE(cpl_snow2D, cpl_evap2D, cpl_tsol2D)
+
+  REAL, ALLOCATABLE, SAVE:: cpl_t_int_2D(:,:)
+  !$OMP THREADPRIVATE(cpl_t_int_2D)
+
   REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_fder2D, cpl_albe2D
   !$OMP THREADPRIVATE(cpl_fder2D, cpl_albe2D)
@@ -122,4 +135,5 @@
     USE geometry_mod, ONLY : longitude_deg, latitude_deg, ind_cell_glo, cell_area
     USE ioipsl_getin_p_mod, ONLY: getin_p
+    use config_ocean_skin_m, only: activate_ocean_skin
 
 ! Input arguments
@@ -217,4 +231,15 @@
     ALLOCATE(read_sit(nbp_lon, jj_nb), stat = error)
     sum_error = sum_error + error
+
+    if (activate_ocean_skin >= 1) then
+       ALLOCATE(read_sss(nbp_lon, jj_nb), stat = error)
+       sum_error = sum_error + error
+    
+       if (activate_ocean_skin == 2) then
+          ALLOCATE(cpl_t_int(klon), stat = error)
+          sum_error = sum_error + error
+       end if
+    end if
+
     ALLOCATE(read_alb_sic(nbp_lon, jj_nb), stat = error)
     sum_error = sum_error + error
@@ -376,4 +401,5 @@
     USE time_phylmdz_mod, ONLY: start_time, itau_phy
     USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
+    use config_ocean_skin_m, only: activate_ocean_skin
 
     INCLUDE "YOMCST.h"
@@ -437,4 +463,5 @@
        read_alb_sic(:,:) = tab_read_flds(:,:,idr_icealw)  ! Albedo at sea ice
        read_sit(:,:)     = tab_read_flds(:,:,idr_icetem)  ! Sea ice temperature
+       if (activate_ocean_skin >= 1) read_sss(:,:) = tab_read_flds(:,:,idr_sss)
 !$OMP END MASTER
 
@@ -494,5 +521,6 @@
 !
 
-  SUBROUTINE cpl_receive_ocean_fields(knon, knindex, tsurf_new, u0_new, v0_new)
+  SUBROUTINE cpl_receive_ocean_fields(knon, knindex, tsurf_new, u0_new, &
+       v0_new, sss)
 !
 ! This routine returns the field for the ocean that has been read from the coupler
@@ -502,4 +530,5 @@
     USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, fco2_ocn_day
     USE indice_sol_mod
+    use config_ocean_skin_m, only: activate_ocean_skin
 
 ! Input arguments
@@ -511,4 +540,8 @@
 !*************************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT)      :: tsurf_new
+
+    REAL, INTENT(OUT):: sss(:) ! (klon)
+    ! bulk salinity of the surface layer of the ocean, in ppt
+
     REAL, DIMENSION(klon), INTENT(OUT)      :: u0_new
     REAL, DIMENSION(klon), INTENT(OUT)      :: v0_new
@@ -525,4 +558,5 @@
 !*************************************************************************************
     CALL cpl2gath(read_sst, tsurf_new, knon, knindex)
+    if (activate_ocean_skin >= 1) CALL cpl2gath(read_sss, sss, knon, knindex)
     CALL cpl2gath(read_sic, sic_new, knon, knindex)
     CALL cpl2gath(read_u0, u0_new, knon, knindex)
@@ -611,14 +645,17 @@
        swdown, lwdown, fluxlat, fluxsens, &
        precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy, windsp,&
-       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
-!
-! This subroutine cumulates some fields for each time-step during a coupling 
-! period. At last time-step in a coupling period the fields are transformed to the 
-! grid accepted by the coupler. No sending to the coupler will be done from here 
-! (it is done in cpl_send_seaice_fields).
-!
+       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, t_int)
+
+    ! This subroutine cumulates some fields for each time-step during
+    ! a coupling period. At last time-step in a coupling period the
+    ! fields are transformed to the grid accepted by the coupler. No
+    ! sending to the coupler will be done from here (it is done in
+    ! cpl_send_seaice_fields). Crucial hypothesis is that the surface
+    ! fractions do not change between coupling time-steps.
+
     USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
     USE indice_sol_mod
     USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
+    use config_ocean_skin_m, only: activate_ocean_skin
 
 ! Input arguments
@@ -634,4 +671,5 @@
     REAL, DIMENSION(klon), INTENT(IN)       :: sens_prec_liq, sens_prec_sol
     REAL, DIMENSION(klon), INTENT(IN)       :: lat_prec_liq, lat_prec_sol
+    real, intent(in):: t_int(:) ! (klon) ocean-air interface temperature, in K
 
 ! Local variables
@@ -669,4 +707,5 @@
        cpl_taumod(1:knon,cpl_index) = 0.0
        IF (carbon_cycle_cpl) cpl_atm_co2(1:knon,cpl_index) = 0.0
+       if (activate_ocean_skin == 2) cpl_t_int(:knon) = 0.
     ENDIF
        
@@ -710,4 +749,7 @@
 !!---OB: this is correct but why knindex ??
        ENDIF
+
+       if (activate_ocean_skin == 2) &
+            cpl_t_int(ig) = cpl_t_int(ig) + t_int(ig) / REAL(nexca)
      ENDDO
 
@@ -755,4 +797,9 @@
              sum_error = sum_error + error
           ENDIF
+
+          if (activate_ocean_skin == 2) then
+             ALLOCATE(cpl_t_int_2D(nbp_lon, jj_nb), stat = error)
+             sum_error = sum_error + error
+          end if
 
           IF (sum_error /= 0) THEN
@@ -810,5 +857,7 @@
        IF (carbon_cycle_cpl) &
             CALL gath2cpl(cpl_atm_co2(:,cpl_index), cpl_atm_co22D(:,:), knon, knindex)
-   ENDIF
+       if (activate_ocean_skin == 2) &
+            CALL gath2cpl(cpl_t_int, cpl_t_int_2D, knon, knindex)
+    ENDIF
 
   END SUBROUTINE cpl_send_ocean_fields
@@ -1145,4 +1194,5 @@
     USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
     USE time_phylmdz_mod, ONLY: start_time, itau_phy
+    use config_ocean_skin_m, only: activate_ocean_skin
 ! Some includes
 !    
@@ -1202,4 +1252,5 @@
     tab_flds(:,:,ids_qraiic) = cpl_sens_rain2D(:,:,2)
     tab_flds(:,:,ids_qsnoic) = cpl_sens_snow2D(:,:,2)
+    if (activate_ocean_skin == 2) tab_flds(:, :, ids_t_int) = cpl_t_int_2D
     
     IF (version_ocean=='nemo') THEN
@@ -1438,4 +1489,6 @@
     ENDIF
 
+    if (activate_ocean_skin == 2) deallocate(cpl_t_int_2d)
+
     IF (sum_error /= 0) THEN
        abort_message='Pb in deallocation of cpl_xxxx2D coupling variables'
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/oasis.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/oasis.F90	(revision 3626)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/oasis.F90	(revision 3627)
@@ -59,7 +59,10 @@
   INTEGER, PARAMETER :: ids_qraiic = 28
   INTEGER, PARAMETER :: ids_qsnoic = 29
-  INTEGER, PARAMETER :: maxsend    = 29  ! Maximum number of fields to send
+  INTEGER, PARAMETER :: ids_t_int = 30
+  
+  INTEGER, PARAMETER :: maxsend    = 30  ! Maximum number of fields to send
   
   ! Id for fields received from ocean
+
   INTEGER, PARAMETER :: idr_sisutw = 1
   INTEGER, PARAMETER :: idr_icecov = 2
@@ -70,5 +73,9 @@
   INTEGER, PARAMETER :: idr_curenz = 7
   INTEGER, PARAMETER :: idr_oceco2 = 8
-  INTEGER, PARAMETER :: maxrecv    = 8  ! Maximum number of fields to receive
+
+  INTEGER, PARAMETER :: idr_sss = 9
+  ! bulk salinity of the surface layer of the ocean, in ppt
+
+  INTEGER, PARAMETER :: maxrecv    = 9  ! Maximum number of fields to receive
   
 
@@ -183,4 +190,9 @@
     infosend(ids_calvin)%action = .TRUE. ; infosend(ids_calvin)%name = 'COCALVIN'
     
+    if (activate_ocean_skin == 2) then
+       infosend(ids_t_int)%action = .TRUE.
+       infosend(ids_t_int)%name = 'T_int'
+    end if
+           
     IF (version_ocean=='nemo') THEN
         infosend(ids_shftot)%action = .TRUE. ; infosend(ids_shftot)%name = 'COQSRMIX'
@@ -219,4 +231,9 @@
    inforecv(idr_icealw)%action = .TRUE. ; inforecv(idr_icealw)%name = 'SIICEALW'
    inforecv(idr_icetem)%action = .TRUE. ; inforecv(idr_icetem)%name = 'SIICTEMW'
+
+   if (activate_ocean_skin >= 1) then
+      inforecv(idr_sss)%action = .TRUE.
+      inforecv(idr_sss)%name = 'salinity'
+   end if
    
    IF (cpl_current ) THEN
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_cpl_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_cpl_mod.F90	(revision 3626)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_cpl_mod.F90	(revision 3627)
@@ -55,5 +55,5 @@
        radsol, snow, agesno, &
        qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l)
+       tsurf_new, dflux_s, dflux_l, sss)
 
 !
@@ -63,9 +63,9 @@
 !
     USE dimphy,           ONLY : klon
-    USE cpl_mod
     USE calcul_fluxs_mod
     USE indice_sol_mod
     USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
-    USE cpl_mod,             ONLY : gath2cpl
+    USE cpl_mod, ONLY : gath2cpl, cpl_receive_ocean_fields, &
+         cpl_send_ocean_fields
     use config_ocean_skin_m, only: activate_ocean_skin
 
@@ -106,4 +106,7 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
+
+    REAL, INTENT(OUT):: sss(:) ! (klon)
+    ! bulk salinity of the surface layer of the ocean, in ppt
   
 
@@ -130,5 +133,6 @@
 !
 !****************************************************************************************
-    CALL cpl_receive_ocean_fields(knon, knindex, tsurf_cpl, u0_cpl, v0_cpl)
+    CALL cpl_receive_ocean_fields(knon, knindex, tsurf_cpl, u0_cpl, v0_cpl, &
+         sss)
 
 !****************************************************************************************
@@ -156,5 +160,12 @@
          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
          sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
-    if (activate_ocean_skin == 2) tsurf_new = tsurf_cpl
+
+    if (activate_ocean_skin == 2) then
+       ! tsurf_new was set to tsurf_in in calcul_flux, correct it to
+       ! the new bulk SST tsurf_cpl:
+       tsurf_new = tsurf_cpl
+    end if
+
+    ! assertion: tsurf_new == tsurf_cpl
     
     do j = 1, knon
@@ -197,5 +208,5 @@
          swnet, lwnet, fluxlat, fluxsens, &
          precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, flux_u1, flux_v1, windsp,&
-         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
+         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, tsurf_in)
     
 
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/surf_ocean_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/surf_ocean_mod.F90	(revision 3626)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/surf_ocean_mod.F90	(revision 3627)
@@ -57,5 +57,5 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: rmu0  
     REAL, DIMENSION(klon), INTENT(IN)        :: fder
-    REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
+    REAL, INTENT(IN):: tsurf_in(klon) ! defined only for subscripts 1:knon
     REAL, DIMENSION(klon), INTENT(IN)        :: p1lay,z1lay ! pression (Pa) et altitude (m) du premier niveau
     REAL, DIMENSION(klon), INTENT(IN)        :: cdragh
@@ -102,6 +102,6 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
 
-    REAL, intent(out):: t_int(:) ! (knon) interface temperature, in K
-    real, intent(out):: s_int(:) ! (knon) interface salinity, in ppt
+    REAL, intent(out):: t_int(:) ! (knon) ocean-air interface temperature, in K
+    real, intent(out):: s_int(:) ! (knon) ocean-air interface salinity, in ppt
 
     REAL, intent(out):: dter(:) ! (knon)
@@ -111,6 +111,6 @@
 
     REAL, intent(out):: dser(:) ! (knon)
-    ! Temperature variation in the diffusive microlayer, that is
-    ! subskin temperature minus ocean-air interface temperature. In K.
+    ! Salinity variation in the diffusive microlayer, that is
+    ! ocean-air interface salinity minus subskin salinity. In ppt.
 
     REAL, intent(out):: tkt(:) ! (knon)
@@ -140,5 +140,9 @@
     real xlv(knon) ! chaleur latente d'évaporation (J / kg)
     real precip_tot(knon) ! rain + snow
-    real S1(knon) ! salinity at depth_1, in ppt
+
+    real sss(klon)
+    ! Bulk salinity of the surface layer of the ocean, in ppt. (Only
+    ! defined for subscripts 1:knon, but we have to declare it with
+    ! size klon because of the coupling machinery.)
 
     ! End definition
@@ -186,5 +190,5 @@
             radsol, snow, agesno, &
             qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-            tsurf_new, dflux_s, dflux_l)
+            tsurf_new, dflux_s, dflux_l, sss)
 
     CASE('slab')
@@ -319,7 +323,7 @@
        rf =  sens_heat_rain(precip_tot, temp_air(:knon), spechum(:knon), rhoa, &
             xlv, tsurf_in(:knon), ps(:knon))
-       s1 = 35.
+       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 = s1, &
+            u = windsp(:knon), t_ocean_1 = tsurf_new(:knon), s1 = sss(:knon), &
             rain = precip_tot, hf = - fluxsens(:knon), hlb = - fluxlat(:knon), &
             rnl = - lwnet(:knon), &
