Index: /LMDZ5/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 2655)
+++ /LMDZ5/trunk/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 2656)
@@ -9,5 +9,5 @@
                      nbp, communicator, &
                      punjours, pdayref,ptimestep, &
-                     rlatu,rlatv,rlonu,rlonv,aire,cu,cv,       &
+                     rlatudyn,rlatvdyn,rlonudyn,rlonvdyn,airedyn,cudyn,cvdyn, &
                      prad,pg,pr,pcpp,iflag_phys)
   USE dimphy, ONLY: init_dimphy
@@ -44,4 +44,6 @@
   USE mod_phys_lmdz_omp_data, ONLY: klon_omp 
 #endif
+  USE ioipsl_getin_p_mod, ONLY: getin_p
+  USE slab_heat_transp_mod, ONLY: ini_slab_transp_geom
   IMPLICIT NONE
 
@@ -52,6 +54,8 @@
 
   include "dimensions.h"
+  include "paramet.h"
   include "iniprint.h"
   include "tracstoke.h"
+  include "comgeom.h"
 
   REAL, INTENT (IN) :: prad ! radius of the planet (m)
@@ -65,11 +69,11 @@
   INTEGER, INTENT(IN) :: nbp ! number of physics columns for this MPI process
   INTEGER, INTENT(IN) :: communicator ! MPI communicator
-  REAL, INTENT (IN) :: rlatu(jj+1) ! latitudes of the physics grid
-  REAL, INTENT (IN) :: rlatv(jj) ! latitude boundaries of the physics grid
-  REAL, INTENT (IN) :: rlonv(ii+1) ! longitudes of the physics grid
-  REAL, INTENT (IN) :: rlonu(ii+1) ! longitude boundaries of the physics grid
-  REAL, INTENT (IN) :: aire(ii+1,jj+1) ! area of the dynamics grid (m2)
-  REAL, INTENT (IN) :: cu((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
-  REAL, INTENT (IN) :: cv((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
+  REAL, INTENT (IN) :: rlatudyn(jj+1) ! latitudes of the physics grid
+  REAL, INTENT (IN) :: rlatvdyn(jj) ! latitude boundaries of the physics grid
+  REAL, INTENT (IN) :: rlonvdyn(ii+1) ! longitudes of the physics grid
+  REAL, INTENT (IN) :: rlonudyn(ii+1) ! longitude boundaries of the physics grid
+  REAL, INTENT (IN) :: airedyn(ii+1,jj+1) ! area of the dynamics grid (m2)
+  REAL, INTENT (IN) :: cudyn((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
+  REAL, INTENT (IN) :: cvdyn((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
   INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation
   REAL, INTENT (IN) :: ptimestep !physics time step (s)
@@ -80,5 +84,8 @@
   CHARACTER (LEN=20) :: modname = 'iniphysiq'
   CHARACTER (LEN=80) :: abort_message
-
+  
+  LOGICAL :: slab_hdiff
+  INTEGER :: slab_ekman
+  CHARACTER (LEN = 6) :: type_ocean 
 
 #ifndef CPP_PARA
@@ -92,7 +99,7 @@
   CALL inigeomphy(ii,jj,nlayer, &
                nbp, communicator, &
-               rlatu,rlatv, &
-               rlonu,rlonv, &
-               aire,cu,cv)
+               rlatudyn,rlatvdyn, &
+               rlonudyn,rlonvdyn, &
+               airedyn,cudyn,cvdyn)
 
   ! --> now initialize things specific to the phylmd physics package
@@ -117,4 +124,18 @@
   ! Copy over "offline" settings
   CALL init_phystokenc(offline,istphy)
+
+  ! Initialization for slab heat transport
+  type_ocean="force"
+  CALL getin_p('type_ocean',type_ocean)
+  slab_hdiff=.FALSE.
+  CALL getin_p('slab_hdiff',slab_hdiff)
+  slab_ekman=0
+  CALL getin_p('slab_ekman',slab_ekman)
+  IF ((type_ocean=='slab').AND.(slab_hdiff.OR.(slab_ekman.GT.0))) THEN
+     CALL ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez,fext,unsaire,&
+                                  cu,cuvsurcv,cv,cvusurcu, &
+                                  aire,apoln,apols, &
+                                  aireu,airev) 
+  END IF
 
   ! Initialize tracer names, numbers, etc. for physics
@@ -159,5 +180,5 @@
 #ifdef INCA
      CALL init_inca_dim(klon_omp,nbp_lev,nbp_lon,nbp_lat - 1, &
-          rlonu,rlatu,rlonv,rlatv)
+          rlonudyn,rlatudyn,rlonvdyn,rlatvdyn)
 #endif
   END IF
Index: /LMDZ5/trunk/libf/phylmd/change_srf_frac_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/change_srf_frac_mod.F90	(revision 2655)
+++ /LMDZ5/trunk/libf/phylmd/change_srf_frac_mod.F90	(revision 2656)
@@ -87,6 +87,11 @@
        CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified)
     CASE ('slab')
+       IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN
+       ! Read fraction from limit.nc
+           CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified)
+       ELSE
        ! Get fraction from slab module
-       CALL ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified)
+           CALL ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified)
+       ENDIF
     CASE ('couple')
        ! Get fraction from the coupler
Index: /LMDZ5/trunk/libf/phylmd/dimphy.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/dimphy.F90	(revision 2655)
+++ /LMDZ5/trunk/libf/phylmd/dimphy.F90	(revision 2656)
@@ -9,7 +9,6 @@
   INTEGER,SAVE :: klevm1
   INTEGER,SAVE :: kflev
-  INTEGER,SAVE :: nslay
 
-!$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon,nslay)
+!$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon)
   REAL,save,allocatable,dimension(:) :: zmasq
 !$OMP THREADPRIVATE(zmasq)   
@@ -24,5 +23,4 @@
     
     klon=klon0
-    nslay=1 ! Slab, provisoire (F. Codron)
     kdlon=klon
     kidia=1
Index: /LMDZ5/trunk/libf/phylmd/limit_slab.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/limit_slab.F90	(revision 2655)
+++ /LMDZ5/trunk/libf/phylmd/limit_slab.F90	(revision 2656)
@@ -164,8 +164,4 @@
   diff_sst(:) = diff_sst_save(:)
   diff_siv(:) = diff_siv_save(:)
-! For Debug purpose
-!  PRINT *,'limit_slab sst',MINVAL(diff_sst(:)),MAXVAL(diff_sst(:))
-!  PRINT *,'limit_slab siv',MINVAL(diff_siv(:)),MAXVAL(diff_siv(:))
-!  PRINT *,'limit_slab bils',MINVAL(lmt_bils(:)),MAXVAL(lmt_bils(:))
 
 END SUBROUTINE limit_slab
Index: /LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 2655)
+++ /LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 2656)
@@ -9,4 +9,6 @@
   USE indice_sol_mod
   USE surface_data
+  USE mod_grid_phy_lmdz, ONLY: klon_glo
+  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
 
   IMPLICIT NONE
@@ -14,43 +16,68 @@
   PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
 
-!****************************************************************************************
+!***********************************************************************************
 ! Global saved variables
-!****************************************************************************************
-
+!***********************************************************************************
+  ! number of slab vertical layers
+  INTEGER, PUBLIC, SAVE :: nslay
+  !$OMP THREADPRIVATE(nslay)
+  ! timestep for coupling (update slab temperature) in timesteps
   INTEGER, PRIVATE, SAVE                           :: cpl_pas
   !$OMP THREADPRIVATE(cpl_pas)
+  ! cyang = 1/heat capacity of top layer (rho.c.H)
   REAL, PRIVATE, SAVE                              :: cyang
   !$OMP THREADPRIVATE(cyang)
+  ! depth of slab layers (1 or 2)
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
   !$OMP THREADPRIVATE(slabh)
+  ! slab temperature
   REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
   !$OMP THREADPRIVATE(tslab)
+  ! heat flux convergence due to Ekman
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_ekman
+  !$OMP THREADPRIVATE(dt_ekman)
+  ! heat flux convergence due to horiz diffusion
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_hdiff
+  !$OMP THREADPRIVATE(dt_hdiff)
+  ! fraction of ocean covered by sea ice (sic / (oce+sic))
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
   !$OMP THREADPRIVATE(fsic)
+  ! temperature of the sea ice
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
   !$OMP THREADPRIVATE(tice)
+  ! sea ice thickness, in kg/m2
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
   !$OMP THREADPRIVATE(seaice)
+  ! net surface heat flux, weighted by open ocean fraction
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
   !$OMP THREADPRIVATE(slab_bils)
+  ! slab_bils accumulated over cpl_pas timesteps
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
   !$OMP THREADPRIVATE(bils_cum)
+  ! net heat flux into the ocean below the ice : conduction + solar radiation
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bilg
   !$OMP THREADPRIVATE(slab_bilg)
+  ! slab_bilg over cpl_pas timesteps
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bilg_cum
   !$OMP THREADPRIVATE(bilg_cum)
-
-!****************************************************************************************
+  ! wind stress saved over cpl_pas timesteps
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: taux_cum
+  !$OMP THREADPRIVATE(taux_cum)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: tauy_cum
+  !$OMP THREADPRIVATE(tauy_cum)
+
+!***********************************************************************************
 ! Parameters (could be read in def file: move to slab_init)
-!****************************************************************************************
+!***********************************************************************************
 ! snow and ice physical characteristics:
     REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
     REAL, PARAMETER :: t_melt=273.15   ! melting ice temp
     REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
-    REAL, PARAMETER :: ice_den=917.
-    REAL, PARAMETER :: sea_den=1025.
-    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity
-    REAL, PARAMETER :: sno_cond=0.31*sno_den
+    REAL, PARAMETER :: ice_den=917. ! ice density
+    REAL, PARAMETER :: sea_den=1025. ! sea water density
+    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice
+    REAL, PARAMETER :: sno_cond=0.31*sno_den ! conductivity of snow
     REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice
+    REAL, PARAMETER :: sea_cap=3995.   ! specific heat capacity, snow and ice
     REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
 
@@ -74,36 +101,90 @@
     REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
     REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
-    REAL, PARAMETER :: pen_frac=0.3 !fraction of penetrating shortwave (no snow)
+    REAL, PARAMETER :: pen_frac=0.3 !fraction of shortwave penetrating into the
+    ! ice (no snow)
     REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
 
-!****************************************************************************************
+! horizontal transport
+   LOGICAL, PUBLIC, SAVE :: slab_hdiff
+   !$OMP THREADPRIVATE(slab_hdiff)
+   REAL, PRIVATE, SAVE    :: coef_hdiff ! coefficient for horizontal diffusion
+   !$OMP THREADPRIVATE(coef_hdiff)
+   INTEGER, PUBLIC, SAVE :: slab_ekman, slab_cadj ! Ekman, conv adjustment
+   !$OMP THREADPRIVATE(slab_ekman)
+   !$OMP THREADPRIVATE(slab_cadj)
+
+!***********************************************************************************
 
 CONTAINS
 !
-!****************************************************************************************
+!***********************************************************************************
 !
   SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
   !, seaice_rst etc
 
-    use IOIPSL
-
-    ! For ok_xxx vars (Ekman...)
-    INCLUDE "clesphys.h"
+    USE ioipsl_getin_p_mod, ONLY : getin_p
+    USE mod_phys_lmdz_transfert_para, ONLY : gather
+    USE slab_heat_transp_mod, ONLY : ini_slab_transp
 
     ! Input variables
-!****************************************************************************************
+!***********************************************************************************
     REAL, INTENT(IN)                         :: dtime
 ! Variables read from restart file
-    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst 
+    ! surface fractions from start file
 
 ! Local variables
-!****************************************************************************************
+!************************************************************************************
     INTEGER                :: error
+    REAL, DIMENSION(klon_glo) :: zmasq_glo
     CHARACTER (len = 80)   :: abort_message
     CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
 
-!****************************************************************************************
+!***********************************************************************************
+! Define some parameters
+!***********************************************************************************
+! Number of slab layers
+    nslay=2
+    CALL getin_p('slab_layers',nslay)
+    print *,'number of slab layers : ',nslay
+! Layer thickness
+    ALLOCATE(slabh(nslay), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation slabh'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+    slabh(1)=50.
+    IF (nslay.GT.1) THEN 
+        slabh(2)=150.
+    END IF
+
+! cyang = 1/heat capacity of top layer (rho.c.H)
+    cyang=1/(slabh(1)*sea_den*sea_cap)
+
+! cpl_pas  coupling period (update of tslab and ice fraction)
+! pour un calcul a chaque pas de temps, cpl_pas=1
+    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
+    CALL getin_p('cpl_pas',cpl_pas)
+    print *,'cpl_pas',cpl_pas
+
+! Horizontal diffusion
+    slab_hdiff=.FALSE.
+    CALL getin_p('slab_hdiff',slab_hdiff)
+    coef_hdiff=25000.
+    CALL getin_p('coef_hdiff',coef_hdiff)
+! Ekman transport
+    slab_ekman=0
+    CALL getin_p('slab_ekman',slab_ekman)
+! Convective adjustment
+    IF (nslay.EQ.1) THEN
+        slab_cadj=0
+    ELSE
+        slab_cadj=1
+    END IF
+    CALL getin_p('slab_cadj',slab_cadj)
+
+!************************************************************************************
 ! Allocate surface fraction read from restart file
-!****************************************************************************************
+!************************************************************************************
     ALLOCATE(fsic(klon), stat = error)
     IF (error /= 0) THEN
@@ -112,11 +193,12 @@
     ENDIF
     fsic(:)=0.
+    !zmasq = continent fraction
     WHERE (1.-zmasq(:)>EPSFRA)
         fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
     END WHERE
 
-!****************************************************************************************
-! Allocate saved variables
-!****************************************************************************************
+!************************************************************************************
+! Allocate saved fields
+!************************************************************************************
     ALLOCATE(tslab(klon,nslay), stat=error)
        IF (error /= 0) CALL abort_physic &
@@ -136,5 +218,5 @@
     bils_cum(:) = 0.0   
 
-    IF (version_ocean=='sicINT') THEN
+    IF (version_ocean=='sicINT') THEN ! interactive sea ice
         ALLOCATE(slab_bilg(klon), stat = error)
         IF (error /= 0) THEN
@@ -161,55 +243,69 @@
     END IF
 
-!****************************************************************************************
-! Define some parameters
-!****************************************************************************************
-! Layer thickness
-    ALLOCATE(slabh(nslay), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation slabh'
-       CALL abort_physic(modname,abort_message,1)
+    IF (slab_hdiff) THEN !horizontal diffusion
+        ALLOCATE(dt_hdiff(klon,nslay), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation dt_hdiff'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        dt_hdiff(:,:) = 0.0   
     ENDIF
-    slabh(1)=50.
-! cyang = 1/heat capacity of top layer (rho.c.H)
-    cyang=1/(slabh(1)*4.228e+06)
-
-! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
-! pour un calcul à chaque pas de temps, cpl_pas=1
-    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
-    CALL getin('cpl_pas',cpl_pas)
-    print *,'cpl_pas',cpl_pas
+
+    IF (slab_ekman.GT.0) THEN ! ekman transport
+        ALLOCATE(dt_ekman(klon,nslay), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation dt_ekman'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        dt_ekman(:,:) = 0.0   
+        ALLOCATE(taux_cum(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation taux_cum'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        taux_cum(:) = 0.0   
+        ALLOCATE(tauy_cum(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation tauy_cum'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        tauy_cum(:) = 0.0   
+    ENDIF
+
+! Initialize transport
+    IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
+      CALL gather(zmasq,zmasq_glo)
+!$OMP MASTER  ! Only master thread
+      IF (is_mpi_root) THEN 
+          CALL ini_slab_transp(zmasq_glo)
+      END IF
+!$OMP END MASTER
+    END IF
 
  END SUBROUTINE ocean_slab_init
 !
-!****************************************************************************************
+!***********************************************************************************
 !
   SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
 
-    USE limit_read_mod
-
-!    INCLUDE "clesphys.h"
+! this routine sends back the sea ice and ocean fraction to the main physics
+! routine. Called only with interactive sea ice
 
 ! Arguments
-!****************************************************************************************
-    INTEGER, INTENT(IN)                        :: itime   ! numero du pas de temps courant
-    INTEGER, INTENT(IN)                        :: jour    ! jour a lire dans l'annee
-    REAL   , INTENT(IN)                        :: dtime   ! pas de temps de la physique (en s)
+!************************************************************************************
+    INTEGER, INTENT(IN)                        :: itime   ! current timestep 
+    INTEGER, INTENT(IN)                        :: jour    !  day in year (not 
+    REAL   , INTENT(IN)                        :: dtime   ! physics timestep (s)
     REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg  ! sub-surface fraction
-    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step
-
-! Local variables
-!****************************************************************************************
-
-    IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN   
-       CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
-    ELSE
+    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is 
+                                                         ! modified at this time step
+
        pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
        pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
        is_modified=.TRUE.
-    END IF
 
   END SUBROUTINE ocean_slab_frac
 !
-!****************************************************************************************
+!************************************************************************************
 !
   SUBROUTINE ocean_slab_noice( & 
@@ -224,58 +320,83 @@
     
     USE calcul_fluxs_mod
+    USE slab_heat_transp_mod, ONLY: divgrad_phy, slab_ekman1, slab_ekman2
+    USE mod_phys_lmdz_para
 
     INCLUDE "clesphys.h"
 
+! This routine 
+! (1) computes surface turbulent fluxes over points with some open ocean    
+! (2) reads additional Q-flux (everywhere)
+! (3) computes horizontal transport (diffusion & Ekman)
+! (4) updates slab temperature every cpl_pas ; creates new ice if needed.
+
+! Note : 
+! klon total number of points
+! knon number of points with open ocean (varies with time)
+! knindex gives position of the knon points within klon.
+! In general, local saved variables have klon values
+! variables exchanged with PBL module have knon.
+
 ! Input arguments
-!****************************************************************************************
-    INTEGER, INTENT(IN)                  :: itime
-    INTEGER, INTENT(IN)                  :: jour
-    INTEGER, INTENT(IN)                  :: knon
-    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
-    REAL, INTENT(IN)                     :: dtime
-    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
-    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragq, cdragm
-    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
-    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
-    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
-    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
-    REAL, DIMENSION(klon), INTENT(IN)    :: ps
-    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
-    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
-    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
+!***********************************************************************************
+    INTEGER, INTENT(IN)                  :: itime ! current timestep INTEGER,
+    INTEGER, INTENT(IN)                  :: jour  ! day in year (for Q-Flux) 
+    INTEGER, INTENT(IN)                  :: knon  ! number of points 
+    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex 
+    REAL, INTENT(IN) :: dtime  ! timestep (s)
+    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay 
+    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragq, cdragm 
+    ! drag coefficients
+    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow 
+    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum ! near surface T, q
+    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ 
+    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV 
+    ! exchange coefficients for boundary layer scheme
+    REAL, DIMENSION(klon), INTENT(IN)    :: ps  ! surface pressure
+    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness ! surface wind
+    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in ! surface temperature
+    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol ! net surface radiative flux
 
 ! In/Output arguments
-!****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(INOUT) :: snow
+!************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2
     
 ! Output arguments
-!****************************************************************************************
+!************************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
     REAL, DIMENSION(klon), INTENT(OUT)   :: evap, fluxsens, fluxlat
     REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1, flux_v1
-    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
+    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new ! new surface tempearture
     REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l      
     REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
 
 ! Local variables
-!****************************************************************************************
-    INTEGER               :: i,ki
-    ! surface fluxes
+!************************************************************************************
+    INTEGER               :: i,ki,k
+    REAL                  :: t_cadj
+    !  for surface heat fluxes
     REAL, DIMENSION(klon) :: cal, beta, dif_grnd
-    ! flux correction
+    ! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes 
     REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
-    ! surface wind stress
+    ! for surface wind stress
     REAL, DIMENSION(klon) :: u0, v0
     REAL, DIMENSION(klon) :: u1_lay, v1_lay
-    ! ice creation
+    ! for new ice creation
     REAL                  :: e_freeze, h_new, dfsic
-
-!****************************************************************************************
-! 1) Flux calculation
-!
-!****************************************************************************************
-    cal(:)      = 0.
-    beta(:)     = 1.
-    dif_grnd(:) = 0.
+    ! horizontal diffusion and Ekman local vars
+    ! dimension = global domain (klon_glo) instead of // subdomains 
+    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo, dt_ekman_glo
+    ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop
+    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp
+    REAL, DIMENSION(klon_glo,nslay) :: tslab_glo
+    REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo
+
+!****************************************************************************************
+! 1) Surface fluxes calculation
+!    
+!****************************************************************************************
+    cal(:)      = 0. ! infinite thermal inertia
+    beta(:)     = 1. ! wet surface
+    dif_grnd(:) = 0. ! no diffusion into ground
     
 ! Suppose zero surface speed
@@ -285,4 +406,5 @@
     v1_lay(:) = v1(:) - v0(:)
 
+! Compute latent & sensible fluxes
     CALL calcul_fluxs(knon, is_oce, dtime, &
          tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, &
@@ -292,12 +414,6 @@
          tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
 
-! - Flux calculation at first modele level for U and V
-    CALL calcul_flux_wind(knon, dtime, &
-         u0, v0, u1, v1, gustiness, cdragm, &
-         AcoefU, AcoefV, BcoefU, BcoefV, &
-         p1lay, temp_air, &
-         flux_u1, flux_v1)  
-
-! Accumulate total fluxes locally
+! save total cumulated heat fluxes locally
+! radiative + turbulent + melt of falling snow 
     slab_bils(:)=0.
     DO i=1,knon
@@ -306,13 +422,28 @@
                       -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki)))
         bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
-! Also taux, tauy, saved vars...
     END DO
 
-!****************************************************************************************
-! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
+!  Compute surface wind stress
+    CALL calcul_flux_wind(knon, dtime, &
+         u0, v0, u1, v1, gustiness, cdragm, &
+         AcoefU, AcoefV, BcoefU, BcoefV, &
+         p1lay, temp_air, &
+         flux_u1, flux_v1)  
+
+! save cumulated wind stress
+    IF (slab_ekman.GT.0) THEN 
+      DO i=1,knon
+          ki=knindex(i)
+          taux_cum(ki)=taux_cum(ki)+flux_u1(i)*(1.-fsic(ki))/cpl_pas
+          tauy_cum(ki)=tauy_cum(ki)+flux_v1(i)*(1.-fsic(ki))/cpl_pas
+      END DO
+    ENDIF
+
+!****************************************************************************************
+! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file limit_slab.nc
 !
 !****************************************************************************************
     lmt_bils(:)=0.
-    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
+    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) 
     ! lmt_bils and diff_sst,siv saved by limit_slab
     qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
@@ -320,71 +451,160 @@
 
 !****************************************************************************************
-! 3) Recalculate new temperature
-!
+! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport)
+!    Bring to freezing temp and make sea ice if necessary
+!  
 !***********************************************o*****************************************
     tsurf_new=tsurf_in
     IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
-        ! Compute transport
-        ! Add read QFlux and SST tendency
-        tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas 
-        ! Add cumulated surface fluxes
-        tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
-        ! Update surface temperature
-        SELECT CASE(version_ocean)
-        CASE('sicNO')
-            DO i=1,knon
-                ki=knindex(i)
-                tsurf_new(i)=tslab(ki,1)
+! ***********************************
+!  Horizontal transport 
+! ***********************************
+      IF (slab_ekman.GT.0) THEN 
+          ! copy wind stress to global var
+          CALL gather(taux_cum,taux_glo)
+          CALL gather(tauy_cum,tauy_glo)
+      END IF
+
+      IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN 
+        CALL gather(tslab,tslab_glo)
+      ! Compute horiz transport on one process only
+!$OMP MASTER  ! Only master thread
+        IF (is_mpi_root) THEN ! Only master processus          
+          IF (slab_hdiff) THEN 
+              dt_hdiff_glo(:,:)=0.
+          END IF
+          IF (slab_ekman.GT.0) THEN 
+              dt_ekman_glo(:,:)=0.
+          END IF
+          DO i=1,cpl_pas ! time splitting for numerical stability
+            IF (slab_ekman.GT.0) THEN
+              SELECT CASE (slab_ekman)
+                CASE (1)
+                  CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
+                CASE (2)
+                  CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
+                CASE DEFAULT
+                  dt_ekman_tmp(:,:)=0.
+              END SELECT
+              dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:)
+              ! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1   
+              DO k=1,nslay
+                dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den)
+              ENDDO
+              tslab_glo=tslab_glo+dt_ekman_tmp*dtime
+            ENDIF
+            IF (slab_hdiff) THEN ! horizontal diffusion
+              ! laplacian of slab T
+              CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp)
+              ! multiply by diff coef and normalize to 50m slab equivalent
+              dt_hdiff_tmp=dt_hdiff_tmp*coef_hdiff*50./SUM(slabh)
+              dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:)
+              tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
+            END IF
+          END DO ! time splitting
+          IF (slab_hdiff) THEN
+            !dt_hdiff_glo saved in W/m2
+            DO k=1,nslay
+              dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)*slabh(k)*sea_den*sea_cap/cpl_pas
             END DO
-        CASE('sicOBS') ! check for sea ice or tslab below freezing
-            DO i=1,knon
-                ki=knindex(i)
-                IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
-                    tslab(ki,1)=t_freeze
-                END IF
-                tsurf_new(i)=tslab(ki,1)
-            END DO
-        CASE('sicINT')
-            DO i=1,knon
-                ki=knindex(i)
-                IF (fsic(ki).LT.epsfra) THEN ! Free of ice
-                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
-                        ! quantity of new ice formed
-                        e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
-                        ! new ice
-                        tice(ki)=t_freeze
-                        fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
-                        IF (fsic(ki).GT.ice_frac_min) THEN
-                            seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
-                            tslab(ki,1)=t_freeze
-                        ELSE
-                            fsic(ki)=0.
-                        END IF
-                        tsurf_new(i)=t_freeze
-                    ELSE
-                        tsurf_new(i)=tslab(ki,1)
-                    END IF 
-                ELSE ! ice present
-                    tsurf_new(i)=t_freeze
-                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
-                        ! quantity of new ice formed over open ocean
-                        e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
-                                 /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
-                        ! new ice height and fraction
-                        h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
-                        dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
-                        h_new=MIN(e_freeze/dfsic,h_ice_max)
-                        ! update tslab to freezing over open ocean only
-                        tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
-                        ! update sea ice
-                        seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
-                                   /(dfsic+fsic(ki))
-                        fsic(ki)=fsic(ki)+dfsic
-                        ! update snow?
-                    END IF !freezing
-                END IF ! sea ice present
-            END DO
-        END SELECT
-        bils_cum(:)=0.0! clear cumulated fluxes
+          END IF
+          IF (slab_ekman.GT.0) THEN 
+            ! dt_ekman_glo saved in W/m2
+            dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas
+          END IF
+        END IF ! is_mpi_root
+!$OMP END MASTER
+!$OMP BARRIER
+        ! Send new fields back to all processes
+        CALL Scatter(tslab_glo,tslab)
+        IF (slab_hdiff) THEN
+          CALL Scatter(dt_hdiff_glo,dt_hdiff)
+        END IF
+        IF (slab_ekman.GT.0) THEN 
+          CALL Scatter(dt_ekman_glo,dt_ekman)
+          ! clear wind stress
+          taux_cum(:)=0.
+          tauy_cum(:)=0.
+        END IF
+      ENDIF ! transport
+
+! ***********************************
+! Other heat fluxes
+! ***********************************
+      ! Add read QFlux
+      tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas 
+      ! Add cumulated surface fluxes
+      tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
+      ! Convective adjustment if 2 layers 
+      IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN
+        DO i=1,klon
+          IF (tslab(i,2).GT.tslab(i,1)) THEN
+            ! mean (mass-weighted) temperature
+            t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:))
+            tslab(i,1)=t_cadj
+            tslab(i,2)=t_cadj
+          END IF
+        END DO
+      END IF
+! ***********************************
+! Update surface temperature and ice 
+! ***********************************
+      SELECT CASE(version_ocean)
+      CASE('sicNO') ! no sea ice even below freezing !
+          DO i=1,knon
+              ki=knindex(i)
+              tsurf_new(i)=tslab(ki,1)
+          END DO
+      CASE('sicOBS') ! "realistic" case, for prescribed sea ice
+        ! tslab cannot be below freezing, or above it if there is sea ice
+          DO i=1,knon
+              ki=knindex(i)
+              IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
+                  tslab(ki,1)=t_freeze
+              END IF
+              tsurf_new(i)=tslab(ki,1)
+          END DO
+      CASE('sicINT') ! interactive sea ice
+          DO i=1,knon
+              ki=knindex(i)
+              IF (fsic(ki).LT.epsfra) THEN ! Free of ice
+                  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
+                      ! quantity of new ice formed
+                      e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
+                      ! new ice
+                      tice(ki)=t_freeze
+                      fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
+                      IF (fsic(ki).GT.ice_frac_min) THEN
+                          seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
+                          tslab(ki,1)=t_freeze
+                      ELSE
+                          fsic(ki)=0.
+                      END IF
+                      tsurf_new(i)=t_freeze
+                  ELSE
+                      tsurf_new(i)=tslab(ki,1)
+                  END IF 
+              ELSE ! ice present
+                  tsurf_new(i)=t_freeze
+                  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
+                      ! quantity of new ice formed over open ocean
+                      e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
+                               /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
+                      ! new ice height and fraction
+                      h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
+                      dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
+                      h_new=MIN(e_freeze/dfsic,h_ice_max)
+                      ! update tslab to freezing over open ocean only
+                      tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
+                      ! update sea ice
+                      seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
+                                 /(dfsic+fsic(ki))
+                      fsic(ki)=fsic(ki)+dfsic
+                      ! update snow?
+                  END IF ! tslab below freezing
+              END IF ! sea ice present
+          END DO
+      END SELECT
+      bils_cum(:)=0.0! clear cumulated fluxes
     END IF ! coupling time
   END SUBROUTINE ocean_slab_noice
@@ -530,7 +750,7 @@
            seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
         ENDIF
-! Melt / Freeze from above if Tsurf>0
+! Melt / Freeze snow from above if Tsurf>0
         IF (tsurf_new(i).GT.t_melt) THEN
-            ! energy available for melting snow (in kg/m2 of snow)
+            ! energy available for melting snow (in kg of melted snow /m2)
             e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
                /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
@@ -672,5 +892,5 @@
     END IF ! coupling time
     
-    !tests fraction glace
+    !tests ice fraction 
     WHERE (fsic.LT.ice_frac_min)
         tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
@@ -691,9 +911,14 @@
     IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
     IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
+    IF (ALLOCATED(tice)) DEALLOCATE(tice)
+    IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
     IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
     IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
     IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
     IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
-    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
+    IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum)
+    IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum)
+    IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman)
+    IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff)
 
   END SUBROUTINE ocean_slab_final
Index: /LMDZ5/trunk/libf/phylmd/phyetat0.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phyetat0.F90	(revision 2655)
+++ /LMDZ5/trunk/libf/phylmd/phyetat0.F90	(revision 2656)
@@ -3,5 +3,5 @@
 SUBROUTINE phyetat0 (fichnom, clesphy0, tabcntr0)
 
-  USE dimphy, only: klon, zmasq, klev, nslay
+  USE dimphy, only: klon, zmasq, klev
   USE iophy, ONLY : init_iophy_new
   USE ocean_cpl_mod,    ONLY : ocean_cpl_init
@@ -25,5 +25,5 @@
   USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
   USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
-  USE ocean_slab_mod, ONLY: tslab, seaice, tice, ocean_slab_init
+  USE ocean_slab_mod, ONLY: nslay, tslab, seaice, tice, ocean_slab_init
   USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
 
@@ -442,5 +442,15 @@
   IF ( type_ocean == 'slab' ) THEN
       CALL ocean_slab_init(dtime, pctsrf)
-      found=phyetat0_get(nslay,tslab,"tslab","tslab",0.)
+      IF (nslay.EQ.1) THEN
+        found=phyetat0_get(1,tslab,"tslab01","tslab",0.)
+        IF (.NOT. found) THEN
+            found=phyetat0_get(1,tslab,"tslab","tslab",0.)
+        END IF
+      ELSE
+          DO i=1,nslay
+            WRITE(str2,'(i2.2)') i
+            found=phyetat0_get(1,tslab(:,i),"tslab"//str2,"tslab",0.)  
+          END DO
+      END IF
       IF (.NOT. found) THEN 
           PRINT*, "phyetat0: Le champ <tslab> est absent"
@@ -452,6 +462,6 @@
 
       ! Sea ice variables
-      found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
       IF (version_ocean == 'sicINT') THEN
+          found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
           IF (.NOT. found) THEN 
               PRINT*, "phyetat0: Le champ <tice> est absent"
@@ -459,4 +469,5 @@
                   tice(:)=ftsol(:,is_sic)
           END IF 
+          found=phyetat0_get(1,seaice,"seaice","seaice",0.)
           IF (.NOT. found) THEN
               PRINT*, "phyetat0: Le champ <seaice> est absent"
Index: /LMDZ5/trunk/libf/phylmd/phyredem.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phyredem.F90	(revision 2655)
+++ /LMDZ5/trunk/libf/phylmd/phyredem.F90	(revision 2656)
@@ -30,5 +30,5 @@
   USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic, epsfra
   USE surface_data, ONLY: type_ocean, version_ocean
-  USE ocean_slab_mod, ONLY : tslab, seaice, tice, fsic
+  USE ocean_slab_mod, ONLY : nslay, tslab, seaice, tice, fsic
   USE time_phylmdz_mod, ONLY: annee_ref, day_end, itau_phy, pdtphys
 
@@ -58,5 +58,5 @@
 
   INTEGER isoil, nsrf,isw
-  CHARACTER (len=7) :: str7
+  CHARACTER (len=2) :: str2
   CHARACTER (len=256) :: nam, lnam
   INTEGER           :: it, iiq
@@ -303,5 +303,12 @@
   ! Restart variables for Slab ocean
   IF (type_ocean == 'slab') THEN
-      CALL put_field("tslab", "Slab ocean temperature", tslab)
+      IF (nslay.EQ.1) THEN
+        CALL put_field("tslab", "Slab ocean temperature", tslab)
+      ELSE
+        DO it=1,nslay
+          WRITE(str2,'(i2.2)') it
+          CALL put_field("tslab"//str2, "Slab ocean temperature", tslab(:,it))
+        END DO
+      END IF
       IF (version_ocean == 'sicINT') THEN
           CALL put_field("seaice", "Slab seaice (kg/m2)", seaice)
Index: /LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2655)
+++ /LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2656)
@@ -677,4 +677,8 @@
   TYPE(ctrl_out), SAVE :: o_slab_sic = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     'seaice', 'Epaisseur banquise slab', 'kg/m2', (/ ('', i=1, 9) /))
+  TYPE(ctrl_out), SAVE :: o_slab_hdiff = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
+    'slab_hdiff', 'Horizontal diffusion', 'W/m2', (/ ('', i=1, 9) /))
+  TYPE(ctrl_out), SAVE :: o_slab_ekman = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
+    'slab_ekman', 'Ekman heat transport', 'W/m2', (/ ('', i=1, 9) /))
   TYPE(ctrl_out), SAVE :: o_ale_bl = ctrl_out((/ 1, 1, 1, 10, 10, 10, 11, 11, 11 /), &
     'ale_bl', 'ALE BL', 'm2/s2', (/ ('', i=1, 9) /))
Index: /LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 2655)
+++ /LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 2656)
@@ -24,5 +24,5 @@
     ! defined and initialised in phys_output_mod.F90
 
-    USE dimphy, only: klon, klev, klevp1, nslay
+    USE dimphy, only: klon, klev, klevp1
     USE mod_phys_lmdz_para, ONLY: is_north_pole_phy,is_south_pole_phy
     USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
@@ -89,4 +89,5 @@
          o_slab_qflux, o_tslab, o_slab_bils, &
          o_slab_bilg, o_slab_sic, o_slab_tice, &
+         o_slab_hdiff, o_slab_ekman, &
          o_weakinv, o_dthmin, o_cldtau, &
          o_cldemi, o_pr_con_l, o_pr_con_i, &
@@ -291,5 +292,6 @@
   
 
-    USE ocean_slab_mod, only: tslab, slab_bils, slab_bilg, tice, seaice
+    USE ocean_slab_mod, only: nslay, tslab, slab_bils, slab_bilg, tice, &
+        seaice, slab_ekman,slab_hdiff, dt_ekman, dt_hdiff
     USE pbl_surface_mod, only: snow
     USE indice_sol_mod, only: nbsrf
@@ -982,5 +984,5 @@
               CALL histwrite_phy(o_tslab, zx_tmp_fi2d)
           ELSE
-              CALL histwrite_phy(o_tslab, tslab)
+              CALL histwrite_phy(o_tslab, tslab(:,1:nslay))
           END IF
           IF (version_ocean=='sicINT') THEN
@@ -988,4 +990,20 @@
               CALL histwrite_phy(o_slab_tice, tice)
               CALL histwrite_phy(o_slab_sic, seaice)
+          END IF
+          IF (slab_hdiff) THEN
+            IF (nslay.EQ.1) THEN
+                zx_tmp_fi2d(:)=dt_hdiff(:,1)
+                CALL histwrite_phy(o_slab_hdiff, zx_tmp_fi2d)
+            ELSE
+                CALL histwrite_phy(o_slab_hdiff, dt_hdiff(:,1:nslay))
+            END IF
+          END IF
+          IF (slab_ekman.GT.0) THEN
+            IF (nslay.EQ.1) THEN
+                zx_tmp_fi2d(:)=dt_ekman(:,1)
+                CALL histwrite_phy(o_slab_ekman, zx_tmp_fi2d)
+            ELSE
+                CALL histwrite_phy(o_slab_ekman, dt_ekman(:,1:nslay))
+            END IF
           END IF
        ENDIF !type_ocean == force/slab
Index: /LMDZ5/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2655)
+++ /LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2656)
@@ -1949,9 +1949,9 @@
           ! Particular case with annual mean insolation
           zzz=real(90) ! could be revisited
-          IF (read_climoz/=-1) THEN
-             abort_message ='read_climoz=-1 is recommended when ' &
-                  // 'solarlong0=1000.'
-             CALL abort_physic (modname,abort_message,1)
-          ENDIF
+!          IF (read_climoz/=-1) THEN
+!             abort_message ='read_climoz=-1 is recommended when ' &
+!                  // 'solarlong0=1000.'
+!             CALL abort_physic (modname,abort_message,1)
+!          ENDIF
        ELSE
           ! Case where the season is imposed with solarlong0
Index: /LMDZ5/trunk/libf/phylmd/slab_heat_transp_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/slab_heat_transp_mod.F90	(revision 2656)
+++ /LMDZ5/trunk/libf/phylmd/slab_heat_transp_mod.F90	(revision 2656)
@@ -0,0 +1,757 @@
+!
+MODULE slab_heat_transp_mod
+!
+! Slab ocean : temperature tendencies due to horizontal diffusion 
+! and / or Ekman transport
+
+USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo
+IMPLICIT NONE
+
+  ! Variables copied over from dyn3d dynamics:
+  REAL,SAVE,ALLOCATABLE :: fext(:) ! Coriolis f times cell area
+  !$OMP THREADPRIVATE(fext)
+  REAL,SAVE,ALLOCATABLE :: unsairez(:) ! 1/cell area
+  !$OMP THREADPRIVATE(unsairez)
+  REAL,SAVE,ALLOCATABLE :: unsaire(:)
+  !$OMP THREADPRIVATE(unsaire)
+  REAL,SAVE,ALLOCATABLE :: cu(:) ! cell longitude dim (m)
+  !$OMP THREADPRIVATE(cu)
+  REAL,SAVE,ALLOCATABLE :: cv(:) ! cell latitude dim (m)
+  !$OMP THREADPRIVATE(cv)
+  REAL,SAVE,ALLOCATABLE :: cuvsurcv(:) ! cu/cv (v points)
+  !$OMP THREADPRIVATE(cuvsurcv)
+  REAL,SAVE,ALLOCATABLE :: cvusurcu(:) ! cv/cu (u points)
+  !$OMP THREADPRIVATE(cvusurcu)
+  REAL,SAVE,ALLOCATABLE :: aire(:) ! cell area
+  !$OMP THREADPRIVATE(aire)
+  REAL,SAVE :: apoln ! area of north pole points
+  !$OMP THREADPRIVATE(apoln)
+  REAL,SAVE :: apols ! area of south pole points
+  !$OMP THREADPRIVATE(apols)
+  REAL,SAVE,ALLOCATABLE :: aireu(:) ! area of u cells 
+  !$OMP THREADPRIVATE(aireu)
+  REAL,SAVE,ALLOCATABLE :: airev(:) ! area of v cells
+  !$OMP THREADPRIVATE(airev)
+
+  ! Local variables for horiz mass flux in slab
+  LOGICAL,SAVE :: alpha_var
+  !$OMP THREADPRIVATE(alpha_var)
+  LOGICAL,SAVE :: slab_upstream
+  !$OMP THREADPRIVATE(slab_upstream)
+
+  REAL,SAVE,ALLOCATABLE :: zmasqu(:) ! continent mask for zonal mass flux
+  !$OMP THREADPRIVATE(zmasqu)
+  REAL,SAVE,ALLOCATABLE :: zmasqv(:) ! continent mask for merid mass flux
+  !$OMP THREADPRIVATE(zmasqv)
+  REAL,SAVE,ALLOCATABLE :: unsfv(:) ! 1/f, v points
+  !$OMP THREADPRIVATE(unsfv)
+  REAL,SAVE,ALLOCATABLE :: unsev(:) ! 1/epsilon (drag)
+  !$OMP THREADPRIVATE(unsev)
+  REAL,SAVE,ALLOCATABLE :: unsfu(:) ! 1/F, u points
+  !$OMP THREADPRIVATE(unsfu)
+  REAL,SAVE,ALLOCATABLE :: unseu(:)
+  !$OMP THREADPRIVATE(unseu)
+
+  ! Routines from dyn3d, valid on global dynamics grid only:
+  PRIVATE :: gr_fi_dyn, gr_dyn_fi ! to go between 1D nd 2D horiz grid
+  PRIVATE :: gr_scal_v,gr_v_scal,gr_scal_u ! change on 2D grid U,V, T points
+  PRIVATE :: grad,diverg
+
+CONTAINS
+ 
+  SUBROUTINE ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez_,fext_,unsaire_,&
+                                  cu_,cuvsurcv_,cv_,cvusurcu_, &
+                                  aire_,apoln_,apols_, &
+                                  aireu_,airev_)
+    ! number of points in lon, lat
+    IMPLICIT NONE
+    ! Routine copies some geometry variables from the dynamical core
+    ! see global vars for meaning
+    INTEGER,INTENT(IN) :: ip1jm
+    INTEGER,INTENT(IN) :: ip1jmp1
+    REAL,INTENT(IN) :: unsairez_(ip1jm)
+    REAL,INTENT(IN) :: fext_(ip1jm)
+    REAL,INTENT(IN) :: unsaire_(ip1jmp1)
+    REAL,INTENT(IN) :: cu_(ip1jmp1)
+    REAL,INTENT(IN) :: cuvsurcv_(ip1jm)
+    REAL,INTENT(IN) :: cv_(ip1jm)
+    REAL,INTENT(IN) :: cvusurcu_(ip1jmp1)
+    REAL,INTENT(IN) :: aire_(ip1jmp1)
+    REAL,INTENT(IN) :: apoln_
+    REAL,INTENT(IN) :: apols_
+    REAL,INTENT(IN) :: aireu_(ip1jmp1)
+    REAL,INTENT(IN) :: airev_(ip1jm)
+
+
+    ! Sanity check on dimensions
+    if ((ip1jm.ne.((nbp_lon+1)*(nbp_lat-1))).or. &
+        (ip1jmp1.ne.((nbp_lon+1)*nbp_lat))) then
+      write(*,*) "ini_slab_transp_geom Error: wrong array sizes"
+      stop
+    endif
+! Allocations could be done only on master process/thread...
+    allocate(unsairez(ip1jm))
+    unsairez(:)=unsairez_(:)
+    allocate(fext(ip1jm))
+    fext(:)=fext_(:)
+    allocate(unsaire(ip1jmp1))
+    unsaire(:)=unsaire_(:)
+    allocate(cu(ip1jmp1))
+    cu(:)=cu_(:)
+    allocate(cuvsurcv(ip1jm))
+    cuvsurcv(:)=cuvsurcv_(:)
+    allocate(cv(ip1jm))
+    cv(:)=cv_(:)
+    allocate(cvusurcu(ip1jmp1))
+    cvusurcu(:)=cvusurcu_(:)
+    allocate(aire(ip1jmp1))
+    aire(:)=aire_(:)
+    apoln=apoln_
+    apols=apols_
+    allocate(aireu(ip1jmp1))
+    aireu(:)=aireu_(:)
+    allocate(airev(ip1jm))
+    airev(:)=airev_(:)
+
+  END SUBROUTINE ini_slab_transp_geom
+
+  SUBROUTINE ini_slab_transp(zmasq)
+
+!    USE ioipsl_getin_p_mod, only: getin_p
+    USE IOIPSL, ONLY : getin 
+    IMPLICIT NONE
+
+    REAL zmasq(klon_glo) ! ocean / continent mask, 1=continent
+    REAL zmasq_2d((nbp_lon+1)*nbp_lat)
+    REAL ff((nbp_lon+1)*(nbp_lat-1)) ! Coriolis parameter
+    REAL eps ! epsilon friction timescale (s-1)
+    INTEGER :: slab_ekman 
+    INTEGER i
+    INTEGER :: iim,iip1,jjp1,ip1jm,ip1jmp1
+
+! Some definition for grid size
+    ip1jm=(nbp_lon+1)*(nbp_lat-1)
+    ip1jmp1=(nbp_lon+1)*nbp_lat
+    iim=nbp_lon
+    iip1=nbp_lon+1
+    jjp1=nbp_lat
+    ip1jm=(nbp_lon+1)*(nbp_lat-1)
+    ip1jmp1=(nbp_lon+1)*nbp_lat
+
+! Define ocean / continent mask (no flux into continent cell)
+    allocate(zmasqu(ip1jmp1))
+    allocate(zmasqv(ip1jm))
+    zmasqu=1.
+    zmasqv=1.
+
+    ! convert mask to 2D grid
+    CALL gr_fi_dyn(1,iip1,jjp1,zmasq,zmasq_2d)
+    ! put flux mask to 0 at boundaries of continent cells
+    DO i=1,ip1jmp1-1
+      IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+1).GT.1e-5) THEN
+              zmasqu(i)=0.
+      ENDIF
+    END DO
+    DO i=iip1,ip1jmp1,iip1
+      zmasqu(i)=zmasqu(i-iim)
+    END DO
+    DO i=1,ip1jm
+      IF (zmasq_2d(i).GT.1e-5 .OR. zmasq_2d(i+iip1).GT.1e-5) THEN
+              zmasqv(i)=0.
+      END IF
+    END DO
+    slab_ekman=2
+    CALL getin("slab_ekman",slab_ekman)
+! Coriolis and friction for Ekman transport
+    IF (slab_ekman.GT.0) THEN
+      allocate(unsfv(ip1jm))
+      allocate(unsev(ip1jm))
+      allocate(unsfu(ip1jmp1))
+      allocate(unseu(ip1jmp1))
+
+      eps=1e-5 ! Drag
+      CALL getin('slab_eps',eps)
+      print *,'epsilon=',eps
+      ff=fext*unsairez ! Coriolis
+      ! coefs to convert tau_x, tau_y to Ekman mass fluxes
+      ! on 2D grid v points
+      DO i=1,ip1jm
+        unsev(i)=eps/(ff(i)*ff(i)+eps**2)
+        unsfv(i)=ff(i)/(ff(i)*ff(i)+eps**2)
+      ENDDO
+      ! compute values on 2D u grid
+      CALL gr_v_scal(1,unsfv,unsfu)
+      CALL gr_v_scal(1,unsev,unseu)
+    END IF
+
+! Other options for Ekman transport
+    ! Alpha variable?
+      alpha_var=.FALSE.
+      CALL getin('slab_alphav',alpha_var)
+      print *,'alpha variable',alpha_var
+!  centered ou upstream scheme for meridional transport
+      slab_upstream=.FALSE.
+      CALL getin('slab_upstream',slab_upstream)
+      print *,'upstream slab scheme',slab_upstream
+  
+  END SUBROUTINE ini_slab_transp
+
+  SUBROUTINE divgrad_phy(nlevs,temp,delta)
+! Computes temperature tendency due to horizontal diffusion :
+! T Laplacian, later multiplied by diffusion coef and time-step
+
+    IMPLICIT NONE
+ 
+    INTEGER, INTENT(IN) :: nlevs ! nlevs : slab layers
+    REAL, INTENT(IN)   :: temp(klon_glo,nlevs) ! slab temperature
+    REAL , INTENT(OUT) :: delta(klon_glo,nlevs) ! temp laplacian (heat flux div.)
+    REAL :: delta_2d((nbp_lon+1)*nbp_lat,nlevs)
+    REAL ghx((nbp_lon+1)*nbp_lat,nlevs), ghy((nbp_lon+1)*(nbp_lat-1),nlevs)
+    INTEGER :: ll,iip1,jjp1
+ 
+    iip1=nbp_lon+1
+    jjp1=nbp_lat
+    
+    ! transpose temp to 2D horiz. grid
+    CALL gr_fi_dyn(nlevs,iip1,jjp1,temp,delta_2d)
+    ! computes gradient (proportional to heat flx)
+    CALL grad(nlevs,delta_2d,ghx,ghy)
+    ! put flux to 0 at ocean / continent boundary
+    DO ll=1,nlevs
+        ghx(:,ll)=ghx(:,ll)*zmasqu
+        ghy(:,ll)=ghy(:,ll)*zmasqv
+    END DO
+    ! flux divergence
+    CALL diverg(nlevs,ghx,ghy,delta_2d)
+    ! laplacian back to 1D grid
+    CALL gr_dyn_fi(nlevs,iip1,jjp1,delta_2d,delta)
+ 
+    RETURN
+  END SUBROUTINE divgrad_phy
+
+  SUBROUTINE slab_ekman1(tx_phy,ty_phy,ts_phy,dt_phy)
+! 1.5 Layer Ekman transport temperature tendency
+! note : tendency dt later multiplied by (delta t)/(rho.H)
+! to convert from divergence of heat fluxes to T
+
+      IMPLICIT NONE
+      
+      ! tx, ty : wind stress (different grids)
+      ! fluxm, fluz : mass *or heat* fluxes
+      ! dt : temperature tendency
+      INTEGER ij
+
+      ! ts surface temp, td deep temp (diagnosed)
+      REAL ts_phy(klon_glo)
+      REAL ts((nbp_lon+1)*nbp_lat),td((nbp_lon+1)*nbp_lat)
+      ! zonal and meridional wind stress 
+      REAL tx_phy(klon_glo),ty_phy(klon_glo)
+      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
+      REAL txv((nbp_lon+1)*(nbp_lat-1)),tyv((nbp_lon+1)*(nbp_lat-1))
+      ! zonal and meridional Ekman mass fluxes at u, v points (2D grid)
+      REAL fluxz((nbp_lon+1)*nbp_lat),fluxm((nbp_lon+1)*(nbp_lat-1))
+      ! vertical  and absolute mass fluxes (to estimate alpha)
+      REAL fluxv((nbp_lon+1)*nbp_lat),fluxt((nbp_lon+1)*(nbp_lat-1))
+      ! temperature tendency
+      REAL dt((nbp_lon+1)*nbp_lat),dt_phy(klon_glo)
+      REAL alpha((nbp_lon+1)*nbp_lat) ! deep temperature coef
+
+      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
+
+! Grid definitions
+      iim=nbp_lon
+      iip1=nbp_lon+1
+      iip2=nbp_lon+2
+      jjp1=nbp_lat
+      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
+      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+
+! Convert taux,y to 2D  scalar grid
+! Note: 2D grid size = iim*jjm. iip1=iim+1
+! First and last points in zonal direction are the same 
+! we use 1 index ij from 1 to (iim+1)*(jjm+1)
+      CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu)
+      CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu)
+! convert to u,v grid (Arakawa C)
+! Multiply by f or eps to get mass flux
+      ! Meridional mass flux
+      CALL gr_scal_v(1,txu,txv) ! wind stress at v points
+      CALL gr_scal_v(1,tyu,tyv)
+      fluxm=tyv*unsev-txv*unsfv ! in kg.s-1.m-1 (zonal distance)
+      ! Zonal mass flux
+      CALL gr_scal_u(1,txu,txu) ! wind stress at u points
+      CALL gr_scal_u(1,tyu,tyu)
+      fluxz=tyu*unsfu+txu*unseu
+            
+! Correct flux: continent mask and horiz grid size
+      ! multiply m-flux by mask and dx: flux in kg.s-1
+      fluxm=fluxm*cv*cuvsurcv*zmasqv 
+      ! multiply z-flux by mask and dy: flux in kg.s-1
+      fluxz=fluxz*cu*cvusurcu*zmasqu
+
+! Compute vertical  and absolute mass flux (for variable alpha)
+      IF (alpha_var) THEN
+        DO ij=iip2,ip1jm
+        fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
+        fluxt(ij)=ABS(fluxz(ij))+ABS(fluxz(ij-1)) &
+               +ABS(fluxm(ij))+ABS(fluxm(ij-iip1))
+        ENDDO
+        DO ij=iip1,ip1jmi1,iip1
+            fluxt(ij+1)=fluxt(ij+iip1)
+            fluxv(ij+1)=fluxv(ij+iip1)
+        END DO
+        fluxt(1)=SUM(ABS(fluxm(1:iim)))
+        fluxt(ip1jmp1)=SUM(ABS(fluxm(ip1jm-iim:ip1jm-1)))
+        fluxv(1)=-SUM(fluxm(1:iim))
+        fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
+        fluxt=MAX(fluxt,1.e-10)
+      ENDIF
+
+! Compute alpha coefficient. 
+! Tdeep = Tsurf * alpha + 271.35 * (1-alpha)
+      IF (alpha_var) THEN
+          ! increase alpha (and Tdeep) in downwelling regions
+          ! and decrease in upwelling regions
+          ! to avoid "hot spots" where there is surface convergence
+          DO ij=iip2,ip1jm
+              alpha(ij)=2./3.-fluxv(ij)/fluxt(ij)/3.
+          ENDDO
+          alpha(1:iip1)=2./3.-fluxv(1)/fluxt(1)/3.
+          alpha(ip1jm+1:ip1jmp1)=2./3.-fluxv(ip1jmp1)/fluxt(ip1jmp1)/3.
+      ELSE
+          alpha(:)=2./3.
+          ! Tsurf-Tdeep ~ 10° in the Tropics
+      ENDIF
+
+! Estimate deep temperature
+      CALL gr_fi_dyn(1,iip1,jjp1,ts_phy,ts)
+      DO ij=1,ip1jmp1
+         td(ij)=271.35+(ts(ij)-271.35)*alpha(ij)
+         td(ij)=MIN(td(ij),ts(ij))
+      END DO
+       
+! Meridional heat flux: multiply mass flux by (ts-td)
+! flux in K.kg.s-1
+      IF (slab_upstream) THEN
+        ! upstream scheme to avoid hot spots
+        DO ij=1,ip1jm
+        IF (fluxm(ij).GE.0.) THEN
+           fluxm(ij)=fluxm(ij)*(ts(ij+iip1)-td(ij))
+        ELSE
+           fluxm(ij)=fluxm(ij)*(ts(ij)-td(ij+iip1))
+        END IF
+        END DO
+      ELSE
+        ! centered scheme better in mid-latitudes
+        DO ij=1,ip1jm
+          fluxm(ij)=fluxm(ij)*(ts(ij+iip1)+ts(ij)-td(ij)-td(ij+iip1))/2.
+        END DO
+      ENDIF
+
+! Zonal heat flux
+! upstream scheme
+      DO ij=iip2,ip1jm
+          fluxz(ij)=fluxz(ij)*(ts(ij)+ts(ij+1)-td(ij+1)-td(ij))/2.
+      END DO
+      DO ij=iip1*2,ip1jmp1,iip1
+           fluxz(ij)=fluxz(ij-iim)
+      END DO
+
+! temperature tendency = divergence of heat fluxes
+! dt in K.s-1.kg.m-2 (T trend times mass/horiz surface)
+      DO ij=iip2,ip1jm
+         dt(ij)=(fluxz(ij-1)-fluxz(ij)+fluxm(ij)-fluxm(ij-iip1)) &
+               /aire(ij) ! aire : grid area 
+      END DO
+      DO ij=iip1,ip1jmi1,iip1
+         dt(ij+1)=dt(ij+iip1) 
+      END DO
+! special treatment at the Poles
+      dt(1)=SUM(fluxm(1:iim))/apoln
+      dt(ip1jmp1)=-SUM(fluxm(ip1jm-iim:ip1jm-1))/apols
+      dt(2:iip1)=dt(1)
+      dt(ip1jm+1:ip1jmp1)=dt(ip1jmp1)
+
+! tendencies back to 1D grid
+      CALL gr_dyn_fi(1,iip1,jjp1,dt,dt_phy)
+
+      RETURN
+  END SUBROUTINE slab_ekman1
+
+  SUBROUTINE slab_ekman2(tx_phy,ty_phy,ts_phy,dt_phy)
+! Temperature tendency for 2-layers slab ocean
+! note : tendency dt later multiplied by (delta time)/(rho.H)
+! to convert from divergence of heat fluxes to T
+      
+      IMPLICIT NONE
+      
+      ! Here, temperature and flux variables are on 2 layers
+      INTEGER ij
+
+      REAL tx_phy(klon_glo),ty_phy(klon_glo)
+      REAL txv((nbp_lon+1)*(nbp_lat-1)), tyv((nbp_lon+1)*(nbp_lat-1))
+      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
+      ! slab temperature on  1D, 2D grid
+      REAL ts_phy(klon_glo,2), ts((nbp_lon+1)*nbp_lat,2)
+      ! zonal and meridional Ekman mass flux at u, v points (2D grid)
+      REAL fluxz((nbp_lon+1)*nbp_lat), fluxm((nbp_lon+1)*(nbp_lat-1))
+      ! vertical mass flux between the 2 layers
+      REAL fluxv((nbp_lon+1)*nbp_lat)
+      ! zonal and meridional heat fluxes
+      REAL fluxtz((nbp_lon+1)*nbp_lat,2)
+      REAL fluxtm((nbp_lon+1)*(nbp_lat-1),2)
+      ! temperature tendency (in K.s-1.kg.m-2)
+      REAL dt((nbp_lon+1)*nbp_lat,2), dt_phy(klon_glo,2)
+
+      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
+
+! Grid definitions
+      iim=nbp_lon
+      iip1=nbp_lon+1
+      iip2=nbp_lon+2
+      jjp1=nbp_lat
+      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
+      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+
+! Convert taux,y to 2D  scalar grid
+      CALL gr_fi_dyn(1,iip1,jjp1,tx_phy,txu)
+      CALL gr_fi_dyn(1,iip1,jjp1,ty_phy,tyu)
+! Multiply taux,y by f or eps, and convert to 2D u,v grids
+! (Arakawa C grid)
+      ! Meridional flux
+      CALL gr_scal_v(1,txu,txv) ! wind stress at v points
+      CALL gr_scal_v(1,tyu,tyv)
+      fluxm=tyv*unsev-txv*unsfv
+      ! Zonal flux
+      CALL gr_scal_u(1,txu,txu) ! wind stress at u points
+      CALL gr_scal_u(1,tyu,tyu)
+      fluxz=tyu*unsfu+txu*unseu
+            
+! Convert temperature to 2D grid 
+      CALL gr_fi_dyn(2,iip1,jjp1,ts_phy,ts)
+       
+!  Mass fluxes (apply continent mask, multiply by horiz grid dimension) 
+      fluxm=fluxm*cv*cuvsurcv*zmasqv
+      fluxz=fluxz*cu*cvusurcu*zmasqu
+!  Vertical mass flux from mass budget (divergence of horiz fluxes)
+      DO ij=iip2,ip1jm
+        fluxv(ij)=fluxz(ij)-fluxz(ij-1)-fluxm(ij)+fluxm(ij-iip1)
+      ENDDO
+      DO ij=iip1,ip1jmi1,iip1
+         fluxv(ij+1)=fluxv(ij+iip1)
+      END DO
+!  vertical mass flux at Poles
+      fluxv(1)=-SUM(fluxm(1:iim))     
+      fluxv(ip1jmp1)=SUM(fluxm(ip1jm-iim:ip1jm-1))
+      fluxv=fluxv
+
+! Meridional heat fluxes 
+      DO ij=1,ip1jm
+          ! centered scheme
+          fluxtm(ij,1)=fluxm(ij)*(ts(ij+iip1,1)+ts(ij,1))/2.
+          fluxtm(ij,2)=-fluxm(ij)*(ts(ij+iip1,2)+ts(ij,2))/2.
+      END DO
+
+! Zonal heat fluxes
+! Schema upstream      
+      DO ij=iip2,ip1jm
+      IF (fluxz(ij).GE.0.) THEN
+             fluxtz(ij,1)=fluxz(ij)*ts(ij,1)
+             fluxtz(ij,2)=-fluxz(ij)*ts(ij+1,2)
+      ELSE
+             fluxtz(ij,1)=fluxz(ij)*ts(ij+1,1)
+             fluxtz(ij,2)=-fluxz(ij)*ts(ij,2)
+      ENDIF
+      ENDDO
+      DO ij=iip1*2,ip1jmp1,iip1
+             fluxtz(ij,:)=fluxtz(ij-iim,:)
+      END DO
+                   
+! Temperature tendency :
+      DO ij=iip2,ip1jm
+! divergence of horizontal heat fluxes
+         dt(ij,:)=fluxtz(ij-1,:)-fluxtz(ij,:) &
+                 +fluxtm(ij,:)-fluxtm(ij-iip1,:)
+! + vertical heat flux (mass flux * T, upstream scheme)
+         IF (fluxv(ij).GT.0.) THEN
+           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,2)
+           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,2)
+         ELSE
+           dt(ij,1)=dt(ij,1)+fluxv(ij)*ts(ij,1)
+           dt(ij,2)=dt(ij,2)-fluxv(ij)*ts(ij,1)
+         ENDIF
+         ! divide by cell area
+         dt(ij,:)=dt(ij,:)/aire(ij)
+      END DO
+      DO ij=iip1,ip1jmi1,iip1
+         dt(ij+1,:)=dt(ij+iip1,:) 
+      END DO
+! Pôles
+      dt(1,:)=SUM(fluxtm(1:iim,:),dim=1)
+        IF (fluxv(1).GT.0.) THEN
+          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,2)
+          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,2)
+        ELSE
+          dt(1,1)=dt(1,1)+fluxv(1)*ts(1,1)
+          dt(1,2)=dt(1,2)-fluxv(1)*ts(1,1)
+        ENDIF
+      dt(1,:)=dt(1,:)/apoln
+      dt(ip1jmp1,:)=-SUM(fluxtm(ip1jm-iim:ip1jm-1,:),dim=1)
+       IF (fluxv(ip1jmp1).GT.0.) THEN
+         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,2)
+         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,2)
+       ELSE
+         dt(ip1jmp1,1)=dt(ip1jmp1,1)+fluxv(ip1jmp1)*ts(ip1jmp1,1)
+         dt(ip1jmp1,2)=dt(ip1jmp1,2)-fluxv(ip1jmp1)*ts(ip1jmp1,1)
+       ENDIF
+      dt(ip1jmp1,:)=dt(ip1jmp1,:)/apols
+      dt(2:iip1,1)=dt(1,1)
+      dt(2:iip1,2)=dt(1,2)
+      dt(ip1jm+1:ip1jmp1,1)=dt(ip1jmp1,1)
+      dt(ip1jm+1:ip1jmp1,2)=dt(ip1jmp1,2)
+
+! T tendency back to 1D grid...
+      CALL gr_dyn_fi(2,iip1,jjp1,dt,dt_phy)
+
+      RETURN
+  END SUBROUTINE slab_ekman2
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  SUBROUTINE gr_fi_dyn(nfield,im,jm,pfi,pdyn)
+  ! Transfer a variable from 1D "physics" grid to 2D "dynamics" grid
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: im,jm,nfield
+  REAL,INTENT(IN) :: pfi(klon_glo,nfield) ! on 1D grid
+  REAL,INTENT(OUT) :: pdyn(im,jm,nfield) ! on 2D grid
+
+  INTEGER :: i,j,ifield,ig
+
+  DO ifield=1,nfield
+    ! Handle poles
+    DO i=1,im
+      pdyn(i,1,ifield)=pfi(1,ifield)
+      pdyn(i,jm,ifield)=pfi(klon_glo,ifield)
+    ENDDO
+    ! Other points
+    DO j=2,jm-1
+      ig=2+(j-2)*(im-1)
+      CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
+      pdyn(im,j,ifield)=pdyn(1,j,ifield)
+    ENDDO
+  ENDDO ! of DO ifield=1,nfield
+
+  END SUBROUTINE gr_fi_dyn
+
+  SUBROUTINE gr_dyn_fi(nfield,im,jm,pdyn,pfi)
+  ! Transfer a variable from 2D "dynamics" grid to 1D "physics" grid
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: im,jm,nfield
+  REAL,INTENT(IN) :: pdyn(im,jm,nfield) ! on 2D grid
+  REAL,INTENT(OUT) :: pfi(klon_glo,nfield) ! on 1D grid
+
+  INTEGER j,ifield,ig
+
+  ! Sanity check:
+  IF(klon_glo.NE.2+(jm-2)*(im-1)) THEN
+    WRITE(*,*) "gr_dyn_fi error, wrong sizes"
+    STOP
+  ENDIF
+
+  ! Handle poles
+  CALL SCOPY(nfield,pdyn,im*jm,pfi,klon_glo)
+  CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(klon_glo,1),klon_glo)
+  ! Other points
+  DO ifield=1,nfield
+    DO j=2,jm-1
+      ig=2+(j-2)*(im-1)
+      CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE gr_dyn_fi
+
+  SUBROUTINE  grad(klevel,pg,pgx,pgy)
+  ! compute the covariant components pgx,pgy of the gradient of pg
+  ! pgx = d(pg)/dx * delta(x) = delta(pg)
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: klevel
+  REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel)
+
+  INTEGER :: l,ij
+  INTEGER :: iim,iip1,ip1jm,ip1jmp1
+
+  iim=nbp_lon
+  iip1=nbp_lon+1
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+
+  DO l=1,klevel
+    DO ij=1,ip1jmp1-1
+      pgx(ij,l)=pg(ij+1,l)-pg(ij,l)
+    ENDDO
+    ! correction for pgx(ip1,j,l) ...
+    ! ... pgx(iip1,j,l)=pgx(1,j,l) ...
+    DO ij=iip1,ip1jmp1,iip1
+      pgx(ij,l)=pgx(ij-iim,l)
+    ENDDO
+    DO ij=1,ip1jm
+      pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l)
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE grad
+
+  SUBROUTINE diverg(klevel,x,y,div)
+  ! computes the divergence of a vector field of components
+  ! x,y. x and y being covariant components
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: klevel
+  REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel)
+  REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel)
+
+  INTEGER :: l,ij
+  INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1
+
+  REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1)
+  REAL :: sumypn,sumyps
+  REAL,EXTERNAL :: SSUM
+
+  iim=nbp_lon
+  iip1=nbp_lon+1
+  iip2=nbp_lon+2
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+  ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
+
+  DO l=1,klevel
+    DO ij=iip2,ip1jm-1
+      div(ij+1,l)= &
+        cvusurcu(ij+1)*x(ij+1,l)-cvusurcu(ij)*x(ij,l)+ &
+        cuvsurcv(ij-iim)*y(ij-iim,l)-cuvsurcv(ij+1)*y(ij+1,l)
+    ENDDO
+    ! correction for div(1,j,l) ...
+    ! ... div(1,j,l)= div(iip1,j,l) ...
+    DO ij=iip2,ip1jm,iip1
+      div(ij,l)=div(ij+iim,l)
+    ENDDO
+    ! at the poles
+    DO ij=1,iim
+      aiy1(ij)=cuvsurcv(ij)*y(ij,l)
+      aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l)
+    ENDDO
+    sumypn=SSUM(iim,aiy1,1)/apoln
+    sumyps=SSUM(iim,aiy2,1)/apols
+    DO ij=1,iip1
+      div(ij,l)=-sumypn
+      div(ij+ip1jm,l)=sumyps
+    ENDDO
+    ! End (poles)
+  ENDDO ! of DO l=1,klevel
+
+  !!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 )
+  DO l=1,klevel
+    DO ij=iip2,ip1jm
+      div(ij,l)=div(ij,l)*unsaire(ij)
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE diverg
+
+  SUBROUTINE gr_v_scal(nx,x_v,x_scal)
+  ! convert values from v points to scalar points on C-grid
+  ! used to  compute unsfu, unseu (u points, but depends only on latitude)
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: nx ! number of levels or fields
+  REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
+  REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx)
+
+  INTEGER :: l,ij
+  INTEGER :: iip1,iip2,ip1jm,ip1jmp1
+
+  iip1=nbp_lon+1
+  iip2=nbp_lon+2
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+
+  DO l=1,nx
+    DO ij=iip2,ip1jm
+      x_scal(ij,l)= &
+                   (airev(ij-iip1)*x_v(ij-iip1,l)+airev(ij)*x_v(ij,l)) &
+                  /(airev(ij-iip1)+airev(ij))
+    ENDDO
+    DO ij=1,iip1
+      x_scal(ij,l)=0.
+    ENDDO
+    DO ij=ip1jm+1,ip1jmp1
+      x_scal(ij,l)=0.
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE gr_v_scal
+ 
+  SUBROUTINE gr_scal_v(nx,x_scal,x_v)
+  ! convert values from scalar points to v points on C-grid
+  ! used to compute wind stress at V points
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: nx ! number of levels or fields
+  REAL,INTENT(OUT) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
+  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
+
+  INTEGER :: l,ij
+  INTEGER :: iip1,ip1jm
+
+  iip1=nbp_lon+1
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+
+      DO l=1,nx
+        DO ij=1,ip1jm
+          x_v(ij,l)= &
+            (cu(ij)*cvusurcu(ij)*x_scal(ij,l)+ &
+            cu(ij+iip1)*cvusurcu(ij+iip1)*x_scal(ij+iip1,l)) &
+            /(cu(ij)*cvusurcu(ij)+cu(ij+iip1)*cvusurcu(ij+iip1))
+        ENDDO
+      ENDDO
+
+  END SUBROUTINE gr_scal_v
+
+  SUBROUTINE gr_scal_u(nx,x_scal,x_u)
+  ! convert values from scalar points to U points on C-grid
+  ! used to compute wind stress at U points
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: nx
+  REAL,INTENT(OUT) :: x_u((nbp_lon+1)*nbp_lat,nx)
+  REAL,INTENT(IN) :: x_scal((nbp_lon+1)*nbp_lat,nx)
+
+  INTEGER :: l,ij
+  INTEGER :: iip1,jjp1,ip1jmp1
+
+  iip1=nbp_lon+1
+  jjp1=nbp_lat
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+
+  DO l=1,nx
+     DO ij=1,ip1jmp1-1
+        x_u(ij,l)= &
+         (aire(ij)*x_scal(ij,l)+aire(ij+1)*x_scal(ij+1,l)) &
+         /(aire(ij)+aire(ij+1))
+     ENDDO
+  ENDDO
+
+  CALL SCOPY(nx*jjp1,x_u(1,1),iip1,x_u(iip1,1),iip1)
+
+  END SUBROUTINE gr_scal_u
+
+END MODULE slab_heat_transp_mod
