Index: LMDZ6/branches/Ocean_skin/libf/phylmd/cpl_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/cpl_mod.F90	(revision 3739)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/cpl_mod.F90	(revision 3740)
@@ -49,6 +49,6 @@
   !$OMP THREADPRIVATE(cpl_snow,cpl_evap,cpl_tsol)
 
-  REAL, ALLOCATABLE, SAVE:: cpl_t_int(:), cpl_s_int(:)
-  !$OMP THREADPRIVATE(cpl_t_int, cpl_s_int)
+  REAL, ALLOCATABLE, SAVE:: cpl_delta_temp(:), cpl_s_int(:)
+  !$OMP THREADPRIVATE(cpl_delta_temp, cpl_s_int)
   
   REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_fder, cpl_albe, cpl_taux, cpl_tauy
@@ -94,6 +94,6 @@
   !$OMP THREADPRIVATE(cpl_snow2D, cpl_evap2D, cpl_tsol2D)
 
-  REAL, ALLOCATABLE, SAVE:: cpl_t_int_2D(:,:), cpl_s_int_2D(:,:)
-  !$OMP THREADPRIVATE(cpl_t_int_2D, cpl_s_int_2D)
+  REAL, ALLOCATABLE, SAVE:: cpl_delta_temp_2D(:,:), cpl_s_int_2D(:,:)
+  !$OMP THREADPRIVATE(cpl_delta_temp_2D, cpl_s_int_2D)
 
   REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cpl_fder2D, cpl_albe2D
@@ -237,5 +237,5 @@
     
        if (activate_ocean_skin == 2) then
-          ALLOCATE(cpl_t_int(klon), cpl_s_int(klon), stat = error)
+          ALLOCATE(cpl_delta_temp(klon), cpl_s_int(klon), stat = error)
           sum_error = sum_error + error
        end if
@@ -645,5 +645,6 @@
        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, t_int, s_int)
+       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, tsurf_in, &
+       sst_nff, s_int)
 
     ! This subroutine cumulates some fields for each time-step during
@@ -671,5 +672,12 @@
     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
+    real, intent(in):: tsurf_in(:) ! (klon)
+    
+    REAL, intent(in):: sst_nff(:) ! (knon)
+    ! SST not used to compute surface fluxes, in K. If
+    ! activate_ocean_skin == 0 then it is not defined; if
+    ! activate_ocean_skin == 1 then it is the ocean-air interface
+    ! temperature; if activate_ocean_skin == 2 then it is the bulk SST.
+
     real, intent(in):: s_int(:) ! (knon) ocean-air interface salinity, in ppt
 
@@ -710,5 +718,5 @@
 
        if (activate_ocean_skin == 2) then
-          cpl_t_int(:knon) = 0.
+          cpl_delta_temp(:knon) = 0.
           cpl_s_int = 0.
        end if
@@ -756,5 +764,6 @@
 
        if (activate_ocean_skin == 2) then
-          cpl_t_int(ig) = cpl_t_int(ig) + t_int(ig) / REAL(nexca)
+          cpl_delta_temp(ig) = cpl_delta_temp(ig) &
+               + (tsurf_in(ig) - sst_nff(ig)) / REAL(nexca)
           cpl_s_int(ig) = cpl_s_int(ig) + s_int(ig) / REAL(nexca)
        end if
@@ -806,5 +815,5 @@
 
           if (activate_ocean_skin == 2) then
-             ALLOCATE(cpl_t_int_2D(nbp_lon, jj_nb), &
+             ALLOCATE(cpl_delta_temp_2D(nbp_lon, jj_nb), &
                   cpl_s_int_2D(nbp_lon, jj_nb), stat = error)
              sum_error = sum_error + error
@@ -866,5 +875,5 @@
             CALL gath2cpl(cpl_atm_co2(:,cpl_index), cpl_atm_co22D(:,:), knon, knindex)
        if (activate_ocean_skin == 2) then
-          CALL gath2cpl(cpl_t_int, cpl_t_int_2D, knon, knindex)
+          CALL gath2cpl(cpl_delta_temp, cpl_delta_temp_2D, knon, knindex)
           CALL gath2cpl(cpl_s_int, cpl_s_int_2D, knon, knindex)
        end if
@@ -1264,5 +1273,5 @@
 
     if (activate_ocean_skin == 2) then
-       tab_flds(:, :, ids_t_int) = cpl_t_int_2D
+       tab_flds(:, :, ids_delta_temp) = cpl_delta_temp_2D
        tab_flds(:, :, ids_s_int) = cpl_s_int_2D
     end if
@@ -1503,5 +1512,5 @@
     ENDIF
 
-    if (activate_ocean_skin == 2) deallocate(cpl_t_int_2d, cpl_s_int_2d)
+    if (activate_ocean_skin == 2) deallocate(cpl_delta_temp_2d, cpl_s_int_2d)
 
     IF (sum_error /= 0) THEN
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/oasis.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/oasis.F90	(revision 3739)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/oasis.F90	(revision 3740)
@@ -59,5 +59,5 @@
   INTEGER, PARAMETER :: ids_qraiic = 28
   INTEGER, PARAMETER :: ids_qsnoic = 29
-  INTEGER, PARAMETER :: ids_t_int = 30, ids_s_int = 31
+  INTEGER, PARAMETER :: ids_delta_temp = 30, ids_s_int = 31
   
   INTEGER, PARAMETER :: maxsend    = 31  ! Maximum number of fields to send
@@ -190,6 +190,6 @@
     
     if (activate_ocean_skin == 2) then
-       infosend(ids_t_int)%action = .TRUE.
-       infosend(ids_t_int)%name = 'COSSTSKN'
+       infosend(ids_delta_temp)%action = .TRUE.
+       infosend(ids_delta_temp)%name = 'CODTEMP'
        infosend(ids_s_int)%action = .TRUE.
        infosend(ids_s_int)%name = 'COSSSSKN'
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_cpl_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_cpl_mod.F90	(revision 3739)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/ocean_cpl_mod.F90	(revision 3740)
@@ -55,5 +55,5 @@
        radsol, snow, agesno, &
        qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l, sens_prec_liq, sss, s_int, rhoa)
+       tsurf_new, dflux_s, dflux_l, sens_prec_liq, sss, s_int, rhoa, sst_nff)
 
 !
@@ -91,7 +91,13 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: ps
     REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
-    REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
+    REAL, INTENT(IN) :: tsurf_in(:) ! (klon)
     real, intent(in):: s_int(:) ! (knon) ocean-air interface salinity, in ppt
     real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
+
+    REAL, intent(in):: sst_nff(:) ! (knon)
+    ! SST not used to compute surface fluxes, in K. If
+    ! activate_ocean_skin == 0 then it is not defined; if
+    ! activate_ocean_skin == 1 then it is the ocean-air interface
+    ! temperature; if activate_ocean_skin == 2 then it is the bulk SST.
 
 ! In/Output arguments
@@ -211,5 +217,5 @@
          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, tsurf_in, s_int)
+         lat_prec_sol, tsurf_in, sst_nff, s_int)
 
   END SUBROUTINE ocean_cpl_noice
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/pbl_surface_mod.F90	(revision 3739)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/pbl_surface_mod.F90	(revision 3740)
@@ -292,6 +292,6 @@
     USE print_control_mod,  ONLY : prt_level,lunout
     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: sst_nff, dter, dser, tkt, tks, taur, sss
+    use phys_state_var_mod, only: s_int, ds_ns, dt_ns, sst_nff
+    use phys_output_var_mod, only: dter, dser, tkt, tks, taur, sss
 #ifdef CPP_XIOS
     USE wxios, ONLY: missing_val
@@ -1429,6 +1429,9 @@
 
        if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
-          if (activate_ocean_skin == 2 .and. type_ocean == "couple") &
-               ys_int(:knon) = s_int(ni(:knon))
+          if (activate_ocean_skin == 2 .and. type_ocean == "couple") then
+             ys_int(:knon) = s_int(ni(:knon))
+             ysst_nff(:knon) = sst_nff(ni(:knon))
+          end if
+          
           yds_ns(:knon) = ds_ns(ni(:knon))
           ydt_ns(:knon) = dt_ns(ni(:knon))
@@ -3194,5 +3197,5 @@
 
     USE indice_sol_mod
-    use phys_state_var_mod, only: s_int, ds_ns, dt_ns
+    use phys_state_var_mod, only: s_int, ds_ns, dt_ns, sst_nff
     use config_ocean_skin_m, only: activate_ocean_skin
 
@@ -3284,5 +3287,9 @@
                    if (activate_ocean_skin >= 1) then
                       if (activate_ocean_skin == 2 &
-                           .and. type_ocean == "couple") s_int(i) = 35.
+                           .and. type_ocean == "couple") then
+                         s_int(i) = 35.
+                         sst_nff(i) = tsurf(i,nsrf)
+                      end if
+                      
                       ds_ns(i) = 0.
                       dt_ns(i) = 0.
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/phyetat0.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/phyetat0.F90	(revision 3739)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/phyetat0.F90	(revision 3740)
@@ -19,5 +19,5 @@
        wake_s, wake_dens, zgam, zmax0, zmea, zpic, zsig, &
        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m, treedrg, &
-       ale_wake, ale_bl_stat, s_int, ds_ns, dt_ns
+       ale_wake, ale_bl_stat, s_int, ds_ns, dt_ns, sst_nff
 !FC
   USE geometry_mod, ONLY : longitude_deg, latitude_deg
@@ -541,7 +541,11 @@
 
   if (activate_ocean_skin >= 1) then
-     if (activate_ocean_skin == 2 .and. type_ocean == 'couple') found &
-          = phyetat0_get(1, s_int, "S_int", "salinity at air-ocean interface", &
-          35.)
+     if (activate_ocean_skin == 2 .and. type_ocean == 'couple') then
+        found = phyetat0_get(1, s_int, "S_int", &
+             "salinity at air-ocean interface", 35.)
+        found = phyetat0_get(1, sst_nff, "T1", "bulk SST of coupled ocean", 0.)
+        if (.not. found) sst_nff = ftsol(:,is_oce)
+     end if
+     
      found = phyetat0_get(1, ds_ns, "dS_ns", "delta salinity near surface", 0.)
      found = phyetat0_get(1, dt_ns, "dT_ns", "delta temperature near surface", &
@@ -552,4 +556,5 @@
         ds_ns = missing_val
         dt_ns = missing_val
+        sst_nff = missing_val
      end where
   end if
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/phyredem.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/phyredem.F90	(revision 3739)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/phyredem.F90	(revision 3740)
@@ -27,5 +27,5 @@
                                 ale_wake, ale_bl_stat,                       &
                                 du_gwd_rando, du_gwd_front, u10m, v10m,      &
-                                treedrg, s_int, ds_ns, dt_ns
+                                treedrg, s_int, ds_ns, dt_ns, sst_nff
   USE geometry_mod, ONLY : longitude_deg, latitude_deg
   USE iostart, ONLY: open_restartphy, close_restartphy, enddef_restartphy, put_field, put_var
@@ -342,7 +342,10 @@
 
     if (activate_ocean_skin >= 1) then
-       if (activate_ocean_skin == 2 .and. type_ocean == 'couple') &
-            CALL put_field(pass, "S_int", "salinity at air-ocean interface", &
-            s_int)
+       if (activate_ocean_skin == 2 .and. type_ocean == 'couple') then
+          CALL put_field(pass, "S_int", "salinity at air-ocean interface", &
+               s_int)
+          CALL put_field(pass, "T1", "bulk SST of coupled ocean", sst_nff)
+       end if
+       
        CALL put_field(pass, "dS_ns", "delta salinity near surface", ds_ns)
        CALL put_field(pass, "dT_ns", "delta temperature near surface", dT_ns)
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_var_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_var_mod.F90	(revision 3739)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_var_mod.F90	(revision 3740)
@@ -135,10 +135,4 @@
   ! Ocean-atmosphere interface, subskin ocean and near-surface ocean:
   
-  REAL, ALLOCATABLE, SAVE:: sst_nff(:)
-  ! SST not used to compute surface fluxes, in K. If
-  ! activate_ocean_skin == 0 then it is not allocated; if
-  ! activate_ocean_skin == 1 then it is the ocean-air interface
-  ! temperature; if activate_ocean_skin == 2 then it is the bulk SST.
-
   REAL, ALLOCATABLE, SAVE:: dter(:)
   ! Temperature variation in the diffusive microlayer, that is
@@ -161,5 +155,5 @@
   ! bulk salinity of the surface layer of the ocean, in ppt
   
-  !$OMP THREADPRIVATE(sst_nff, dter, dser, tkt, tks, taur, sss)
+  !$OMP THREADPRIVATE(dter, dser, tkt, tks, taur, sss)
 
 CONTAINS
@@ -222,6 +216,6 @@
     IF (ok_gwd_rando) allocate(zustr_gwd_rando(klon), zvstr_gwd_rando(klon))
 
-    if (activate_ocean_skin >= 1) allocate(sst_nff(klon), dter(klon), &
-         dser(klon), tkt(klon), tks(klon), taur(klon), sss(klon))
+    if (activate_ocean_skin >= 1) allocate(dter(klon), dser(klon), tkt(klon), &
+         tks(klon), taur(klon), sss(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 3739)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/phys_output_write_mod.F90	(revision 3740)
@@ -249,5 +249,6 @@
          ulevSTD, vlevSTD, wlevSTD, philevSTD, qlevSTD, tlevSTD, &
          rhlevSTD, O3STD, O3daySTD, uvSTD, vqSTD, vTSTD, wqSTD, vphiSTD, &
-         wTSTD, u2STD, v2STD, T2STD, missing_val_nf90, s_int, ds_ns, dt_ns
+         wTSTD, u2STD, v2STD, T2STD, missing_val_nf90, s_int, ds_ns, dt_ns, &
+         sst_nff
 
     USE phys_local_var_mod, ONLY: zxfluxlat, slp, ptstar, pt0, zxtsol, zt2m, &
@@ -360,5 +361,5 @@
          alt_tropo, &
 !Ionela
-         ok_4xCO2atm, sst_nff, dter, dser, tkt, tks, taur, sss
+         ok_4xCO2atm, dter, dser, tkt, tks, taur, sss
 
     USE ocean_slab_mod, ONLY: nslay, tslab, slab_bilg, tice, seaice, &
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/phys_state_var_mod.F90	(revision 3739)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/phys_state_var_mod.F90	(revision 3740)
@@ -444,5 +444,11 @@
       ! minus foundation temperature. (Can be negative.) In K.
       
-      !$OMP THREADPRIVATE(s_int, ds_ns, dt_ns)
+      REAL, ALLOCATABLE, SAVE:: sst_nff(:)
+      ! SST not used to compute surface fluxes, in K. If
+      ! activate_ocean_skin == 0 then it is not allocated; if
+      ! activate_ocean_skin == 1 then it is the ocean-air interface
+      ! temperature; if activate_ocean_skin == 2 then it is the bulk SST.
+
+      !$OMP THREADPRIVATE(s_int, ds_ns, dt_ns, sst_nff)
 
     CONTAINS
@@ -653,5 +659,5 @@
       ENDIF
       if (activate_ocean_skin >= 1) ALLOCATE(s_int(klon), ds_ns(klon), &
-           dt_ns(klon))
+           dt_ns(klon), sst_nff(klon))
 
     END SUBROUTINE phys_state_var_init
@@ -791,5 +797,5 @@
 !!! fin nrlmd le 10/04/2012
 
-      if (activate_ocean_skin >= 1) deALLOCATE(s_int, ds_ns, dt_ns)
+      if (activate_ocean_skin >= 1) deALLOCATE(s_int, ds_ns, dt_ns, sst_nff)
 
       is_initialized=.FALSE.
Index: LMDZ6/branches/Ocean_skin/libf/phylmd/surf_ocean_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/surf_ocean_mod.F90	(revision 3739)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/surf_ocean_mod.F90	(revision 3740)
@@ -88,4 +88,10 @@
     ! minus foundation temperature. (Can be negative.) In K.
 
+    REAL, intent(inout):: sst_nff(:) ! (knon)
+    ! SST not used to compute surface fluxes, in K. If
+    ! activate_ocean_skin == 0 then it is not defined; if
+    ! activate_ocean_skin == 1 then it is the ocean-air interface
+    ! temperature; if activate_ocean_skin == 2 then it is the bulk SST.
+
     ! Output variables
     !******************************************************************************
@@ -103,10 +109,4 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
 
-    REAL, intent(out):: sst_nff(:) ! (knon)
-    ! SST not used to compute surface fluxes, in K. If
-    ! activate_ocean_skin == 0 then it is not defined; if
-    ! activate_ocean_skin == 1 then it is the ocean-air interface
-    ! temperature; if activate_ocean_skin == 2 then it is the bulk SST.
-
     REAL, intent(out):: dter(:) ! (knon)
     ! Temperature variation in the diffusive microlayer, that is
@@ -190,5 +190,6 @@
             radsol, snow, agesno, &
             qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-            tsurf_new, dflux_s, dflux_l, sens_prec_liq, sss, s_int, rhoa)
+            tsurf_new, dflux_s, dflux_l, sens_prec_liq, sss, s_int, rhoa, &
+            sst_nff)
 
     CASE('slab')
