Index: LMDZ6/branches/LMDZ-tracers/libf/phylmd/cpl_mod.F90
===================================================================
--- LMDZ6/branches/LMDZ-tracers/libf/phylmd/cpl_mod.F90	(revision 3725)
+++ LMDZ6/branches/LMDZ-tracers/libf/phylmd/cpl_mod.F90	(revision 3851)
@@ -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
@@ -153,4 +167,5 @@
     cpl_old_calving=.FALSE.
     CALL getin_p("cpl_old_calving",cpl_old_calving)
+    WRITE(lunout,*)' cpl_old_calving = ', cpl_old_calving
 
 
@@ -217,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
@@ -240,5 +266,4 @@
     ALLOCATE(cell_area2D(nbp_lon, jj_nb), stat = error)    
     sum_error = sum_error + error
-
 
     CALL gather_omp(longitude_deg,rlon_mpi)
@@ -251,4 +276,10 @@
       CALL Grid1DTo2D_mpi(rlat_mpi,lat2D)
       CALL Grid1DTo2D_mpi(cell_area_mpi,cell_area2D)
+      !--the next line is required for lat-lon grid and should have no impact
+      !--for an unstructured grid for which nbp_lon=1
+      !--if north pole in process mpi then divide cell area of pole cell by number of replicates
+      IF (is_north_pole_dyn) cell_area2D(:,1)=cell_area2D(:,1)/FLOAT(nbp_lon)
+      !--if south pole in process mpi then divide cell area of pole cell by number of replicates
+      IF (is_south_pole_dyn) cell_area2D(:,jj_nb)=cell_area2D(:,jj_nb)/FLOAT(nbp_lon)
       mask_calving(:,:,:) = 0 
       WHERE ( lat2D >= 40) mask_calving(:,:,1) = 1
@@ -278,5 +309,4 @@
     ENDIF
     
-           
     IF (sum_error /= 0) THEN
        abort_message='Pb allocation variables couplees'
@@ -349,5 +379,4 @@
     ENDIF    ! is_sequential
     
-
 !*************************************************************************************
 ! compatibility test
@@ -376,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"
@@ -437,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
 
@@ -494,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
@@ -502,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
@@ -511,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
@@ -525,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)
@@ -611,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
@@ -632,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
@@ -669,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
        
@@ -710,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
 
@@ -755,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
@@ -810,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
@@ -846,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
@@ -1145,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
 !    
@@ -1202,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
@@ -1438,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'
