Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/ce0l.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/ce0l.F90	(revision 3814)
+++ LMDZ6/trunk/libf/dynphy_lonlat/phylmd/ce0l.F90	(revision 3815)
@@ -29,11 +29,16 @@
   USE iniphysiq_mod,  ONLY: iniphysiq
   USE mod_const_mpi,  ONLY: comm_lmdz
+
 #ifdef CPP_PARA
   USE mod_const_mpi,  ONLY: init_const_mpi
-  USE parallel_lmdz,  ONLY: init_parallel, mpi_rank, omp_rank
+  USE parallel_lmdz,  ONLY: init_parallel, mpi_rank, omp_rank, using_mpi
   USE bands,          ONLY: read_distrib, distrib_phys
   USE mod_hallo,      ONLY: init_mod_hallo
   USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys
-#endif
+#ifdef CPP_XIOS
+  USE xios, only: xios_finalize
+#endif
+#endif
+
   USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, g, kappa, omeg, r, rad, &
                           pi, jmp1
@@ -50,4 +55,8 @@
   include "comgeom2.h"
   include "iniprint.h"
+#ifdef CPP_MPI
+  include 'mpif.h'
+#endif
+  
   REAL               :: masque(iip1,jjp1)             !--- CONTINENTAL MASK
   REAL               :: phis  (iip1,jjp1)             !--- GROUND GEOPOTENTIAL
@@ -67,5 +76,7 @@
   REAL, ALLOCATABLE  :: masktmp(:)
 
-#ifndef CPP_PARA
+#ifdef CPP_PARA
+  integer ierr
+#else
 ! for iniphysiq in serial mode
   INTEGER,PARAMETER :: mpi_rank=0
@@ -168,5 +179,7 @@
     ocemask = ocetmp
     IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
-      DO j=1,jjp1; ocemask(:,j) = ocetmp(:,jjp1-j+1); END DO
+       DO j=1,jjp1
+          ocemask(:,j) = ocetmp(:,jjp1-j+1)
+       END DO
     END IF
     DEALLOCATE(ocetmp,lon_omask,lat_omask,dlon_omask,dlat_omask)
@@ -239,4 +252,10 @@
 #ifdef CPP_PARA
   END IF
+#ifdef CPP_XIOS
+  CALL xios_finalize
+#endif
+#ifdef CPP_MPI
+  IF (using_mpi) call MPI_FINALIZE(ierr)
+#endif
 #endif
 
Index: LMDZ6/trunk/libf/misc/handle_err_m.F90
===================================================================
--- LMDZ6/trunk/libf/misc/handle_err_m.F90	(revision 3814)
+++ LMDZ6/trunk/libf/misc/handle_err_m.F90	(revision 3815)
@@ -27,4 +27,5 @@
 
     if (ncerr /= nf90_noerr) then
+       print *, "NetCDF95 handle_err:"
        print *, message, ":"
        if (present(varid)) print *, "varid = ", varid
@@ -39,5 +40,5 @@
           end if
        end if
-       call abort_physic("NetCDF95 handle_err", "", 1)
+       call abort_physic("NetCDF95 handle_err", "see above", 1)
     end if
 
Index: LMDZ6/trunk/libf/phylmd/YOMCST.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/YOMCST.h	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/YOMCST.h	(revision 3815)
@@ -22,5 +22,5 @@
       REAL RMO3,RMCO2,RMC,RMCH4,RMN2O,RMCFC11,RMCFC12
       REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV
-      REAL RKAPPA,RETV
+      REAL RKAPPA,RETV, eps_w
 ! A1.5,6 Thermodynamic liquid,solid phases
       REAL RCW,RCS
@@ -38,5 +38,5 @@
      &      ,R     ,RMD   ,RMV   ,RD    ,RV    ,RCPD                    &
      &      ,RMO3  ,RMCO2 ,RMC   ,RMCH4 ,RMN2O ,RMCFC11 ,RMCFC12        &
-     &      ,RCPV  ,RCVD  ,RCVV  ,RKAPPA,RETV                           &
+     &      ,RCPV  ,RCVD  ,RCVV  ,RKAPPA,RETV, eps_w                    &
      &      ,RCW   ,RCS                                                 &
      &      ,RLVTT ,RLSTT ,RLMLT ,RTT   ,RATM                           &
Index: LMDZ6/trunk/libf/phylmd/calcul_fluxs_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/calcul_fluxs_mod.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/calcul_fluxs_mod.F90	(revision 3815)
@@ -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.
@@ -259,18 +261,12 @@
 !
 ! calcul de l'enthalpie des precipitations liquides et solides
-!
-!       if (PRESENT(enth_prec_liq))                   &
-!       enth_prec_liq(i) = rcw * (t1lay(i) - tsurf(i)) * &
-!                          precip_rain(i)
-!       if (PRESENT(enth_prec_sol))                  &
-!       enth_prec_sol(i) = rcs * (t1lay(i) - tsurf(i)) * &
-!                          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.
+       ! On calcule par rapport a T=0
+       !! sens_prec_liq(i) = rcw * (t1lay(i) - RTT) * precip_rain(i)
+       !! sens_prec_sol(i) = rcs * (t1lay(i) - RTT) * precip_snow(i)
+
        if (PRESENT(lat_prec_liq))                    &
          lat_prec_liq(i) =  precip_rain(i) * (RLVTT - RLVTT) 
@@ -280,13 +276,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/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 3815)
@@ -29,4 +29,5 @@
     USE mod_grid_phy_lmdz, ONLY: klon_glo
     USE print_control_mod, ONLY: lunout
+    use config_ocean_skin_m, only: config_ocean_skin
     USE phys_state_var_mod, ONLY: phys_tstep
 
@@ -2820,6 +2821,5 @@
 
     !$OMP END MASTER
-
-    RETURN
+    call config_ocean_skin
 
   END SUBROUTINE conf_phys
Index: LMDZ6/trunk/libf/phylmd/cpl_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/cpl_mod.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/cpl_mod.F90	(revision 3815)
@@ -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_delta_sst(:), cpl_delta_sal(:)
+  !$OMP THREADPRIVATE(cpl_delta_sst, cpl_delta_sal)
+  
   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_delta_sst_2D(:,:), cpl_delta_sal_2D(:,:)
+  !$OMP THREADPRIVATE(cpl_delta_sst_2D, cpl_delta_sal_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
@@ -218,4 +232,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_delta_sst(klon), cpl_delta_sal(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
@@ -380,4 +405,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"
@@ -441,4 +467,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
 
@@ -498,5 +525,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
@@ -506,4 +534,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
@@ -515,4 +544,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
@@ -529,4 +562,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)
@@ -615,14 +649,18 @@
        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, delta_sst, &
+       delta_sal)
+
+    ! 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
@@ -636,6 +674,13 @@
     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):: delta_sst(:) ! (knon)
+    ! Ocean-air interface temperature minus bulk SST, in
+    ! K. Defined only if activate_ocean_skin >= 1.
+
+    real, intent(in):: delta_sal(:) ! (knon)
+    ! Ocean-air interface salinity minus bulk salinity, in ppt.
 
 ! Local variables
@@ -673,4 +718,9 @@
        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) then
+          cpl_delta_sst = 0.
+          cpl_delta_sal = 0.
+       end if
     ENDIF
        
@@ -714,4 +764,9 @@
 !!---OB: this is correct but why knindex ??
        ENDIF
+
+       if (activate_ocean_skin == 2) then
+          cpl_delta_sst(ig) = cpl_delta_sst(ig) + delta_sst(ig) / REAL(nexca)
+          cpl_delta_sal(ig) = cpl_delta_sal(ig) + delta_sal(ig) / REAL(nexca)
+       end if
      ENDDO
 
@@ -759,4 +814,10 @@
              sum_error = sum_error + error
           ENDIF
+
+          if (activate_ocean_skin == 2) then
+             ALLOCATE(cpl_delta_sst_2D(nbp_lon, jj_nb), &
+                  cpl_delta_sal_2D(nbp_lon, jj_nb), stat = error)
+             sum_error = sum_error + error
+          end if
 
           IF (sum_error /= 0) THEN
@@ -814,5 +875,9 @@
        IF (carbon_cycle_cpl) &
             CALL gath2cpl(cpl_atm_co2(:,cpl_index), cpl_atm_co22D(:,:), knon, knindex)
-   ENDIF
+       if (activate_ocean_skin == 2) then
+          CALL gath2cpl(cpl_delta_sst, cpl_delta_sst_2D, knon, knindex)
+          CALL gath2cpl(cpl_delta_sal, cpl_delta_sal_2D, knon, knindex)
+       end if
+    ENDIF
 
   END SUBROUTINE cpl_send_ocean_fields
@@ -850,5 +915,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
@@ -1149,4 +1214,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
 !    
@@ -1206,4 +1272,9 @@
     tab_flds(:,:,ids_qraiic) = cpl_sens_rain2D(:,:,2)
     tab_flds(:,:,ids_qsnoic) = cpl_sens_snow2D(:,:,2)
+
+    if (activate_ocean_skin == 2) then
+       tab_flds(:, :, ids_delta_sst) = cpl_delta_sst_2D
+       tab_flds(:, :, ids_delta_sal) = cpl_delta_sal_2D
+    end if
     
     IF (version_ocean=='nemo') THEN
@@ -1442,4 +1513,6 @@
     ENDIF
 
+    if (activate_ocean_skin == 2) deallocate(cpl_delta_sst_2d, cpl_delta_sal_2d)
+
     IF (sum_error /= 0) THEN
        abort_message='Pb in deallocation of cpl_xxxx2D coupling variables'
Index: LMDZ6/trunk/libf/phylmd/oasis.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/oasis.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/oasis.F90	(revision 3815)
@@ -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_delta_sst = 30, ids_delta_sal = 31
+  
+  INTEGER, PARAMETER :: maxsend    = 31  ! 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
   
 
@@ -110,6 +117,5 @@
     USE geometry_mod, ONLY: ind_cell_glo                    
     USE mod_phys_lmdz_mpi_data, ONLY: klon_mpi_para_nb
-    
-
+    use config_ocean_skin_m, only: activate_ocean_skin
 
 ! Local variables
@@ -183,4 +189,11 @@
     infosend(ids_calvin)%action = .TRUE. ; infosend(ids_calvin)%name = 'COCALVIN'
     
+    if (activate_ocean_skin == 2) then
+       infosend(ids_delta_sst)%action = .TRUE.
+       infosend(ids_delta_sst)%name = 'CODELSST'
+       infosend(ids_delta_sal)%action = .TRUE.
+       infosend(ids_delta_sal)%name = 'CODELSSS'
+    end if
+           
     IF (version_ocean=='nemo') THEN
         infosend(ids_shftot)%action = .TRUE. ; infosend(ids_shftot)%name = 'COQSRMIX'
@@ -219,4 +232,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 = 'SISUSALW'
+   end if
    
    IF (cpl_current ) THEN
Index: LMDZ6/trunk/libf/phylmd/ocean_cpl_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ocean_cpl_mod.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/ocean_cpl_mod.F90	(revision 3815)
@@ -52,8 +52,9 @@
        AcoefH, AcoefQ, BcoefH, BcoefQ, &
        AcoefU, AcoefV, BcoefU, BcoefV, &
-       ps, u1, v1, gustiness, &
+       ps, u1, v1, gustiness, tsurf_in, &
        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, sss, delta_sal, rhoa, &
+       delta_sst)
 
 !
@@ -63,9 +64,10 @@
 !
     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
 
     INCLUDE "YOMCST.h"
@@ -90,4 +92,14 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: ps
     REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
+    REAL, INTENT(IN) :: tsurf_in(:) ! (klon)
+
+    real, intent(in):: delta_sal(:) ! (knon)
+    ! ocean-air interface salinity minus bulk salinity, in ppt
+
+    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
+
+    REAL, intent(in):: delta_sst(:) ! (knon)
+    ! Ocean-air interface temperature minus bulk SST, in K. Defined
+    ! only if activate_ocean_skin >= 1.
 
 ! In/Output arguments
@@ -104,4 +116,8 @@
     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)
+    ! bulk salinity of the surface layer of the ocean, in ppt
   
 
@@ -116,5 +132,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    
 
@@ -128,5 +144,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)
 
 !****************************************************************************************
@@ -138,5 +155,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.
     
 
@@ -147,10 +164,20 @@
 
     CALL calcul_fluxs(knon, is_oce, dtime, &
-         tsurf_cpl, p1lay, cal, beta, cdragh, cdragq, ps, &
+         merge(tsurf_in, tsurf_cpl, activate_ocean_skin == 2), p1lay, cal, &
+         beta, cdragh, cdragq, ps, &
          precip_rain, precip_snow, snow, qsurf,  &
          radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
          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
+       ! 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
       i = knindex(j)
@@ -189,9 +216,8 @@
 !****************************************************************************************
 
-    CALL cpl_send_ocean_fields(itime, knon, knindex, &
-         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)
-    
+    CALL cpl_send_ocean_fields(itime, knon, knindex, 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, delta_sst, delta_sal)
 
   END SUBROUTINE ocean_cpl_noice
@@ -210,5 +236,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 
@@ -245,4 +271,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
@@ -272,5 +299,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    
 
@@ -280,5 +307,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.
 
 !****************************************************************************************
@@ -313,5 +340,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/trunk/libf/phylmd/ocean_forced_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 3815)
@@ -19,8 +19,8 @@
        AcoefH, AcoefQ, BcoefH, BcoefQ, &
        AcoefU, AcoefV, BcoefU, BcoefV, &
-       ps, u1, v1, gustiness, &
+       ps, u1, v1, gustiness, tsurf_in, &
        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
@@ -35,4 +35,5 @@
     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 config_ocean_skin_m, only: activate_ocean_skin
 
     INCLUDE "YOMCST.h"
@@ -53,4 +54,6 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: ps
     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
@@ -67,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
@@ -76,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    
 
@@ -122,5 +126,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
@@ -132,10 +136,12 @@
 ! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf
     CALL calcul_fluxs(knon, is_oce, dtime, &
-         tsurf_lim, p1lay, cal, beta, cdragh, cdragq, ps, &
+         merge(tsurf_in, tsurf_lim, activate_ocean_skin == 2), p1lay, cal, &
+         beta, cdragh, cdragq, ps, &
          precip_rain, precip_snow, snow, qsurf,  &
          radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
          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
 
     do j = 1, knon
@@ -167,5 +173,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. 
@@ -201,4 +207,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
@@ -230,5 +237,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    
 
@@ -264,5 +271,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
@@ -277,5 +284,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/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 3815)
@@ -27,4 +27,5 @@
 !!                                  wx_pbl_fuse, wx_pbl_split
                                   wx_pbl0_fuse, wx_pbl0_split
+  use config_ocean_skin_m, only: activate_ocean_skin
 
   IMPLICIT NONE
@@ -292,4 +293,11 @@
     USE print_control_mod,  ONLY : prt_level,lunout
     USE ioipsl_getin_p_mod, ONLY : getin_p
+    use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal
+    use phys_output_var_mod, only: dter, dser, tkt, tks, taur, sss
+#ifdef CPP_XIOS
+    USE wxios, ONLY: missing_val
+#else
+    use netcdf, only: missing_val => nf90_fill_real
+#endif
 
     IMPLICIT NONE
@@ -824,4 +832,9 @@
     ! Martin
 
+    real, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, ydser, &
+         ytkt, ytks, ytaur, ysss
+    ! compression of delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, tkt, tks,
+    ! taur, sss on ocean points
+
 !****************************************************************************************
 ! End of declarations
@@ -959,5 +972,4 @@
 !!! jyg le 10/02/2012
     rh2m_x(:) = 0. ; qsat2m_x(:) = 0. ; rh2m_w(:) = 0. ; qsat2m_w(:) = 0.
-!!!
 
 ! 2b) Initialization of all local variables that will be compressed later
@@ -1426,4 +1438,14 @@
           ENDDO
        ENDIF
+
+       if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
+          if (activate_ocean_skin == 2 .and. type_ocean == "couple") then
+             ydelta_sal(:knon) = delta_sal(ni(:knon))
+             ydelta_sst(:knon) = delta_sst(ni(:knon))
+          end if
+          
+          yds_ns(:knon) = ds_ns(ni(:knon))
+          ydt_ns(:knon) = dt_ns(ni(:knon))
+       end if
        
 !****************************************************************************************
@@ -1949,5 +1971,7 @@
                yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
                ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
-               y_flux_u1, y_flux_v1)
+               y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &
+               yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
+               ytkt(:knon), ytks(:knon), ytaur(:knon), ysss)
       IF (prt_level >=10) THEN
           print *,'arg de surf_ocean: ycdragh ',ycdragh
@@ -2563,4 +2587,28 @@
        ENDIF
 
+       if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
+          delta_sal = missing_val
+          ds_ns = missing_val
+          dt_ns = missing_val
+          delta_sst = missing_val
+          dter = missing_val
+          dser = missing_val
+          tkt = missing_val
+          tks = missing_val
+          taur = missing_val
+          sss = missing_val
+          
+          delta_sal(ni(:knon)) = ydelta_sal(:knon)
+          ds_ns(ni(:knon)) = yds_ns(:knon)
+          dt_ns(ni(:knon)) = ydt_ns(:knon)
+          delta_sst(ni(:knon)) = ydelta_sst(:knon)
+          dter(ni(:knon)) = ydter(:knon)
+          dser(ni(:knon)) = ydser(:knon)
+          tkt(ni(:knon)) = ytkt(:knon)
+          tks(ni(:knon)) = ytks(:knon)
+          taur(ni(:knon)) = ytaur(:knon)
+          sss(ni(:knon)) = ysss(:knon)
+       end if
+
 !****************************************************************************************
 ! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m 
@@ -3166,11 +3214,13 @@
 
 !albedo SB >>>
-SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
-     evap, z0m, z0h, agesno,                                  &
-     tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke)  
-!albedo SB <<<
+  SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
+       evap, z0m, z0h, agesno,                                  &
+       tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke)  
+    !albedo SB <<<
     ! Give default values where new fraction has appread
 
     USE indice_sol_mod
+    use phys_state_var_mod, only: delta_sal, ds_ns, dt_ns, delta_sst
+    use config_ocean_skin_m, only: activate_ocean_skin
 
     INCLUDE "dimsoil.h"
@@ -3259,4 +3309,14 @@
                       alb_dif(i,k,nsrf) = 0.06
                    ENDDO
+                   if (activate_ocean_skin >= 1) then
+                      if (activate_ocean_skin == 2 &
+                           .and. type_ocean == "couple") then
+                         delta_sal(i) = 0.
+                         delta_sst(i) = 0.
+                      end if
+                      
+                      ds_ns(i) = 0.
+                      dt_ns(i) = 0.
+                   end if
                 ELSE IF (nsrf.EQ.is_sic) THEN
                    tsurf(i,nsrf) = 271.35
Index: LMDZ6/trunk/libf/phylmd/phyetat0.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phyetat0.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/phyetat0.F90	(revision 3815)
@@ -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
+       ale_wake, ale_bl_stat, ds_ns, dt_ns, delta_sst, delta_sal
 !FC
   USE geometry_mod, ONLY : longitude_deg, latitude_deg
@@ -29,4 +29,10 @@
   USE ocean_slab_mod, ONLY: nslay, tslab, seaice, tice, ocean_slab_init
   USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
+#ifdef CPP_XIOS
+  USE wxios, ONLY: missing_val
+#else
+  use netcdf, only: missing_val => nf90_fill_real
+#endif
+  use config_ocean_skin_m, only: activate_ocean_skin
 
   IMPLICIT none
@@ -35,5 +41,4 @@
   ! Objet: Lecture de l'etat initial pour la physique
   !======================================================================
-  include "netcdf.inc"
   include "dimsoil.h"
   include "clesphys.h"
@@ -519,4 +524,24 @@
   ENDIF ! Slab        
 
+  if (activate_ocean_skin >= 1) then
+     if (activate_ocean_skin == 2 .and. type_ocean == 'couple') then
+        found = phyetat0_get(1, delta_sal, "delta_sal", &
+             "ocean-air interface salinity minus bulk salinity", 0.)
+        found = phyetat0_get(1, delta_sst, "delta_SST", &
+             "ocean-air interface temperature minus bulk SST", 0.)
+     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", &
+          0.)
+
+     where (pctsrf(:, is_oce) == 0.)
+        ds_ns = missing_val
+        dt_ns = missing_val
+        delta_sst = missing_val
+        delta_sal = missing_val
+     end where
+  end if
+
   ! on ferme le fichier
   CALL close_startphy
Index: LMDZ6/trunk/libf/phylmd/phyredem.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phyredem.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/phyredem.F90	(revision 3815)
@@ -26,6 +26,8 @@
                                 detr_therm, ale_bl, ale_bl_trig, alp_bl,     &
                                 ale_wake, ale_bl_stat,                       &
-                                du_gwd_rando, du_gwd_front, u10m, v10m,      &
-                                treedrg, solswfdiff
+                                du_gwd_rando, du_gwd_front, u10m, v10m, &
+                                treedrg, solswfdiff, delta_sal, ds_ns, dt_ns, &
+                                delta_sst
+
   USE geometry_mod, ONLY : longitude_deg, latitude_deg
   USE iostart, ONLY: open_restartphy, close_restartphy, enddef_restartphy, put_field, put_var
@@ -37,4 +39,5 @@
   USE ocean_slab_mod, ONLY : nslay, tslab, seaice, tice, fsic
   USE time_phylmdz_mod, ONLY: annee_ref, day_end, itau_phy, pdtphys
+  use config_ocean_skin_m, only: activate_ocean_skin  
 
   IMPLICIT none
@@ -342,7 +345,19 @@
          "tendency on zonal wind due to acama gravity waves", du_gwd_front)
 
+    if (activate_ocean_skin >= 1) then
+       if (activate_ocean_skin == 2 .and. type_ocean == 'couple') then
+          CALL put_field(pass, "delta_sal", &
+               "ocean-air interface salinity minus bulk salinity", delta_sal)
+          CALL put_field(pass, "delta_SST", &
+               "ocean-air interface temperature minus bulk SST", delta_sst)
+       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)
+    end if
+    
     IF (pass==1) CALL enddef_restartphy
     IF (pass==2) CALL close_restartphy
-  ENDDO
+ ENDDO
   
   !$OMP BARRIER
Index: LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3815)
@@ -1954,3 +1954,43 @@
 #endif
 
+   type(ctrl_out), save:: o_delta_sst &
+        = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'delta_SST', &
+        "ocean-air interface temperature minus bulk SST", "K", '')
+
+   type(ctrl_out), save:: o_delta_sal &
+        = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'delta_sal', &
+        "ocean-air interface salinity minus bulk salinity", "ppt", '')
+
+   type(ctrl_out), save:: o_ds_ns &
+        = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'dS_ns', &
+        "subskin salinity minus foundation salinity", "ppt", '')
+
+   type(ctrl_out), save:: o_dt_ns &
+        = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'dT_ns', &
+        "subskin temperature minus foundation temperature", "K", '')
+
+   type(ctrl_out), save:: o_dter &
+        = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'dTer', &
+        "ocean-air interface temperature minus sub-skin temperature", "K", '')
+
+   type(ctrl_out), save:: o_dser &
+        = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'dSer', &
+        "ocean-air interface salinity minus sub-skin salinity", "ppt", '')
+
+   type(ctrl_out), save:: o_tkt &
+        = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'tkt', &
+        "thickness of thermal microlayer", "m", '')
+
+   type(ctrl_out), save:: o_tks &
+        = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'tks', &
+        "thickness of salinity microlayer", "m", '')
+
+   type(ctrl_out), save:: o_taur &
+        = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'taur', &
+        "momentum flux due to rain", "Pa", '')
+
+   type(ctrl_out), save:: o_sss &
+        = ctrl_out([1, 10, 10, 1, 10, 10, 11, 11, 11, 11], 'SSS', &
+        "bulk sea-surface salinity", "ppt", '')
+
 END MODULE phys_output_ctrlout_mod
Index: LMDZ6/trunk/libf/phylmd/phys_output_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_var_mod.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/phys_output_var_mod.F90	(revision 3815)
@@ -133,4 +133,28 @@
  !$OMP THREADPRIVATE(sens_prec_liq_o, sens_prec_sol_o,lat_prec_liq_o,lat_prec_sol_o)
 
+  ! Ocean-atmosphere interface, subskin ocean and near-surface ocean:
+  
+  REAL, ALLOCATABLE, SAVE:: dter(:)
+  ! Temperature variation in the diffusive microlayer, that is
+  ! ocean-air interface temperature minus subskin temperature. In K.
+      
+  REAL, SAVE, ALLOCATABLE:: dser(:)
+  ! Temperature variation in the diffusive microlayer, that is
+  ! subskin temperature minus ocean-air interface temperature. In K.
+
+  REAL, SAVE, ALLOCATABLE:: tkt(:)
+  ! épaisseur (m) de la couche de diffusion thermique (microlayer)
+  ! cool skin thickness
+
+  REAL, SAVE, ALLOCATABLE:: tks(:)
+  ! épaisseur (m) de la couche de diffusion de masse (microlayer)
+  
+  REAL, SAVE, ALLOCATABLE:: taur(:) ! momentum flux due to rain, in Pa
+
+  REAL, SAVE, ALLOCATABLE:: sss(:)
+  ! bulk salinity of the surface layer of the ocean, in ppt
+  
+  !$OMP THREADPRIVATE(dter, dser, tkt, tks, taur, sss)
+
 CONTAINS
 
@@ -138,4 +162,5 @@
   SUBROUTINE phys_output_var_init
     use dimphy
+    use config_ocean_skin_m, only: activate_ocean_skin
 
     IMPLICIT NONE
@@ -191,4 +216,7 @@
     IF (ok_gwd_rando) allocate(zustr_gwd_rando(klon), zvstr_gwd_rando(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/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3815)
@@ -199,6 +199,7 @@
          o_col_O3_strato, o_col_O3_tropo,                 & 
 !--interactive CO2
-         o_flx_co2_ocean, o_flx_co2_land, o_flx_co2_ff, o_flx_co2_bb
-
+         o_flx_co2_ocean, o_flx_co2_land, o_flx_co2_ff, o_flx_co2_bb, &
+         o_delta_sst, o_delta_sal, o_ds_ns, o_dt_ns, o_dter, o_dser, o_tkt, &
+         o_tks, o_taur, o_sss
 
 #ifdef CPP_StratAer
@@ -245,6 +246,7 @@
          T2sumSTD, nlevSTD, du_gwd_rando, du_gwd_front, &
          ulevSTD, vlevSTD, wlevSTD, philevSTD, qlevSTD, tlevSTD, &
-         rhlevSTD, O3STD, O3daySTD, uvSTD, vqSTD, vTSTD, wqSTD, &
-         vphiSTD, wTSTD, u2STD, v2STD, T2STD, missing_val_nf90
+         rhlevSTD, O3STD, O3daySTD, uvSTD, vqSTD, vTSTD, wqSTD, vphiSTD, &
+         wTSTD, u2STD, v2STD, T2STD, missing_val_nf90, delta_sal, ds_ns, &
+         dt_ns, delta_sst
 
     USE phys_local_var_mod, ONLY: zxfluxlat, slp, ptstar, pt0, zxtsol, zt2m, &
@@ -353,5 +355,5 @@
          alt_tropo, &
 !Ionela
-         ok_4xCO2atm
+         ok_4xCO2atm, dter, dser, tkt, tks, taur, sss
 
     USE ocean_slab_mod, ONLY: nslay, tslab, slab_bilg, tice, seaice, &
@@ -380,4 +382,5 @@
 #endif
     USE tracinca_mod, ONLY: config_inca
+    use config_ocean_skin_m, only: activate_ocean_skin
 
     USE vertical_layers_mod, ONLY: presnivs
@@ -2385,4 +2388,17 @@
        ENDIF   !(iflag_phytrac==1)
 
+       if (activate_ocean_skin >= 1) then
+          CALL histwrite_phy(o_delta_sst, delta_sst)
+          CALL histwrite_phy(o_delta_sal, delta_sal)
+          CALL histwrite_phy(o_ds_ns, ds_ns)
+          CALL histwrite_phy(o_dt_ns, dt_ns)
+          CALL histwrite_phy(o_dter, dter)
+          CALL histwrite_phy(o_dser, dser)
+          CALL histwrite_phy(o_tkt, tkt)
+          CALL histwrite_phy(o_tks, tks)
+          CALL histwrite_phy(o_taur, taur)
+          CALL histwrite_phy(o_sss, sss)
+       end if
+
        IF (.NOT.vars_defined) THEN
           !$OMP MASTER
Index: LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3815)
@@ -430,5 +430,26 @@
 !$OMP THREADPRIVATE(is_initialized)    
 
-CONTAINS
+      ! Ocean-atmosphere interface:
+
+      REAL, ALLOCATABLE, SAVE:: ds_ns(:) ! (klon)
+      ! "delta salinity near surface". Salinity variation in the
+      ! near-surface turbulent layer. That is subskin salinity minus
+      ! foundation salinity. In ppt.
+
+      REAL, ALLOCATABLE, SAVE:: dt_ns(:) ! (klon)
+      ! "delta temperature near surface". Temperature variation in the
+      ! near-surface turbulent layer. That is subskin temperature
+      ! minus foundation temperature. (Can be negative.) In K.
+      
+      REAL, ALLOCATABLE, SAVE:: delta_sst(:) ! (klon)
+      ! Ocean-air interface temperature minus bulk SST, in
+      ! K. Allocated and defined only if activate_ocean_skin >= 1.
+
+      REAL, ALLOCATABLE, SAVE:: delta_sal(:) ! (klon)
+      ! Ocean-air interface salinity minus bulk salinity, in ppt
+      
+      !$OMP THREADPRIVATE(delta_sal, ds_ns, dt_ns, delta_sst)
+
+    CONTAINS
 
 !======================================================================
@@ -438,4 +459,5 @@
 USE infotrac_phy, ONLY : nbtr
 USE indice_sol_mod
+use config_ocean_skin_m, only: activate_ocean_skin
 IMPLICIT NONE
 
@@ -638,10 +660,15 @@
         du_gwd_front(:,:) = 0 !ym missing init   
       ENDIF
-END SUBROUTINE phys_state_var_init
+      if (activate_ocean_skin >= 1) ALLOCATE(delta_sal(klon), ds_ns(klon), &
+           dt_ns(klon), delta_sst(klon))
+
+    END SUBROUTINE phys_state_var_init
 
 !======================================================================
-SUBROUTINE phys_state_var_end
+    SUBROUTINE phys_state_var_end
+      ! Useful only for lmdz1d.
 !USE dimphy
 USE indice_sol_mod
+use config_ocean_skin_m, only: activate_ocean_skin
 IMPLICIT NONE
 include "clesphys.h"
@@ -771,5 +798,10 @@
       DEALLOCATE(ale_bl_trig)
 !!! fin nrlmd le 10/04/2012
+
+      if (activate_ocean_skin >= 1) deALLOCATE(delta_sal, ds_ns, dt_ns, &
+           delta_sst)
+
       is_initialized=.FALSE.
+      
 END SUBROUTINE phys_state_var_end
 
Index: LMDZ6/trunk/libf/phylmd/suphel.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/suphel.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/suphel.F90	(revision 3815)
@@ -131,4 +131,5 @@
   rcvv = rcpv - rv
   rkappa = rd/rcpd
+  eps_w = rmv / rmd
   retv = rv/rd - 1.
   WRITE (UNIT=6, FMT='('' *** Thermodynamic, gas     ***'')')
Index: LMDZ6/trunk/libf/phylmd/surf_ocean_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/surf_ocean_mod.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/surf_ocean_mod.F90	(revision 3815)
@@ -20,7 +20,9 @@
        z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, & 
        tsurf_new, dflux_s, dflux_l, lmt_bils, &
-       flux_u1, flux_v1)
+       flux_u1, flux_v1, delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, tkt, tks, &
+       taur, sss)
 
     use albedo, only: alboc, alboc_cd
+    use bulk_flux_m, only: bulk_flux
     USE dimphy, ONLY: klon, zmasq
     USE surface_data, ONLY     : type_ocean
@@ -30,4 +32,5 @@
     USE indice_sol_mod, ONLY : nbsrf, is_oce
     USE limit_read_mod
+    use config_ocean_skin_m, only: activate_ocean_skin
     !
     ! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force, 
@@ -50,8 +53,8 @@
     REAL, DIMENSION(klon), INTENT(IN)        :: lwnet  ! net longwave radiation at surface  
     REAL, DIMENSION(klon), INTENT(IN)        :: alb1   ! albedo in visible SW interval
-    REAL, DIMENSION(klon), INTENT(IN)        :: windsp
+    REAL, DIMENSION(klon), INTENT(IN)        :: windsp ! wind at 10 m, in m s-1
     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
@@ -73,6 +76,24 @@
     REAL, DIMENSION(klon), INTENT(inOUT):: z0h
 
+    REAL, intent(inout):: delta_sst(:) ! (knon)
+    ! Ocean-air interface temperature minus bulk SST, in K. Defined
+    ! only if activate_ocean_skin >= 1.
+
+    real, intent(inout):: delta_sal(:) ! (knon)
+    ! Ocean-air interface salinity minus bulk salinity, in ppt. Defined
+    ! only if activate_ocean_skin >= 1.
+
+    REAL, intent(inout):: ds_ns(:) ! (knon)
+    ! "delta salinity near surface". Salinity variation in the
+    ! near-surface turbulent layer. That is subskin salinity minus
+    ! foundation salinity. In ppt.
+
+    REAL, intent(inout):: dt_ns(:) ! (knon)
+    ! "delta temperature near surface". Temperature variation in the
+    ! near-surface turbulent layer. That is subskin temperature
+    ! minus foundation temperature. (Can be negative.) In K.
+
     ! Output variables
-    !******************************************************************************
+    !**************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT)       :: z0m
     !albedo SB >>>
@@ -83,11 +104,35 @@
     !albedo SB <<<     
     REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
-    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
+    REAL, INTENT(OUT):: tsurf_new(klon) ! sea surface temperature, in K
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
     REAL, DIMENSION(klon), INTENT(OUT)       :: lmt_bils
     REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
 
+    REAL, intent(out):: dter(:) ! (knon)
+    ! Temperature variation in the diffusive microlayer, that is
+    ! ocean-air interface temperature minus subskin temperature. In
+    ! K.
+
+    REAL, intent(out):: dser(:) ! (knon)
+    ! Salinity variation in the diffusive microlayer, that is
+    ! ocean-air interface salinity minus subskin salinity. In ppt.
+
+    REAL, intent(out):: tkt(:) ! (knon)
+    ! épaisseur (m) de la couche de diffusion thermique (microlayer)
+    ! cool skin thickness
+
+    REAL, intent(out):: tks(:) ! (knon)
+    ! épaisseur (m) de la couche de diffusion de masse (microlayer)
+
+    REAL, intent(out):: taur(:) ! (knon)
+    ! momentum flux due to rain, in Pa
+
+    real, intent(out):: 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.)
+
     ! Local variables
-    !******************************************************************************
+    !*************************************************************************
     INTEGER               :: i, k
     REAL                  :: tmp
@@ -97,7 +142,11 @@
     REAL, DIMENSION(klon) :: cdragq ! Cdrag pour l'evaporation
     CHARACTER(len=20),PARAMETER :: modname="surf_ocean"
-
-    ! End definition
-    !******************************************************************************
+    real rhoa(knon) ! density of moist air  (kg / m3)
+    REAL sens_prec_liq(knon)
+
+    REAL t_int(knon) ! ocean-air interface temperature, in K
+    real s_int(knon) ! ocean-air interface salinity, in ppt
+
+    !**************************************************************************
 
 
@@ -126,5 +175,5 @@
     ENDIF
 
-
+    rhoa = PS(:KNON) / (Rd * temp_air(:knon) * (1. + retv * spechum(:knon)))
     !******************************************************************************
     ! Switch according to type of ocean (couple, slab or forced)
@@ -139,8 +188,9 @@
             AcoefH, AcoefQ, BcoefH, BcoefQ, &
             AcoefU, AcoefV, BcoefU, BcoefV, &
-            ps, u1, v1, gustiness, &
+            ps, u1, v1, gustiness, tsurf_in, &
             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, sss, delta_sal, rhoa, &
+            delta_sst)
 
     CASE('slab')
@@ -162,8 +212,8 @@
             AcoefH, AcoefQ, BcoefH, BcoefQ, &
             AcoefU, AcoefV, BcoefU, BcoefV, &
-            ps, u1, v1, gustiness, &
+            ps, u1, v1, gustiness, tsurf_in, &
             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
 
@@ -268,8 +318,22 @@
        CALL abort_physic(modname,'version non prevue',1)
     ENDIF
-    !
-    !******************************************************************************
+
+    if (activate_ocean_skin >= 1) then
+       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_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 = [(rlvtt, i = 1, knon)], rf = - sens_prec_liq, dtime = dtime, &
+            rns = swnet(:knon))
+       delta_sst = t_int - tsurf_new(:knon)
+       delta_sal = s_int - sss(:knon)
+       if (activate_ocean_skin >= 2) tsurf_new(:knon) = t_int
+    end if
+    
   END SUBROUTINE surf_ocean
-  !******************************************************************************
+  !****************************************************************************
   !
 END MODULE surf_ocean_mod
Index: LMDZ6/trunk/libf/phylmd/surf_seaice_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/surf_seaice_mod.F90	(revision 3814)
+++ LMDZ6/trunk/libf/phylmd/surf_seaice_mod.F90	(revision 3815)
@@ -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
