Index: /LMDZ4/trunk/libf/phylmd/ocean_forced_mod.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 781)
+++ /LMDZ4/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 781)
@@ -0,0 +1,365 @@
+!
+! $Header$
+!
+MODULE ocean_forced_mod
+!
+! This module is used for both the sub-surfaces ocean and sea-ice for the case of a 
+! forced ocean,  "ocean=force".
+!
+  USE surface_data,     ONLY : calice, calsno, tau_gl
+  USE fonte_neige_mod,  ONLY : fonte_neige
+  USE calcul_fluxs_mod, ONLY : calcul_fluxs
+  USE dimphy
+
+  IMPLICIT NONE
+
+  REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE   :: tmp_flux_o, tmp_flux_g
+  !$OMP THREADPRIVATE(tmp_flux_o,tmp_flux_g)
+
+CONTAINS
+!
+!****************************************************************************************
+! 
+  SUBROUTINE ocean_forced_init
+! Allocate fields needed for this module
+!    
+    INTEGER              :: error
+    CHARACTER (len = 80) :: abort_message
+    CHARACTER (len = 20) :: modname = 'ocean_forced_init'
+!****************************************************************************************
+
+    ALLOCATE(tmp_flux_o(1:klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation tmp_flux_o'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+
+    ALLOCATE(tmp_flux_g(1:klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation tmp_flux_g'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF    
+
+
+  END SUBROUTINE ocean_forced_init
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_forced_noice(itime, dtime, jour, knon, knindex, &
+       debut, &
+       p1lay, tq_cdrag, precip_rain, precip_snow, &
+       temp_air, spechum, &
+       petAcoef, peqAcoef, petBcoef, peqBcoef, &
+       ps, u1_lay, v1_lay, &
+       radsol, snow, qsurf, &
+       agesno, &
+       evap, fluxsens, fluxlat, &
+       tsurf_new, dflux_s, dflux_l, pctsrf_oce)
+!
+! This subroutine treats the "open ocean", all grid points that are not entierly covered
+! by ice.
+! The routine reads data from climatologie file and does some calculations at the 
+! surface. 
+!
+    INCLUDE "indicesol.inc"
+    INCLUDE "YOMCST.inc"
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                      :: itime, jour, knon
+    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
+    LOGICAL, INTENT(IN)                      :: debut
+    REAL, INTENT(IN)                         :: dtime
+    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
+    REAL, DIMENSION(klon), INTENT(IN)        :: tq_cdrag
+    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)        :: petAcoef, peqAcoef
+    REAL, DIMENSION(klon), INTENT(IN)        :: petBcoef, peqBcoef
+    REAL, DIMENSION(klon), INTENT(IN)        :: ps
+    REAL, DIMENSION(klon), INTENT(IN)        :: u1_lay, v1_lay
+
+! In/Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT)     :: radsol
+    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
+    REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno !? put to 0 in ocean
+  
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
+    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
+    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_oce
+
+! Local variables
+!****************************************************************************************
+    INTEGER                     :: i
+    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd
+    REAL, DIMENSION(klon)       :: alb_neig, tsurf_lim, zx_sl
+    LOGICAL                     :: check=.FALSE.
+    REAL, DIMENSION(klon,nbsrf) :: pctsrf_lim
+
+!****************************************************************************************
+! Start calculation
+!****************************************************************************************
+    IF (check) WRITE(*,*)' Entering ocean_forced_noice'
+
+!****************************************************************************************
+! 1)    
+! Read from climatologie file SST and fraction of sub-surfaces
+!
+!****************************************************************************************
+! Get from file tsurf_lim and pctsrf_lim
+    CALL interfoce_lim(itime, dtime, jour, & 
+         knon, knindex, &
+         debut, &
+         tsurf_lim, pctsrf_lim)
+    
+!****************************************************************************************
+! 2)
+! Flux calculation
+!
+!****************************************************************************************
+! Set some variables for calcul_fluxs
+    cal = 0.
+    beta = 1.
+    dif_grnd = 0.
+    alb_neig(:) = 0.
+    agesno(:) = 0.
+  
+! 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, tq_cdrag, ps, &
+         precip_rain, precip_snow, snow, qsurf,  &
+         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+         petAcoef, peqAcoef, petBcoef, peqBcoef, &
+         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
+!****************************************************************************************
+! 3)
+! Calculate tmp_flux_o
+!
+!****************************************************************************************
+!IM: flux ocean-atmosphere utile pour le "slab" ocean
+! The flux are written to hist file
+    tmp_flux_o(:) = 0.0
+    DO i=1, knon
+       zx_sl(i) = RLVTT
+       IF (tsurf_new(i) .LT. RTT) zx_sl(i) = RLSTT
+
+!IM     flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i)
+!       flux_o(i) = fluxsens(i) + fluxlat(i)
+       IF (pctsrf_lim(knindex(i),is_oce) .GT. epsfra) THEN
+          tmp_flux_o(knindex(i)) = fluxsens(i) + fluxlat(i)
+       ENDIF
+    ENDDO
+
+!****************************************************************************************
+! 4)
+! Return the new values for the ocean fraction
+!
+!****************************************************************************************
+
+    pctsrf_oce(:) = pctsrf_lim(:,is_oce)
+
+  END SUBROUTINE ocean_forced_noice
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_forced_ice(itime, dtime, jour, knon, knindex, &
+       debut, &
+       tsurf_in, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
+       petAcoef, peqAcoef, petBcoef, peqBcoef, &
+       ps, u1_lay, v1_lay, &
+       radsol, snow, qsurf, qsol, agesno, &
+       tsoil, alblw, evap, fluxsens, fluxlat, &
+       tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_sic)
+!
+! This subroutine treats the ocean where there is ice. 
+! The routine reads data from climatologie file and does flux calculations at the 
+! surface.
+!    
+    INCLUDE "indicesol.inc"
+    INCLUDE "dimsoil.h"
+    INCLUDE "YOMCST.inc"
+    INCLUDE "clesphys.inc"
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                  :: itime, jour, knon
+    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
+    LOGICAL, INTENT(IN)                  :: debut
+    REAL, INTENT(IN)                     :: dtime
+    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
+    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
+    REAL, DIMENSION(klon), INTENT(IN)    :: tq_cdrag
+    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
+    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
+    REAL, DIMENSION(klon), INTENT(IN)    :: ps
+    REAL, DIMENSION(klon), INTENT(IN)    :: u1_lay, v1_lay
+
+
+! In/Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
+    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+
+
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
+    REAL, DIMENSION(klon), INTENT(OUT)            :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
+    REAL, DIMENSION(klon), INTENT(OUT)            :: pctsrf_sic
+
+
+! Local variables
+!****************************************************************************************
+    LOGICAL                     :: check=.FALSE.
+    INTEGER                     :: i
+    REAL                        :: zfra
+    REAL, PARAMETER             :: t_grnd=271.35
+    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd, capsol
+    REAL, DIMENSION(klon)       :: alb_neig, tsurf_tmp
+    REAL, DIMENSION(klon)       :: soilcap, soilflux
+    REAL, DIMENSION(klon,nbsrf) :: pctsrf_lim
+
+!****************************************************************************************
+! Start calculation
+!****************************************************************************************
+    IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon 
+
+!****************************************************************************************
+! 1)    
+! Read from climatologie file SST and fraction of sub-surfaces
+!
+!****************************************************************************************
+    CALL interfoce_lim(itime, dtime, jour, & 
+         knon, knindex, &
+         debut, &
+         tsurf_tmp, pctsrf_lim)
+  
+    DO i = 1, knon
+       tsurf_tmp(i) = tsurf_in(i)
+       IF (pctsrf_lim(knindex(i),is_sic) < EPSFRA) THEN
+          snow(i) = 0.0
+          tsurf_tmp(i) = RTT - 1.8
+          IF (soil_model) tsoil(i,:) = RTT -1.8
+       ENDIF
+    ENDDO
+
+!****************************************************************************************
+! 2) 
+! Flux calculation : tsurf_new, evap, fluxlat, fluxsens, 
+!                    dflux_s, dflux_l and qsurf
+!****************************************************************************************
+
+! calculate the parameters cal, beta, capsol and dif_grnd
+
+    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
+    
+    IF (soil_model) THEN 
+! update tsoil and calculate soilcap and soilflux
+       CALL soil(dtime, is_sic, knon,snow, tsurf_tmp, tsoil,soilcap, soilflux)
+       cal(1:knon) = RCPD / soilcap(1:knon)
+       radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
+       dif_grnd = 1.0 / tau_gl
+    ELSE 
+       dif_grnd = 1.0 / tau_gl
+       cal = RCPD * calice
+       WHERE (snow > 0.0) cal = RCPD * calsno 
+    ENDIF
+
+    beta = 1.0
+    CALL calcul_fluxs(knon, is_sic, dtime, &
+         tsurf_tmp, p1lay, cal, beta, tq_cdrag, ps, &
+         precip_rain, precip_snow, snow, qsurf,  &
+         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+         petAcoef, peqAcoef, petBcoef, peqBcoef, &
+         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
+!****************************************************************************************
+! 3)
+! Calculations due to snow and runoff
+!
+!****************************************************************************************
+    CALL fonte_neige( knon, is_sic, knindex, dtime, &
+         tsurf_tmp, precip_rain, precip_snow, &
+         snow, qsol, tsurf_new, evap)
+    
+! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)
+! 
+    CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:))  
+
+    WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.
+
+    alb_new(:) = 0.0
+    DO i=1, knon
+       zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0)))
+       alb_new(i) = alb_neig(i) * zfra +  0.6 * (1.0-zfra)
+    ENDDO
+!!      alb_new(1 : knon) = 0.6
+
+!****************************************************************************************
+! 4)
+! Calculate tmp_flux_g
+!
+!****************************************************************************************
+!
+    tmp_flux_g(:) = 0.0
+    DO i = 1, knon
+!IM: faire dependre le coefficient de conduction de la glace de mer
+!    de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
+!    actuel correspond a 3m de glace de mer, cf. L.Li
+!
+       IF (cal(i) .GT. 1.0e-15 .AND. pctsrf_lim(knindex(i),is_sic) .GT. epsfra) &
+            tmp_flux_g(knindex(i)) = (tsurf_new(i)-t_grnd) * dif_grnd(i) *RCPD/cal(i)
+       
+    ENDDO
+
+!!$    z0_new = 0.002
+!!$    z0_new = SQRT(z0_new**2+rugoro**2)
+    alblw(1:knon) = alb_new(1:knon)
+
+!****************************************************************************************
+! 5)
+! Return the new values for the seaice fraction
+!
+!****************************************************************************************
+
+    pctsrf_sic(:) = pctsrf_lim(:,is_sic)
+
+  END SUBROUTINE ocean_forced_ice
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_forced_get_vars(flux_o, flux_g)
+! Get some variables from module oceanforced.
+! This subroutine returns variables to a external routine
+
+    REAL, DIMENSION(klon), INTENT(OUT) :: flux_o
+    REAL, DIMENSION(klon), INTENT(OUT) :: flux_g
+
+! Initialize the output variables
+    flux_o(:) = tmp_flux_o(:)
+    flux_g(:) = tmp_flux_g(:)
+
+  END SUBROUTINE ocean_forced_get_vars
+!
+!****************************************************************************************
+!
+END MODULE ocean_forced_mod
+
+
+
+
+
+
Index: /LMDZ4/trunk/libf/phylmd/ocean_slab_mod.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 781)
+++ /LMDZ4/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 781)
@@ -0,0 +1,648 @@
+!
+! $Header$
+!
+MODULE ocean_slab_mod
+!
+! This module is used for both surface ocean and sea-ice when using the slab ocean,
+! "ocean=slab".
+!
+  USE surface_data,     ONLY : tau_gl, calice, calsno
+  USE fonte_neige_mod,  ONLY : fonte_neige
+  USE calcul_fluxs_mod, ONLY : calcul_fluxs
+  USE dimphy
+  
+  IMPLICIT NONE
+
+  INTEGER, PRIVATE, SAVE                           :: lmt_pas, julien, idayvrai
+  !$OMP THREADPRIVATE(lmt_pas,julien,idayvrai)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_seaice
+  !$OMP THREADPRIVATE(tmp_seaice)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_tslab_loc
+  !$OMP THREADPRIVATE(tmp_tslab_loc)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slab_bils
+  !$OMP THREADPRIVATE(slab_bils)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE , SAVE  :: lmt_bils
+  !$OMP THREADPRIVATE(lmt_bils)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: tmp_pctsrf_slab
+  !$OMP THREADPRIVATE(tmp_pctsrf_slab)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_tslab
+  !$OMP THREADPRIVATE(tmp_tslab)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_radsol
+  !$OMP THREADPRIVATE(tmp_radsol)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_flux_o, tmp_flux_g
+  !$OMP THREADPRIVATE(tmp_flux_o,tmp_flux_g)
+  LOGICAL, PRIVATE, SAVE                           :: check = .FALSE.
+  !$OMP THREADPRIVATE(check)
+
+CONTAINS
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_slab_init(dtime, tslab_rst, seaice_rst, pctsrf_rst)
+
+    INCLUDE "indicesol.inc"
+    INCLUDE "iniprint.h"
+
+! Input variables
+!****************************************************************************************
+    REAL, INTENT(IN)                         :: dtime
+! Variables read from restart file
+    REAL, DIMENSION(klon), INTENT(IN)        :: tslab_rst          
+    REAL, DIMENSION(klon), INTENT(IN)        :: seaice_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
+
+
+! Local variables
+!****************************************************************************************
+    INTEGER                :: error
+    CHARACTER (len = 80)   :: abort_message
+    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
+
+
+    WRITE(lunout,*) '************************'
+    WRITE(lunout,*) 'SLAB OCEAN est actif, prenez precautions !'
+    WRITE(lunout,*) '************************'    
+
+! Allocate variables initialize from restart fields
+    ALLOCATE(tmp_tslab(klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation tmp_tslab'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    tmp_tslab(:) = tslab_rst(:)
+
+    ALLOCATE(tmp_tslab_loc(klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation tmp_tslab_loc'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    tmp_tslab_loc(:) = tslab_rst(:)
+
+    ALLOCATE(tmp_seaice(klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation tmp_seaice'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    tmp_seaice(:) = seaice_rst(:)
+
+    ALLOCATE(tmp_pctsrf_slab(klon,nbsrf), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation tmp_pctsrf_slab'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    tmp_pctsrf_slab(:,:) = pctsrf_rst(:,:)
+    
+! Allocate some other variables internal in module mod_oceanslab
+    ALLOCATE(tmp_radsol(klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation tmp_radsol'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+
+    ALLOCATE(tmp_flux_o(klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation tmp_flux_o'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    
+    ALLOCATE(tmp_flux_g(klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation tmp_flux_g'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+
+! a mettre un slab_bils aussi en force !!!
+    ALLOCATE(slab_bils(klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation slab_bils'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    slab_bils(:) = 0.0   
+
+    ALLOCATE(lmt_bils(klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation lmt_bils'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+
+
+! pour une lecture une fois par jour    
+    lmt_pas = NINT(86400./dtime * 1.0) 
+
+  END SUBROUTINE ocean_slab_init
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_slab_noice( & 
+       dtime, knon, knindex, &
+       p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
+       petAcoef, peqAcoef, petBcoef, peqBcoef, &
+       ps, u1_lay, v1_lay, &
+       radsol, snow, qsurf, agesno, &
+       evap, fluxsens, fluxlat, &
+       tsurf_new, &
+       dflux_s, dflux_l, pctsrf_oce)
+
+    INCLUDE "indicesol.inc"
+    INCLUDE "iniprint.h"
+
+! Input arguments
+!****************************************************************************************
+    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)    :: tq_cdrag
+    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
+    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
+    REAL, DIMENSION(klon), INTENT(IN)    :: ps
+    REAL, DIMENSION(klon), INTENT(IN)    :: u1_lay, v1_lay
+
+! In/Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
+    REAL, DIMENSION(klon), INTENT(INOUT) :: snow
+    REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
+    
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
+    REAL, DIMENSION(klon), INTENT(OUT)   :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
+    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l      
+    REAL, DIMENSION(klon), INTENT(OUT)   :: pctsrf_oce
+
+! Local variables
+!****************************************************************************************
+    INTEGER                :: i
+    REAL, DIMENSION(klon)  :: cal, beta, dif_grnd
+    REAL, DIMENSION(klon)  :: alb_neig, tsurf_temp
+
+!****************************************************************************************
+    IF (check) WRITE(*,*)' Entering ocean_slab_noice'    
+
+    tsurf_new(1:knon) = tmp_tslab(knindex(1:knon))
+    pctsrf_oce(:)   = tmp_pctsrf_slab(:,is_oce)
+    
+    tsurf_temp(:) = tsurf_new(:)
+    cal = 0.
+    beta = 1.
+    dif_grnd = 0.
+    alb_neig(:) = 0.
+    agesno(:) = 0.
+    
+    CALL calcul_fluxs(knon, is_oce, dtime, &
+         tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
+         precip_rain, precip_snow, snow, qsurf,  &
+         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+         petAcoef, peqAcoef, petBcoef, peqBcoef, &
+         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+    
+    tmp_flux_o(:) = 0.0
+    tmp_radsol(:) = 0.0
+
+    DO i=1, knon
+       tmp_radsol(knindex(i))=radsol(i)
+       
+       IF (pctsrf_oce(knindex(i)) .GT. epsfra) &
+            tmp_flux_o(knindex(i)) = fluxsens(i) + fluxlat(i)
+    ENDDO
+    
+  END SUBROUTINE ocean_slab_noice
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_slab_ice(   &
+       itime, dtime, jour, knon, knindex, &
+       debut, &
+       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
+       petAcoef, peqAcoef, petBcoef, peqBcoef, &
+       ps, u1_lay, v1_lay, &
+       radsol, snow, qsurf, qsol, agesno, &
+       tsoil, &
+       alblw, evap, fluxsens, fluxlat, &
+       tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_sic)
+
+    INCLUDE "indicesol.inc"
+    INCLUDE "dimsoil.h"
+    INCLUDE "YOMCST.inc"
+    INCLUDE "iniprint.h"
+    INCLUDE "clesphys.inc"
+
+! Input arguments  
+!****************************************************************************************
+    INTEGER, INTENT(IN)                  :: itime, jour, knon
+    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
+    REAL, INTENT(IN)                     :: dtime
+    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf
+    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
+    REAL, DIMENSION(klon), INTENT(IN)    :: tq_cdrag
+    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
+    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
+    REAL, DIMENSION(klon), INTENT(IN)    :: ps
+    REAL, DIMENSION(klon), INTENT(IN)    :: u1_lay, v1_lay
+    LOGICAL, INTENT(IN)                  :: debut
+
+!In/Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
+    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
+    REAL, DIMENSION(klon), INTENT(OUT)            :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
+    REAL, DIMENSION(klon), INTENT(OUT)            :: pctsrf_sic
+
+! Local variables
+!****************************************************************************************
+    INTEGER                              :: i
+    REAL, DIMENSION(klon)                :: cal, beta, dif_grnd, capsol
+    REAL, DIMENSION(klon)                :: alb_neig, tsurf_temp
+    REAL, DIMENSION(klon)                :: soilcap, soilflux
+    REAL, DIMENSION(klon)                :: zfra
+    REAL, PARAMETER                      :: t_grnd=271.35
+    REAL                                 :: amn, amx 
+    REAL, DIMENSION(klon)                :: tslab
+    REAL, DIMENSION(klon)                :: seaice ! glace de mer (kg/m2)
+    REAL, DIMENSION(klon,nbsrf)          :: pctsrf_new
+
+!****************************************************************************************
+
+    IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon 
+
+! Initialization of output variables
+    alb_new(:) = 0.0
+
+!****************************************************************************************
+! 
+!
+!****************************************************************************************
+    IF ( ok_slab_sicOBS) THEN   
+       ! glace de mer observee, lecture conditions limites
+       CALL interfoce_lim(itime, dtime, jour, & 
+            knon, knindex, &
+            debut, &
+            tsurf_new, pctsrf_new)
+
+       tmp_pctsrf_slab(:,:) = pctsrf_new(:,:)
+       WRITE(lunout,*) 'jour lecture pctsrf_new sic =',jour
+
+    ELSE
+       pctsrf_new=tmp_pctsrf_slab
+    ENDIF
+
+    DO i = 1, knon
+       tsurf_new(i) = tsurf(i)
+       IF (pctsrf_new(knindex(i),is_sic) < EPSFRA) THEN
+          snow(i) = 0.0
+          tsurf_new(i) = RTT - 1.8
+          IF (soil_model) tsoil(i,:) = RTT -1.8
+       ENDIF
+    ENDDO
+    
+    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
+    
+    IF (soil_model) THEN 
+       CALL soil(dtime, is_sic, knon, snow, tsurf_new, tsoil, soilcap, soilflux)
+       cal(1:knon) = RCPD / soilcap(1:knon)
+       radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
+    ELSE 
+       dif_grnd = 1.0 / tau_gl
+       cal = RCPD * calice
+       WHERE (snow > 0.0) cal = RCPD * calsno 
+    ENDIF
+    tsurf_temp = tsurf_new
+    beta = 1.0
+
+!****************************************************************************************
+!
+!
+!****************************************************************************************
+    CALL calcul_fluxs(knon, is_sic, dtime, &
+         tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
+         precip_rain, precip_snow, snow, qsurf,  &
+         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+         petAcoef, peqAcoef, petBcoef, peqBcoef, &
+         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+    
+    CALL fonte_neige( knon, is_sic, knindex, dtime, &
+         tsurf_temp, precip_rain, precip_snow, &
+         snow, qsol, tsurf_new, evap)
+
+!****************************************************************************************
+!     calcul albedo
+!
+!****************************************************************************************
+    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
+    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
+    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
+    alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + & 
+         0.6 * (1.0-zfra(1:knon))
+    !!      alb_new(1 : knon) = 0.6
+
+!
+!IM: flux entre l'ocean et la glace de mer pour le "slab" ocean
+    tmp_flux_g(:) = 0.0
+    DO i = 1, knon
+!
+!IM: faire dependre le coefficient de conduction de la glace de mer
+!    de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
+!    actuel correspond a 3m de glace de mer, cf. L.Li
+!
+       IF ((cal(i).GT.1.0e-15) .AND. (pctsrf_new(knindex(i),is_sic) .GT. epsfra)) THEN
+          tmp_flux_g(knindex(i))=(tsurf_new(i)-t_grnd) &
+               * dif_grnd(i) *RCPD/cal(i)
+       ENDIF
+!
+!IM: Attention: ne pas initialiser le tmp_radsol puisque c'est deja fait sur is_oce;
+!IM:            tmp_radsol doit etre le flux solaire qui arrive sur l'ocean 
+!IM:            et non pas celui qui arrive sur la glace de mer
+    ENDDO
+    
+!
+! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
+!
+
+    IF (check) THEN
+       amn=MIN(tmp_tslab(1),1000.)
+       amx=MAX(tmp_tslab(1),-1000.)
+       DO i=2, klon
+          amn=MIN(tmp_tslab(i),amn)
+          amx=MAX(tmp_tslab(i),amx)
+       ENDDO
+       
+       WRITE(lunout,*) ' debut avant interfoce_slab min max tmp_tslab',amn,amx
+    ENDIF !(check) THEN
+
+    tslab = tmp_tslab   
+
+    CALL interfoce_slab(klon, debut, itime, dtime, jour, &
+         tmp_radsol, tmp_flux_o, tmp_flux_g, tmp_pctsrf_slab, &
+         tslab, seaice, pctsrf_new)
+    
+    tmp_pctsrf_slab=pctsrf_new
+
+    DO i=1, knon
+       tmp_tslab(knindex(i))=tslab(knindex(i))
+    ENDDO
+     
+!!$    z0_new = 0.002
+!!$    z0_new = SQRT(z0_new**2+rugoro**2)
+    alblw(1:knon) = alb_new(1:knon)
+        
+!****************************************************************************************
+! Return the fraction of sea-ice
+! NB! jg : Peut-etre un probleme, faut-il prend aussi tmp_pctsrf_slab(:,is_oce)???
+!****************************************************************************************
+    pctsrf_sic(:) =  tmp_pctsrf_slab(:,is_sic)
+
+
+  END SUBROUTINE ocean_slab_ice
+! 
+!****************************************************************************************
+!
+  SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour, &
+       radsol, fluxo, fluxg, pctsrf, &
+       tslab, seaice, pctsrf_slab)
+!
+! Cette routine calcule la temperature d'un slab ocean, la glace de mer 
+! et les pourcentages de la maille couverte par l'ocean libre et/ou 
+! la glace de mer pour un "slab" ocean de 50m
+!
+! Conception: Laurent Li
+! Re-ecriture + adaptation LMDZ4: I. Musat
+!
+! input:
+!   klon         nombre total de points de grille
+!   debut        logical: 1er appel a la physique
+!   itap         numero du pas de temps
+!   dtime        pas de temps de la physique (en s)
+!   ijour        jour dans l'annee en cours
+!   radsol       rayonnement net au sol (LW + SW)
+!   fluxo        flux turbulent (sensible + latent) sur les mailles oceaniques 
+!   fluxg        flux de conduction entre la surface de la glace de mer et l'ocean
+!   pctsrf       tableau des pourcentages de surface de chaque maille
+! output: 
+!   tslab        temperature de l'ocean libre
+!   seaice       glace de mer (kg/m2)
+!   pctsrf_slab  "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab
+
+    INCLUDE "indicesol.inc"
+    INCLUDE "YOMCST.inc"
+    INCLUDE "iniprint.h"
+    INCLUDE "clesphys.inc"
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                       :: klon
+    LOGICAL, INTENT(IN)                       :: debut    ! not used
+    INTEGER, INTENT(IN)                       :: itap
+    REAL, INTENT(IN)                          :: dtime       ! not used
+    INTEGER, INTENT(IN)                       :: ijour
+    REAL, DIMENSION(klon), INTENT(IN)         :: radsol
+    REAL, DIMENSION(klon), INTENT(IN)         :: fluxo
+    REAL, DIMENSION(klon), INTENT(IN)         :: fluxg
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN)  :: pctsrf
+
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)        :: tslab
+    REAL, DIMENSION(klon), INTENT(OUT)        :: seaice ! glace de mer (kg/m2)
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: pctsrf_slab
+
+! Local variables
+!****************************************************************************************
+    REAL                    :: amn, amx
+    REAL, PARAMETER         :: unjour=86400.
+    REAL, PARAMETER         :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)
+    REAL, PARAMETER         :: cbing=0.334e+05        ! J/kg
+    REAL, DIMENSION(klon)   :: siceh !hauteur de la glace de mer (m)
+    INTEGER                 :: i
+    REAL                    :: zz, za, zb
+!
+!****************************************************************************************
+!
+    julien = MOD(ijour,360)
+
+    IF (check ) THEN
+       amn=MIN(tmp_tslab_loc(1),1000.)
+       amx=MAX(tmp_tslab_loc(1),-1000.)
+       DO i=2, klon
+          amn=MIN(tmp_tslab_loc(i),amn)
+          amx=MAX(tmp_tslab_loc(i),amx)
+       ENDDO
+
+       WRITE(lunout,*) ' debut min max tslab',amn,amx 
+       WRITE(lunout,*) ' itap,lmt_pas unjour',itap,lmt_pas,unjour
+    ENDIF
+
+    pctsrf_slab(1:klon,1:nbsrf) = pctsrf(1:klon,1:nbsrf)
+!
+! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee
+! 
+    IF (MOD(itap,lmt_pas) .EQ. 1) THEN 
+       ! 1er pas de temps de la journee
+       idayvrai = ijour
+       CALL condsurf(julien,idayvrai, lmt_bils)
+    ENDIF
+
+    DO i = 1, klon
+       IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
+            (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
+!
+! fabriquer de la glace si congelation atteinte:
+!          
+          IF (tmp_tslab_loc(i).LT.(RTT-1.8)) THEN
+             zz =  (RTT-1.8)-tmp_tslab_loc(i)
+             tmp_seaice(i) = tmp_seaice(i) + cyang/cbing * zz
+             seaice(i) = tmp_seaice(i)
+             tmp_tslab_loc(i) = RTT-1.8
+          ENDIF
+!
+! faire fondre de la glace si temperature est superieure a 0:
+!
+          IF ((tmp_tslab_loc(i).GT.RTT) .AND. (tmp_seaice(i).GT.0.0)) THEN
+             zz = cyang/cbing * (tmp_tslab_loc(i)-RTT)
+             zz = MIN(zz,tmp_seaice(i))
+             tmp_seaice(i) = tmp_seaice(i) - zz
+             seaice(i) = tmp_seaice(i)
+             tmp_tslab_loc(i) = tmp_tslab_loc(i) - zz*cbing/cyang
+          ENDIF
+!
+! limiter la glace de mer a 10 metres (10000 kg/m2)
+!
+          IF(tmp_seaice(i).GT.45.) THEN
+             tmp_seaice(i) = MIN(tmp_seaice(i),10000.0)
+          ELSE
+             tmp_seaice(i) = 0. 
+          ENDIF
+          seaice(i) = tmp_seaice(i)
+          siceh(i)=tmp_seaice(i)/1000. !en metres
+!
+! determiner la nature du sol (glace de mer ou ocean libre):
+!
+! on fait dependre la fraction de seaice "pctsrf(i,is_sic)" 
+! de l'epaisseur de seaice :
+! pctsrf(i,is_sic)=1. si l'epaisseur de la glace de mer est >= a 20cm
+! et pctsrf(i,is_sic) croit lineairement avec seaice de 0. a 20cm d'epaisseur
+!
+
+          IF(.NOT.ok_slab_sicOBS) THEN
+             pctsrf_slab(i,is_sic)=MIN(siceh(i)/0.20, &
+                  1.-(pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)))
+             pctsrf_slab(i,is_oce)=1.0 - &
+                  (pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)+pctsrf_slab(i,is_sic))
+          ELSE
+             IF (i.EQ.1) WRITE(lunout,*) 'cas ok_slab_sicOBS TRUE : passe sur la modif.'
+          ENDIF !(.NOT.ok_slab_sicOBS) then
+       ENDIF !pctsrf
+    ENDDO
+!
+! Calculer le bilan du flux de chaleur au sol :
+!
+    DO i = 1, klon
+       za = radsol(i) + fluxo(i)
+       zb = fluxg(i)
+       IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
+            (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
+          slab_bils(i)=slab_bils(i)+(za*pctsrf_slab(i,is_oce) &
+               +zb*pctsrf_slab(i,is_sic))/ FLOAT(lmt_pas)
+       ENDIF
+    ENDDO !klon
+!
+! calcul tslab 
+!
+    IF (MOD(itap,lmt_pas).EQ.0) THEN !fin de journee
+!
+! calcul tslab 
+       amn=MIN(tmp_tslab_loc(1),1000.)
+       amx=MAX(tmp_tslab_loc(1),-1000.)
+       DO i = 1, klon
+          IF ((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
+               (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
+             tmp_tslab_loc(i) = tmp_tslab_loc(i) + &
+                  (slab_bils(i)-lmt_bils(i)) &
+                  /cyang*unjour
+
+! on remet l'accumulation a 0
+             slab_bils(i) = 0.
+          ENDIF !pctsrf
+!
+          IF (check) THEN
+             IF(i.EQ.1) THEN  
+                WRITE(lunout,*) 'i tmp_tslab_loc MOD(itap,lmt_pas).EQ.0 sicINTER',i,itap, &
+                     tmp_tslab_loc(i), tmp_seaice(i)
+             ENDIF
+             
+             amn=MIN(tmp_tslab_loc(i),amn)
+             amx=MAX(tmp_tslab_loc(i),amx)
+          ENDIF
+       ENDDO !klon
+    ENDIF !(MOD(itap,lmt_pas).EQ.0) THEN
+
+    IF ( check ) THEN
+       WRITE(lunout,*) 'fin min max tslab',amn,amx 
+    ENDIF
+
+    tslab  = tmp_tslab_loc
+    seaice = tmp_seaice
+
+  END SUBROUTINE interfoce_slab
+!
+!****************************************************************************************  
+!
+  SUBROUTINE ocean_slab_final(tslab_rst, seaice_rst)
+
+! This subroutine will send to phyredem the variables concerning the slab 
+! ocean that should be written to restart file.
+
+!****************************************************************************************
+
+    REAL, DIMENSION(klon), INTENT(OUT) :: tslab_rst
+    REAL, DIMENSION(klon), INTENT(OUT) :: seaice_rst
+
+!****************************************************************************************
+! Set the output variables
+    tslab_rst(:)  = tmp_tslab(:)
+!    tslab_rst(:)  = tmp_tslab_loc(:)
+    seaice_rst(:) = tmp_seaice(:)
+
+! Deallocation of all variables in module
+    DEALLOCATE(tmp_tslab, tmp_tslab_loc, tmp_pctsrf_slab)
+    DEALLOCATE(tmp_seaice, tmp_radsol, tmp_flux_o, tmp_flux_g)
+    DEALLOCATE(slab_bils, lmt_bils)
+
+  END SUBROUTINE ocean_slab_final
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_slab_get_vars(tslab_loc, seaice_loc, flux_o_loc, flux_g_loc)
+! "Get some variables from module ocean_slab_mod"
+! This subroutine prints variables to a external routine
+
+    REAL, DIMENSION(klon), INTENT(OUT) :: tslab_loc
+    REAL, DIMENSION(klon), INTENT(OUT) :: seaice_loc
+    REAL, DIMENSION(klon), INTENT(OUT) :: flux_o_loc
+    REAL, DIMENSION(klon), INTENT(OUT) :: flux_g_loc
+
+! Set the output variables
+    tslab_loc(:)  = tmp_tslab(:)
+    seaice_loc(:) = tmp_seaice(:)
+    flux_o_loc(:) = tmp_flux_o(:)
+    flux_g_loc(:) = tmp_flux_g(:)
+
+  END SUBROUTINE ocean_slab_get_vars
+!
+!****************************************************************************************
+!
+END MODULE ocean_slab_mod
Index: /LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 781)
+++ /LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 781)
@@ -0,0 +1,1279 @@
+!
+! $Header$
+!
+MODULE pbl_surface_mod
+!
+! Planetary Boundary Layer and Surface module
+!
+! This module manage the calculation of turbulent diffusion in the boundary layer 
+! and all interactions towards the differents sub-surfaces.
+!
+!
+  USE dimphy
+  USE mod_phys_lmdz_para,  ONLY : mpi_size
+  USE ioipsl
+  USE surface_data,        ONLY : ocean, ok_veget, debug_surf, newwind
+  USE surf_land_mod,       ONLY : surf_land
+  USE surf_landice_mod,    ONLY : surf_landice
+  USE surf_ocean_mod,      ONLY : surf_ocean
+  USE surf_seaice_mod,     ONLY : surf_seaice
+  USE cpl_mod,             ONLY : gath2cpl
+  USE climb_hq_mod,        ONLY : climb_hq_down, climb_hq_up
+  USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
+  USE coef_diff_turb_mod,  ONLY : coef_diff_turb
+
+!jg+ temporary
+  USE mod_clvent, ONLY : clvent, save_flux
+!jg-
+
+  IMPLICIT NONE
+
+! Declaration of variables saved in restart file
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: qsol
+  !$OMP THREADPRIVATE(qsol)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder
+  !$OMP THREADPRIVATE(fder)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: snow
+  !$OMP THREADPRIVATE(snow)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf
+  !$OMP THREADPRIVATE(qsurf)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: evap
+  !$OMP THREADPRIVATE(evap)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: rugos
+  !$OMP THREADPRIVATE(rugos)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: agesno
+  !$OMP THREADPRIVATE(agesno)
+  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: ftsoil  
+  !$OMP THREADPRIVATE(ftsoil)
+
+CONTAINS
+!
+!****************************************************************************************
+!
+  SUBROUTINE pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst,&
+       evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
+
+! This routine should be called after the restart file has been read.
+! This routine initialize the restart variables and does some validation tests
+! for the index of the different surfaces and tests the choice of type of ocean.
+
+    INCLUDE "indicesol.inc"
+    INCLUDE "dimsoil.h"
+    INCLUDE "iniprint.h"
+ 
+! Input variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(IN)                 :: qsol_rst
+    REAL, DIMENSION(klon), INTENT(IN)                 :: fder_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: snow_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: qsurf_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: evap_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: rugos_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: agesno_rst
+    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst
+
+  
+! Local variables
+!****************************************************************************************
+    INTEGER                       :: ierr
+    CHARACTER(len=80)             :: abort_message
+    CHARACTER(len = 20)           :: modname = 'pbl_surface_init'
+    
+
+!****************************************************************************************
+! Allocate and initialize module variables with fields read from restart file.
+!
+!****************************************************************************************    
+    ALLOCATE(qsol(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(fder(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(snow(klon,nbsrf), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(qsurf(klon,nbsrf), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(evap(klon,nbsrf), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(rugos(klon,nbsrf), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(agesno(klon,nbsrf), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
+
+
+    qsol(:)       = qsol_rst(:)
+    fder(:)       = fder_rst(:)
+    snow(:,:)     = snow_rst(:,:)
+    qsurf(:,:)    = qsurf_rst(:,:)
+    evap(:,:)     = evap_rst(:,:)
+    rugos(:,:)    = rugos_rst(:,:)
+    agesno(:,:)   = agesno_rst(:,:)
+    ftsoil(:,:,:) = ftsoil_rst(:,:,:)
+
+
+!****************************************************************************************
+! Test for sub-surface indices
+!
+!****************************************************************************************
+    IF (is_ter /= 1) THEN 
+      WRITE(lunout,*)" *** Warning ***"
+      WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter
+      WRITE(lunout,*)"or on doit commencer par les surfaces continentales"
+      abort_message="voir ci-dessus"
+      CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+
+    IF ( is_oce > is_sic ) THEN
+      WRITE(lunout,*)' *** Warning ***'
+      WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
+      WRITE(lunout,*)' l''ocean doit etre traite avant la banquise'
+      WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
+      abort_message='voir ci-dessus'
+      CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+
+    IF ( is_lic > is_sic ) THEN
+      WRITE(lunout,*)' *** Warning ***'
+      WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
+      WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer'
+      WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic
+      abort_message='voir ci-dessus'
+      CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+
+!****************************************************************************************
+! Validation of ocean mode
+!
+!****************************************************************************************
+
+    IF (ocean /= 'slab  ' .AND. ocean /= 'force ' .AND. ocean /= 'couple') THEN
+      WRITE(lunout,*)' *** Warning ***'
+      WRITE(lunout,*)'Option couplage pour l''ocean = ', ocean
+      abort_message='option pour l''ocean non valable'
+      CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+
+!****************************************************************************************
+! Test of coherence between variable ok_veget and cpp key CPP_VEGET
+!
+!****************************************************************************************
+    IF (ok_veget) THEN
+#ifndef CPP_VEGET
+       abort_message='Pb de coherence: ok_veget = .true. mais CPP_VEGET = .false.'
+       CALL abort_gcm(modname,abort_message,1)
+#endif
+    ENDIF
+
+
+  END SUBROUTINE pbl_surface_init
+!  
+!****************************************************************************************
+!  
+
+  SUBROUTINE pbl_surface( &
+       dtime,     date0,     itap,     jour,          &
+       debut,     lafin,                              &
+       rlon,      rlat,      rugoro,   rmu0,          &
+       rain_f,    snow_f,    solsw_m,  sollw_m,       &
+       t,         q,         u,        v,             &
+       pplay,     paprs,     pctsrf,                  &
+       ts,        albe,      alblw,    u10m,   v10m,  &
+       sollwdown, cdragh,    cdragm,   zu1,    zv1,   &
+       albsol,    albsollw,  zxsens,   zxevap,        &
+       zxtsol,    zxfluxlat, zt2m,     qsat2m,        &
+       d_t,       d_q,       d_u,      d_v,           & 
+       zcoefh,    pctsrf_new,                         &
+       qsol_d,    zq2m,      s_pblh,   s_plcl,        &
+       s_capCL,   s_oliqCL,  s_cteiCL, s_pblT,        &
+       s_therm,   s_trmb1,   s_trmb2,  s_trmb3,       &
+       zxrugs,    zu10m,     zv10m,    fder_print,    &
+       zxqsurf,   rh2m,      zxfluxu,  zxfluxv,       &
+       rugos_d,   agesno_d,  sollw,    solsw,         &
+       d_ts,      evap_d,    fluxlat,  t2m,           &
+       wfbils,    wfbilo,    flux_t,   flux_u, flux_v,&
+       dflux_t,   dflux_q,   zxsnow,                  &
+       zxfluxt,   zxfluxq,   q2m,      flux_q )
+!****************************************************************************************
+! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
+! Objet: interface de "couche limite" (diffusion verticale)
+!
+!AA REM:
+!AA-----
+!AA Tout ce qui a trait au traceurs est dans phytrac maintenant
+!AA pour l'instant le calcul de la couche limite pour les traceurs
+!AA se fait avec cltrac et ne tient pas compte de la differentiation
+!AA des sous-fraction de sol.
+!AA REM bis :
+!AA----------
+!AA Pour pouvoir extraire les coefficient d'echanges et le vent 
+!AA dans la premiere couche, 3 champs supplementaires ont ete crees
+!AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
+!AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir 
+!AA si les informations des subsurfaces doivent etre prises en compte
+!AA il faudra sortir ces memes champs en leur ajoutant une dimension, 
+!AA c'est a dire nbsrf (nbre de subsurface).
+!
+! Arguments:
+!
+! dtime----input-R- interval du temps (secondes)
+! itap-----input-I- numero du pas de temps
+! date0----input-R- jour initial
+! t--------input-R- temperature (K)
+! q--------input-R- vapeur d'eau (kg/kg)
+! u--------input-R- vitesse u
+! v--------input-R- vitesse v
+! ts-------input-R- temperature du sol (en Kelvin)
+! paprs----input-R- pression a intercouche (Pa)
+! pplay----input-R- pression au milieu de couche (Pa)
+! rlat-----input-R- latitude en degree
+! rugos----input-R- longeur de rugosite (en m)
+!
+! d_t------output-R- le changement pour "t"
+! d_q------output-R- le changement pour "q"
+! d_u------output-R- le changement pour "u"
+! d_v------output-R- le changement pour "v"
+! d_ts-----output-R- le changement pour "ts"
+! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
+!                    (orientation positive vers le bas)
+! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
+! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
+! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
+! dflux_t--output-R- derive du flux sensible
+! dflux_q--output-R- derive du flux latent
+! zu1------output-R- le vent dans la premiere couche
+! zv1------output-R- le vent dans la premiere couche
+! trmb1----output-R- deep_cape
+! trmb2----output-R- inhibition 
+! trmb3----output-R- Point Omega
+! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL
+! plcl-----output-R- Niveau de condensation
+! pblh-----output-R- HCL
+! pblT-----output-R- T au nveau HCL
+!
+    INCLUDE "indicesol.inc"
+    INCLUDE "dimsoil.h"
+    INCLUDE "YOMCST.inc"
+    INCLUDE "iniprint.h"
+    INCLUDE "FCTTRE.inc"
+    INCLUDE "clesphys.inc"
+    INCLUDE "compbl.h"
+    INCLUDE "dimensions90.h"
+    INCLUDE "YOETHF.inc"
+    INCLUDE "temps.inc"
+    INCLUDE "control.inc"
+
+! Input variables
+!****************************************************************************************
+    REAL,                         INTENT(IN)        :: dtime
+    REAL,                         INTENT(IN)        :: date0
+    INTEGER,                      INTENT(IN)        :: itap
+    INTEGER,                      INTENT(IN)        :: jour    ! jour de l'annee en cours
+    LOGICAL,                      INTENT(IN)        :: debut, lafin
+    REAL, DIMENSION(klon),        INTENT(IN)        :: rlon, rlat
+    REAL, DIMENSION(klon),        INTENT(IN)        :: rugoro
+    REAL, DIMENSION(klon),        INTENT(IN)        :: rmu0    ! cosinus de l'angle solaire zenithal
+    REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f, snow_f
+    REAL, DIMENSION(klon),        INTENT(IN)        :: solsw_m ! mean value
+    REAL, DIMENSION(klon),        INTENT(IN)        :: sollw_m ! mean value
+    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t, q
+    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: u, v
+    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pplay
+    REAL, DIMENSION(klon,klev+1), INTENT(IN)        :: paprs
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN)        :: pctsrf
+
+! Input/Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: albe
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: alblw
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: u10m, v10m
+
+! Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: sollwdown
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh, cdragm
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu1
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv1
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: albsol
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: albsollw
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens, zxevap
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxtsol
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zt2m
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsat2m
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t, d_q
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_u, d_v
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zcoefh
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: pctsrf_new
+
+! Output only for diagnostics
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsol_d
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zq2m
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_capCL
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_oliqCL
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_cteiCL
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblT
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_therm
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb1
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb2
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb3
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxrugs
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu10m
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv10m
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: fder_print
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxqsurf
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: rh2m
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxu, zxfluxv
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: rugos_d
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: agesno_d
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: sollw, solsw
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: d_ts
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: evap_d
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: fluxlat
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfbils, wfbilo
+    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t
+    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u, flux_v
+
+! Output not needed
+    REAL, DIMENSION(klon),       INTENT(OUT)        :: dflux_t, dflux_q
+    REAL, DIMENSION(klon),       INTENT(OUT)        :: zxsnow
+    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxt, zxfluxq
+    REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m
+    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q
+
+
+! Local variables with attribute SAVE
+!****************************************************************************************
+    INTEGER, SAVE                            :: nhoridbg, nidbg
+!$OMP THREADPRIVATE(nhoridbg, nidbg)
+    LOGICAL, SAVE                            :: debugindex=.FALSE.
+!$OMP THREADPRIVATE(debugindex)
+    LOGICAL, SAVE                            :: first_call=.TRUE.
+!$OMP THREADPRIVATE(first_call)
+    CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
+!$OMP THREADPRIVATE(cl_surf)
+
+! Other local variables
+!****************************************************************************************
+    INTEGER                            :: i, k, nsrf 
+    INTEGER                            :: knon, j
+    INTEGER                            :: idayref
+    INTEGER , DIMENSION(klon)          :: ni
+    REAL                               :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
+    REAL                               :: amn, amx
+    REAL, DIMENSION(klon)              :: r_co2_ppm     ! taux CO2 atmosphere
+    REAL, DIMENSION(klon)              :: yts, yrugos, ypct, yz0_new
+    REAL, DIMENSION(klon)              :: yalb
+    REAL, DIMENSION(klon)              :: yalblw
+    REAL, DIMENSION(klon)              :: yu1, yv1
+    REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
+    REAL, DIMENSION(klon)              :: yrain_f, ysnow_f
+    REAL, DIMENSION(klon)              :: ysollw, ysolsw, ysollwdown
+    REAL, DIMENSION(klon)              :: yfder
+    REAL, DIMENSION(klon)              :: yrads,yrugoro
+    REAL, DIMENSION(klon)              :: yfluxlat
+    REAL, DIMENSION(klon)              :: y_d_ts
+    REAL, DIMENSION(klon)              :: y_flux_t1, y_flux_q1
+    REAL, DIMENSION(klon)              :: y_dflux_t, y_dflux_q
+    REAL, DIMENSION(klon)              :: u1lay, v1lay
+    REAL, DIMENSION(klon)              :: yt2m, yq2m, yu10m
+    REAL, DIMENSION(klon)              :: yustar
+    REAL, DIMENSION(klon)              :: yu10mx
+    REAL, DIMENSION(klon)              :: yu10my
+    REAL, DIMENSION(klon)              :: ywindsp
+    REAL, DIMENSION(klon)              :: yt10m, yq10m
+    REAL, DIMENSION(klon)              :: ypblh
+    REAL, DIMENSION(klon)              :: ylcl
+    REAL, DIMENSION(klon)              :: ycapCL
+    REAL, DIMENSION(klon)              :: yoliqCL
+    REAL, DIMENSION(klon)              :: ycteiCL
+    REAL, DIMENSION(klon)              :: ypblT
+    REAL, DIMENSION(klon)              :: ytherm
+    REAL, DIMENSION(klon)              :: ytrmb1
+    REAL, DIMENSION(klon)              :: ytrmb2
+    REAL, DIMENSION(klon)              :: ytrmb3
+    REAL, DIMENSION(klon)              :: uzon, vmer
+    REAL, DIMENSION(klon)              :: tair1, qair1, tairsol
+    REAL, DIMENSION(klon)              :: psfce, patm
+    REAL, DIMENSION(klon)              :: qairsol, zgeo1
+    REAL, DIMENSION(klon)              :: rugo1
+    REAL, DIMENSION(klon)              :: yfluxsens, swdown
+    REAL, DIMENSION(klon)              :: petAcoef, peqAcoef, petBcoef, peqBcoef
+    REAL, DIMENSION(klon)              :: ypsref, epot_air
+    REAL, DIMENSION(klon)              :: yevap, ytsurf_new, yalb_new
+    REAL, DIMENSION(klon)              :: pctsrf_nsrf
+    REAL, DIMENSION(klon)              :: ztsol
+    REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q
+    REAL, DIMENSION(klon,klev)         :: y_d_u, y_d_v
+    REAL, DIMENSION(klon,klev)         :: y_flux_t, y_flux_q
+    REAL, DIMENSION(klon,klev)         :: y_flux_u, y_flux_v
+    REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm
+    REAL, DIMENSION(klon,klev)         :: yu, yv
+    REAL, DIMENSION(klon,klev)         :: yt, yq
+    REAL, DIMENSION(klon,klev)         :: ypplay, ydelp
+    REAL, DIMENSION(klon,klev)         :: delp
+    REAL, DIMENSION(klon,klev+1)       :: ypaprs
+    REAL, DIMENSION(klon,nsoilmx)      :: ytsoil
+    REAL, DIMENSION(klon,nbsrf)        :: pctsrf_pot
+    CHARACTER(len=80)                  :: abort_message
+    CHARACTER(len=20)                  :: modname = 'pbl_surface'
+    LOGICAL, PARAMETER                 :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
+    LOGICAL, PARAMETER                 :: check=.FALSE.
+
+! For debugging with IOIPSL
+    INTEGER, DIMENSION(iim*(jjm+1))    :: ndexbg
+    REAL                               :: zjulian
+    REAL, DIMENSION(klon)              :: tabindx
+    REAL, DIMENSION(iim,jjm+1)         :: zx_lon, zx_lat
+    REAL, DIMENSION(iim,jjm+1)         :: debugtab
+
+
+    REAL, DIMENSION(klon,nbsrf)        :: pblh
+    REAL, DIMENSION(klon,nbsrf)        :: plcl
+    REAL, DIMENSION(klon,nbsrf)        :: capCL
+    REAL, DIMENSION(klon,nbsrf)        :: oliqCL
+    REAL, DIMENSION(klon,nbsrf)        :: cteiCL
+    REAL, DIMENSION(klon,nbsrf)        :: pblT
+    REAL, DIMENSION(klon,nbsrf)        :: therm
+    REAL, DIMENSION(klon,nbsrf)        :: trmb1
+    REAL, DIMENSION(klon,nbsrf)        :: trmb2
+    REAL, DIMENSION(klon,nbsrf)        :: trmb3
+    REAL, DIMENSION(klon,nbsrf)        :: zx_rh2m, zx_qsat2m
+    REAL, DIMENSION(klon,nbsrf)        :: zx_qs1, zx_t1
+    REAL, DIMENSION(klon,nbsrf)        :: zdelta1, zcor1
+
+
+!jg+ temporary test
+    REAL, DIMENSION(klon,klev)         :: y_flux_u_old, y_flux_v_old
+    REAL, DIMENSION(klon,klev)         :: y_d_u_old, y_d_v_old
+!jg-
+    
+!****************************************************************************************
+! End of declarations
+!****************************************************************************************
+
+
+!****************************************************************************************
+! 1) Initialisation and validation tests 
+!    Only done first time entering this subroutine
+!
+!****************************************************************************************
+
+
+    IF (first_call) THEN
+       first_call=.FALSE.
+      
+       ! Initilize debug IO
+       IF (debugindex .AND. mpi_size==1) THEN 
+          ! initialize IOIPSL output
+          idayref = day_ini
+          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
+          CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
+          DO i = 1, iim
+             zx_lon(i,1) = rlon(i+1)
+             zx_lon(i,jjm+1) = rlon(i+1)
+          ENDDO
+          CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
+          CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), &
+               1,iim,1,jjm+1, &
+               itau_phy,zjulian,dtime,nhoridbg,nidbg) 
+          ! no vertical axis
+          cl_surf(1)='ter'
+          cl_surf(2)='lic'
+          cl_surf(3)='oce'
+          cl_surf(4)='sic'
+          DO nsrf=1,nbsrf
+             CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, &
+                  jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime) 
+          END DO
+
+          CALL histend(nidbg)
+          CALL histsync(nidbg)
+
+       END IF
+       
+    ENDIF
+          
+!****************************************************************************************
+! 2) Initialization to zero 
+!    Done for all local variables that will be compressed later
+!    and argument with INTENT(OUT)
+!****************************************************************************************
+    cdragh = 0.0  ; cdragm = 0.0     ; dflux_t = 0.0   ; dflux_q = 0.0
+    ypct = 0.0    ; yts = 0.0        ; ysnow = 0.0     ; zu1 = 0.0        
+    zv1 = 0.0     ; yqsurf = 0.0     ; yalb = 0.0      ; yalblw = 0.0    
+    yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0    
+    ysollw = 0.0  ; ysollwdown = 0.0 ; yrugos = 0.0    ; yu1 = 0.0    
+    yv1 = 0.0     ; yrads = 0.0      ; ypaprs = 0.0    ; ypplay = 0.0
+    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
+    yq = 0.0      ; pctsrf_new = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0 
+    yrugoro = 0.0 ; yu10mx = 0.0     ; yu10my = 0.0    ; ywindsp = 0.0   
+    d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
+    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0      
+    d_u = 0.0     ; d_v = 0.0        ; zcoefh = 0.0    ; yqsol = 0.0    
+    ytherm = 0.0
+     
+    ytsoil = 999999. 
+
+!****************************************************************************************
+! 3) - Calculate pressure thickness of each layer
+!    - Calculate the wind at first layer
+!
+!****************************************************************************************
+    DO k = 1, klev
+       DO i = 1, klon
+          delp(i,k) = paprs(i,k)-paprs(i,k+1)
+       ENDDO
+    ENDDO
+    DO i = 1, klon
+       zx_alf1 = 1.0
+       zx_alf2 = 1.0 - zx_alf1
+       u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2
+       v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2
+    ENDDO
+
+
+!****************************************************************************************
+! Test for rugos........ from physiq.. A la fin plutot???
+! Calcul de l'abedo moyen par maille
+!****************************************************************************************
+
+    zxrugs(:) = 0.0
+    DO nsrf = 1, nbsrf
+       DO i = 1, klon
+          rugos(i,nsrf) = MAX(rugos(i,nsrf),0.000015)
+          zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf)
+       ENDDO
+    ENDDO
+
+! Calcul de l'abedo moyen par maille
+    albsol(:)   = 0.0
+    albsollw(:) = 0.0
+    DO nsrf = 1, nbsrf
+       DO i = 1, klon
+          albsol(i)   = albsol(i)   + albe(i,nsrf)  * pctsrf(i,nsrf)
+          albsollw(i) = albsollw(i) + alblw(i,nsrf) * pctsrf(i,nsrf)
+       ENDDO
+    ENDDO
+
+
+
+! Calcule de ztsol (aussi fait dans physiq.F, pourrait etre un argument)
+    ztsol(:) = 0.0
+    DO nsrf = 1, nbsrf
+       DO i = 1, klon
+          ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
+       ENDDO
+    ENDDO
+
+
+! Repartition du longwave par sous-surface linearisee
+    DO nsrf = 1, nbsrf
+       DO i = 1, klon
+          sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))
+          solsw(i,nsrf) = solsw_m(i)*(1.-albe(i,nsrf))/(1.-albsol(i))
+       ENDDO
+    ENDDO
+
+
+    DO i = 1, klon
+       sollwdown(i) = sollw_m(i) + RSIGMA*ztsol(i)**4
+    ENDDO
+
+!****************************************************************************************
+! 4) Loop over different surfaces
+!
+! All points with a possibility of the current surface are used. This is done
+! to allow the sea-ice to appear or disappear. It is considered here that the 
+! entier domaine of the ocean possibly can contain sea-ice.
+! 
+!****************************************************************************************
+
+    pctsrf_pot = pctsrf
+    pctsrf_pot(:,is_oce) = 1. - zmasq(:)
+    pctsrf_pot(:,is_sic) = 1. - zmasq(:)
+      
+    loop_nbsrf: DO nsrf = 1, nbsrf
+
+! Search for index(ni) and size(knon) of domaine to treat
+       ni(:) = 0
+       knon  = 0
+       DO i = 1, klon
+          IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN
+             knon = knon + 1
+             ni(knon) = i
+          ENDIF
+       ENDDO
+
+       ! write index, with IOIPSL
+       IF (debugindex .AND. mpi_size==1) THEN 
+          tabindx(:)=0.
+          DO i=1,knon
+             tabindx(i)=FLOAT(i)
+          END DO
+          debugtab(:,:) = 0.
+          ndexbg(:) = 0
+          CALL gath2cpl(tabindx,debugtab,knon,ni)
+          CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg)
+       ENDIF
+       
+!****************************************************************************************
+! 5) Compress variables 
+!
+!****************************************************************************************
+
+       DO j = 1, knon
+          i = ni(j)
+          ypct(j) = pctsrf(i,nsrf)
+          yts(j) = ts(i,nsrf)
+          ysnow(j) = snow(i,nsrf)
+          yqsurf(j) = qsurf(i,nsrf)
+          yalb(j) = albe(i,nsrf)
+          yalblw(j) = alblw(i,nsrf)
+          yrain_f(j) = rain_f(i)
+          ysnow_f(j) = snow_f(i)
+          yagesno(j) = agesno(i,nsrf)
+          yfder(j) = fder(i)
+          ysolsw(j) = solsw(i,nsrf)
+          ysollw(j) = sollw(i,nsrf)
+          ysollwdown(j) = sollwdown(i)
+          yrugos(j) = rugos(i,nsrf)
+          yrugoro(j) = rugoro(i)
+          yu1(j) = u1lay(i)
+          yv1(j) = v1lay(i)
+          yrads(j) =  ysolsw(j)+ ysollw(j)
+          ypaprs(j,klev+1) = paprs(i,klev+1)
+          yu10mx(j) = u10m(i,nsrf)
+          yu10my(j) = v10m(i,nsrf)
+          ywindsp(j) = SQRT(yu10mx(j)*yu10mx(j) + yu10my(j)*yu10my(j) )
+       END DO
+
+       DO k = 1, klev
+          DO j = 1, knon
+             i = ni(j)
+             ypaprs(j,k) = paprs(i,k)
+             ypplay(j,k) = pplay(i,k)
+             ydelp(j,k) = delp(i,k)
+             yu(j,k) = u(i,k)
+             yv(j,k) = v(i,k)
+             yt(j,k) = t(i,k)
+             yq(j,k) = q(i,k)
+          ENDDO
+       ENDDO
+       
+       DO k = 1, nsoilmx
+          DO j = 1, knon
+             i = ni(j)
+             ytsoil(j,k) = ftsoil(i,k,nsrf)
+          END DO
+       END DO
+       
+       ! qsol(water height in soil) only for bucket continental model
+       IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN 
+          DO j = 1, knon
+             i = ni(j)
+             yqsol(j) = qsol(i)
+          END DO
+       ENDIF
+       
+!****************************************************************************************
+! 6) Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the 
+!    atmosphere and coefficients for turbulent diffusion at surface(Cdrag).
+!
+!****************************************************************************************
+
+       CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
+            ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf,  &
+            ycoefm, ycoefh)
+       
+!jg+
+!****************************************************************************************
+!    => Old method
+! Calculer la diffusion des vitesses "u" et "v"   
+! Output can be used : y_d_u_old, y_flux_u_old, y_d_v_old, y_flux_v_old
+!
+!****************************************************************************************
+       CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp, &
+            y_d_u_old, y_flux_u_old)
+       CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp, &
+            y_d_v_old, y_flux_v_old)
+       ! save_flux est utile pour pouvoir utilise calcul_flux_vent plus tard
+       CALL save_flux(klon, y_flux_u_old(:,1), y_flux_v_old(:,1))
+!jg-
+
+!****************************************************************************************
+! 
+! 8) "La descente" - "The downhill"
+!  
+!  climb_hq_down and climb_wind_down calculate the coefficients
+!  Ccoef_X et Dcoef_X for X=[H, Q, U, V].
+!  Only the coefficients at surface for H and Q are returned.
+!
+!****************************************************************************************
+
+! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q 
+       CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
+            ydelp, yt, yq, dtime, &
+            petAcoef, peqAcoef, petBcoef, peqBcoef)
+
+! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
+       CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv)
+      
+
+!****************************************************************************************
+! 9) Small calculations
+!
+!****************************************************************************************
+          
+       ypsref(:) = ypaprs(:,1)  
+       epot_air(:) = 0.0
+       epot_air(1:knon) = RCPD*yt(1:knon,1)*(ypsref(1:knon)/ypplay(1:knon,1))**RKAPPA
+
+       swdown(:) = 0.0
+       IF (nsrf .EQ. is_ter) THEN 
+          swdown(1:knon) = ysolsw(1:knon)/(1-yalb(1:knon)) 
+       ELSE 
+          swdown(1:knon) = ysolsw(1:knon)
+       ENDIF
+
+       ! constant CO2
+       r_co2_ppm(:) = co2_ppm
+
+!****************************************************************************************
+!
+! 10) Switch selon current surface
+!     It is necessary to start with the continental surfaces because the ocean
+!     needs their run-off.
+!
+!****************************************************************************************
+       SELECT CASE(nsrf)
+     
+       CASE(is_ter)
+          CALL surf_land(itap, dtime, date0, jour, knon, ni,&
+               rlon, rlat, &
+               debut, lafin, ydelp(:,1), epot_air, r_co2_ppm, ysollwdown, ysolsw, swdown, &
+               yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+               petAcoef, peqAcoef, petBcoef, peqBcoef, & 
+               ypsref, yu1, yv1, yrugoro, pctsrf, &
+               yrads, ysnow, yqsurf, yqsol, yagesno, &
+               ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
+               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
+     
+       CASE(is_lic)
+          CALL surf_landice(itap, dtime, knon, ni, &
+               yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+               petAcoef, peqAcoef, petBcoef, peqBcoef, &
+               ypsref, yu1, yv1, yrugoro, pctsrf, &
+               yrads, ysnow, yqsurf, yqsol, yagesno, &
+               ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
+               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
+          
+       CASE(is_oce)
+          CALL surf_ocean(rlon, rlat, ysollw, yalb, &
+               yrugos, ywindsp, rmu0, & 
+               yfder, &
+               itap, dtime, jour, knon, ni, &
+               debut, swdown, &
+               ypplay(:,1), ycoefh(:,1), ycoefm(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+               petAcoef, peqAcoef, petBcoef, peqBcoef, &
+               ypsref, yu1, yv1, yrugoro, pctsrf, &
+               yrads, ysnow, yqsurf, yagesno, &
+               yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
+               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
+          
+       CASE(is_sic)
+          CALL surf_seaice( &
+               rlon, rlat, ysollw, yalb, &
+               yfder, &
+               itap, dtime, jour, knon, ni, &
+               debut, lafin, swdown, &
+               yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+               petAcoef, peqAcoef, petBcoef, peqBcoef, &
+               ypsref, yu1, yv1, yrugoro, pctsrf, &
+               yrads, ysnow, yqsurf, yqsol, yagesno, &
+               ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, &
+               ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
+          
+
+       CASE DEFAULT
+          WRITE(lunout,*) 'Surface index = ', nsrf
+          abort_message = 'Surface index not valid'
+          CALL abort_gcm(modname,abort_message,1)
+       END SELECT
+
+!****************************************************************************************
+! Save the fraction of this sub-surface 
+!
+!****************************************************************************************
+       pctsrf_new(:,nsrf) = pctsrf_nsrf(:)
+
+!****************************************************************************************
+! 11) - Calcul the increment of surface temperature
+!     - Update albedo
+!
+!****************************************************************************************
+       y_d_ts(1:knon)   = ytsurf_new(1:knon) - yts(1:knon)
+ 
+       yalb(1:knon) = yalb_new(1:knon)
+
+!****************************************************************************************
+!
+! 12) "La remontee" - "The uphill"
+!
+!  The fluxes (y_flux_X) and tendancy (y_d_X) are calculated 
+!  for X=H, Q, U and V, for all vertical levels.
+!
+!****************************************************************************************
+! H and Q
+       y_flux_t1(:) =  yfluxsens(:)
+       y_flux_q1(:) = -yevap(:)
+
+       CALL climb_hq_up(knon, dtime, yt, yq, &
+            y_flux_q1, y_flux_t1, ypaprs, ypplay, &
+            y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))    
+       
+! U and V
+       CALL climb_wind_up(knon, dtime, yu, yv, &
+            y_flux_u, y_flux_v, y_d_u, y_d_v)
+
+!jg+ temporary for testing 
+! Use the results from old method
+       IF (.NOT. newwind) THEN
+          y_flux_u(:,:) = y_flux_u_old(:,:)
+          y_flux_v(:,:) = y_flux_v_old(:,:)
+          y_d_u(:,:)    = y_d_u_old(:,:)
+          y_d_v(:,:)    = y_d_v_old(:,:)
+       ENDIF
+!jg-
+   
+       DO j = 1, knon
+          y_dflux_t(j) = y_dflux_t(j) * ypct(j)
+          y_dflux_q(j) = y_dflux_q(j) * ypct(j)
+          yu1(j) = yu1(j) *  ypct(j)
+          yv1(j) = yv1(j) *  ypct(j)
+       ENDDO
+
+!****************************************************************************************
+! 13) Transform variables for output format : 
+!     - Decompress
+!     - Multiply with pourcentage of current surface
+!     - Cumulate in global variable
+!
+!****************************************************************************************
+
+       DO k = 1, klev
+          DO j = 1, knon
+             i = ni(j)
+             ycoefh(j,k) = ycoefh(j,k) * ypct(j)
+             ycoefm(j,k) = ycoefm(j,k) * ypct(j)
+             y_d_t(j,k) = y_d_t(j,k) * ypct(j)
+             y_d_q(j,k) = y_d_q(j,k) * ypct(j)
+             y_d_u(j,k) = y_d_u(j,k) * ypct(j)
+             y_d_v(j,k) = y_d_v(j,k) * ypct(j)
+
+             flux_t(i,k,nsrf) = y_flux_t(j,k)
+             flux_q(i,k,nsrf) = y_flux_q(j,k)
+             flux_u(i,k,nsrf) = y_flux_u(j,k)
+             flux_v(i,k,nsrf) = y_flux_v(j,k)
+          ENDDO
+       ENDDO
+       
+       evap(:,nsrf) = - flux_q(:,1,nsrf)
+       
+       albe(:, nsrf) = 0.
+       alblw(:, nsrf) = 0.
+       snow(:, nsrf) = 0.
+       qsurf(:, nsrf) = 0.
+       rugos(:, nsrf) = 0.
+       fluxlat(:,nsrf) = 0.
+       DO j = 1, knon
+          i = ni(j)
+          d_ts(i,nsrf) = y_d_ts(j)
+          albe(i,nsrf) = yalb(j)  
+          alblw(i,nsrf) = yalblw(j)
+          snow(i,nsrf) = ysnow(j)  
+          qsurf(i,nsrf) = yqsurf(j)
+          rugos(i,nsrf) = yz0_new(j)
+          fluxlat(i,nsrf) = yfluxlat(j)
+          agesno(i,nsrf) = yagesno(j)  
+          cdragh(i) = cdragh(i) + ycoefh(j,1)
+          cdragm(i) = cdragm(i) + ycoefm(j,1)
+          dflux_t(i) = dflux_t(i) + y_dflux_t(j)
+          dflux_q(i) = dflux_q(i) + y_dflux_q(j)
+          zu1(i) = zu1(i) + yu1(j)
+          zv1(i) = zv1(i) + yv1(j)
+       END DO
+
+       IF ( nsrf .EQ. is_ter ) THEN 
+          DO j = 1, knon
+             i = ni(j)
+             qsol(i) = yqsol(j)
+          END DO
+       END IF
+       
+       ftsoil(:,:,nsrf) = 0.
+       DO k = 1, nsoilmx
+          DO j = 1, knon
+             i = ni(j)
+             ftsoil(i, k, nsrf) = ytsoil(j,k)
+          END DO
+       END DO
+       
+       
+#ifdef CRAY
+       DO k = 1, klev
+          DO j = 1, knon
+             i = ni(j)
+#else
+       DO j = 1, knon
+          i = ni(j)
+          DO k = 1, klev
+#endif
+             d_t(i,k) = d_t(i,k) + y_d_t(j,k)
+             d_q(i,k) = d_q(i,k) + y_d_q(j,k)
+             d_u(i,k) = d_u(i,k) + y_d_u(j,k)
+             d_v(i,k) = d_v(i,k) + y_d_v(j,k)
+             zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
+#ifdef CRAY
+          END DO
+       END DO
+#else
+          END DO
+       END DO
+#endif
+
+!****************************************************************************************
+! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m 
+!     Call HBTM
+!
+!****************************************************************************************
+       t2m(:,nsrf)    = 0.
+       q2m(:,nsrf)    = 0.
+       u10m(:,nsrf)   = 0.
+       v10m(:,nsrf)   = 0.
+
+       pblh(:,nsrf)   = 0.        ! Hauteur de couche limite
+       plcl(:,nsrf)   = 0.        ! Niveau de condensation de la CLA
+       capCL(:,nsrf)  = 0.        ! CAPE de couche limite
+       oliqCL(:,nsrf) = 0.        ! eau_liqu integree de couche limite
+       cteiCL(:,nsrf) = 0.        ! cloud top instab. crit. couche limite
+       pblt(:,nsrf)   = 0.        ! T a la Hauteur de couche limite
+       therm(:,nsrf)  = 0.
+       trmb1(:,nsrf)  = 0.        ! deep_cape
+       trmb2(:,nsrf)  = 0.        ! inhibition 
+       trmb3(:,nsrf)  = 0.        ! Point Omega
+
+#undef T2m     
+#define T2m     
+#ifdef T2m
+! diagnostic t,q a 2m et u, v a 10m
+
+       DO j=1, knon
+          i = ni(j)
+          uzon(j) = yu(j,1) + y_d_u(j,1)
+          vmer(j) = yv(j,1) + y_d_v(j,1)
+          tair1(j) = yt(j,1) + y_d_t(j,1)
+          qair1(j) = yq(j,1) + y_d_q(j,1)
+          zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
+               * (ypaprs(j,1)-ypplay(j,1))
+          tairsol(j) = yts(j) + y_d_ts(j)
+          rugo1(j) = yrugos(j)
+          IF(nsrf.EQ.is_oce) THEN
+             rugo1(j) = rugos(i,nsrf)
+          ENDIF
+          psfce(j)=ypaprs(j,1)
+          patm(j)=ypplay(j,1)
+          qairsol(j) = yqsurf(j)
+       END DO
+       
+
+! Calculate the temperature et relative humidity at 2m and the wind at 10m 
+       CALL stdlevvar(klon, knon, nsrf, zxli, &
+            uzon, vmer, tair1, qair1, zgeo1, &
+            tairsol, qairsol, rugo1, psfce, patm, &
+            yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
+
+       DO j=1, knon
+          i = ni(j)
+          t2m(i,nsrf)=yt2m(j)
+          
+          q2m(i,nsrf)=yq2m(j)
+
+! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
+          u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
+          v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
+
+       END DO
+
+
+       CALL HBTM(knon, ypaprs, ypplay, &
+            yt2m,yt10m,yq2m,yq10m,yustar, &
+            y_flux_t,y_flux_q,yu,yv,yt,yq, &
+            ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
+            ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
+       
+       DO j=1, knon
+          i = ni(j)
+          pblh(i,nsrf)   = ypblh(j)
+          plcl(i,nsrf)   = ylcl(j)
+          capCL(i,nsrf)  = ycapCL(j)
+          oliqCL(i,nsrf) = yoliqCL(j)
+          cteiCL(i,nsrf) = ycteiCL(j)
+          pblT(i,nsrf)   = ypblT(j)
+          therm(i,nsrf)  = ytherm(j)
+          trmb1(i,nsrf)  = ytrmb1(j)
+          trmb2(i,nsrf)  = ytrmb2(j)
+          trmb3(i,nsrf)  = ytrmb3(j)
+       END DO
+       
+#else 
+! not defined T2m
+! No calculation
+! Set output variables to zero
+       DO j = 1, knon
+          i = ni(j)
+          pblh(i,nsrf)   = 0.
+          plcl(i,nsrf)   = 0.
+          capCL(i,nsrf)  = 0.
+          oliqCL(i,nsrf) = 0.
+          cteiCL(i,nsrf) = 0.
+          pblT(i,nsrf)   = 0.
+          therm(i,nsrf)  = 0.
+          trmb1(i,nsrf)  = 0.
+          trmb2(i,nsrf)  = 0.
+          trmb3(i,nsrf)  = 0.
+       END DO
+       DO j = 1, knon
+          i = ni(j) 
+          t2m(i,nsrf)=0.
+          q2m(i,nsrf)=0.
+          u10m(i,nsrf)=0.
+          v10m(i,nsrf)=0.
+       END DO
+#endif
+
+!****************************************************************************************
+! 15) End of loop over different surfaces
+!
+!****************************************************************************************
+    END DO loop_nbsrf
+
+!****************************************************************************************
+! 16) Calculate the mean value over all sub-surfaces for som variables
+!
+!     NB!!! jg : Pour garder la convergence numerique j'utilise pctsrf_new comme c'etait 
+!     fait dans physiq.F mais ca devrait plutot etre pctsrf???!!!!! A verifier!
+!****************************************************************************************
+    
+    zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
+    zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
+    DO nsrf = 1, nbsrf
+       DO k = 1, klev
+          DO i = 1, klon
+             zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf_new(i,nsrf)
+             zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf_new(i,nsrf)
+             zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf_new(i,nsrf)
+             zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf_new(i,nsrf)
+          END DO
+       END DO
+    END DO
+
+    DO i = 1, klon
+       zxsens(i)     = - zxfluxt(i,1) ! flux de chaleur sensible au sol
+       zxevap(i)     = - zxfluxq(i,1) ! flux d'evaporation au sol
+       fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
+    ENDDO
+   
+
+    DO i = 1, klon
+       IF ( ABS( pctsrf_new(i, is_ter) + pctsrf_new(i, is_lic) + &
+            pctsrf_new(i, is_oce) + pctsrf_new(i, is_sic)  - 1.) .GT. EPSFRA) & 
+            THEN 
+          WRITE(*,*) 'physiq : pb sous surface au point ', i, &
+               pctsrf_new(i, 1 : nbsrf)
+       ENDIF
+    ENDDO
+
+!
+! Incrementer la temperature du sol
+!
+    zxtsol(:) = 0.0  ; zxfluxlat(:) = 0.0
+    zt2m(:) = 0.0    ; zq2m(:) = 0.0 
+    zu10m(:) = 0.0   ; zv10m(:) = 0.0
+    s_pblh(:) = 0.0  ; s_plcl(:) = 0.0 
+    s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
+    s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
+    s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
+    s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
+    
+    
+    DO nsrf = 1, nbsrf
+       DO i = 1, klon          
+          ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
+          
+          wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
+               + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf_new(i,nsrf)
+          wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
+               pctsrf_new(i,nsrf)
+
+          zxtsol(i)    = zxtsol(i)    + ts(i,nsrf)      * pctsrf_new(i,nsrf)
+          zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf_new(i,nsrf)
+          
+          zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf_new(i,nsrf)
+          zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf_new(i,nsrf)
+          zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf_new(i,nsrf)
+          zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf_new(i,nsrf)
+
+          s_pblh(i)   = s_pblh(i)   + pblh(i,nsrf)  * pctsrf_new(i,nsrf)
+          s_plcl(i)   = s_plcl(i)   + plcl(i,nsrf)  * pctsrf_new(i,nsrf)
+          s_capCL(i)  = s_capCL(i)  + capCL(i,nsrf) * pctsrf_new(i,nsrf)
+          s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf_new(i,nsrf)
+          s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf_new(i,nsrf)
+          s_pblT(i)   = s_pblT(i)   + pblT(i,nsrf)  * pctsrf_new(i,nsrf)
+          s_therm(i)  = s_therm(i)  + therm(i,nsrf) * pctsrf_new(i,nsrf)
+          s_trmb1(i)  = s_trmb1(i)  + trmb1(i,nsrf) * pctsrf_new(i,nsrf)
+          s_trmb2(i)  = s_trmb2(i)  + trmb2(i,nsrf) * pctsrf_new(i,nsrf)
+          s_trmb3(i)  = s_trmb3(i)  + trmb3(i,nsrf) * pctsrf_new(i,nsrf)
+       END DO
+    END DO
+
+    IF (check) THEN
+       amn=MIN(ts(1,is_ter),1000.)
+       amx=MAX(ts(1,is_ter),-1000.)
+       DO i=2, klon
+          amn=MIN(ts(i,is_ter),amn)
+          amx=MAX(ts(i,is_ter),amx)
+       ENDDO
+       PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
+    ENDIF
+!
+! If a sub-surface does not exsist for a grid point, the mean value for all 
+! sub-surfaces is distributed.
+!
+    DO nsrf = 1, nbsrf
+       DO i = 1, klon
+          IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
+             ts(i,nsrf)     = zxtsol(i)
+             t2m(i,nsrf)    = zt2m(i)
+             q2m(i,nsrf)    = zq2m(i)
+             u10m(i,nsrf)   = zu10m(i)
+             v10m(i,nsrf)   = zv10m(i)
+
+! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
+             pblh(i,nsrf)   = s_pblh(i)
+             plcl(i,nsrf)   = s_plcl(i)
+             capCL(i,nsrf)  = s_capCL(i)
+             oliqCL(i,nsrf) = s_oliqCL(i) 
+             cteiCL(i,nsrf) = s_cteiCL(i)
+             pblT(i,nsrf)   = s_pblT(i)
+             therm(i,nsrf)  = s_therm(i)
+             trmb1(i,nsrf)  = s_trmb1(i)
+             trmb2(i,nsrf)  = s_trmb2(i)
+             trmb3(i,nsrf)  = s_trmb3(i)
+          ENDIF
+       ENDDO
+    ENDDO
+
+
+    DO i = 1, klon
+       fder(i) = - 4.0*RSIGMA*zxtsol(i)**3 
+    ENDDO
+    
+    zxqsurf(:) = 0.0
+    zxsnow(:)  = 0.0
+    DO nsrf = 1, nbsrf
+       DO i = 1, klon
+          zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf_new(i,nsrf)
+          zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf_new(i,nsrf)
+       END DO
+    END DO
+
+!
+!IM Calculer l'humidite relative a 2m (rh2m) pour diagnostique
+!IM ajout dependance type surface
+    rh2m(:)   = 0.0
+    qsat2m(:) = 0.0
+
+    DO i = 1, klon
+       DO nsrf=1, nbsrf
+          zx_t1(i,nsrf) = t2m(i,nsrf)
+          IF (thermcep) THEN
+             zdelta1(i,nsrf) = MAX(0.,SIGN(1.,rtt-zx_t1(i,nsrf)))
+             zx_qs1(i,nsrf)  = r2es * &
+                  FOEEW(zx_t1(i,nsrf),zdelta1(i,nsrf))/paprs(i,1)
+             zx_qs1(i,nsrf)  = MIN(0.5,zx_qs1(i,nsrf))
+             zcor1(i,nsrf)   = 1./(1.-retv*zx_qs1(i,nsrf))
+             zx_qs1(i,nsrf)  = zx_qs1(i,nsrf)*zcor1(i,nsrf)
+          END IF
+          zx_rh2m(i,nsrf) = q2m(i,nsrf)/zx_qs1(i,nsrf)
+          zx_qsat2m(i,nsrf)=zx_qs1(i,nsrf)
+          rh2m(i) = rh2m(i)+zx_rh2m(i,nsrf)*pctsrf_new(i,nsrf)
+          qsat2m(i)=qsat2m(i)+zx_qsat2m(i,nsrf)*pctsrf_new(i,nsrf)
+       END DO
+    END DO
+
+! Some of the module declared variables are returned for printing in physiq.F
+    qsol_d(:)     = qsol(:)
+    evap_d(:,:)   = evap(:,:)
+    rugos_d(:,:)  = rugos(:,:) 
+    agesno_d(:,:) = agesno(:,:)
+
+
+  END SUBROUTINE pbl_surface
+!
+!****************************************************************************************
+!
+  SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, &
+       evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
+
+    INCLUDE "indicesol.inc"
+    INCLUDE "dimsoil.h"
+
+! Ouput variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)                 :: qsol_rst
+    REAL, DIMENSION(klon), INTENT(OUT)                 :: fder_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: snow_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: qsurf_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: evap_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: rugos_rst
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: agesno_rst
+    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
+
+ 
+!****************************************************************************************
+! Return module variables for writing to restart file
+!
+!****************************************************************************************    
+    qsol_rst(:)       = qsol(:)
+    fder_rst(:)       = fder(:)
+    snow_rst(:,:)     = snow(:,:)
+    qsurf_rst(:,:)    = qsurf(:,:)
+    evap_rst(:,:)     = evap(:,:)
+    rugos_rst(:,:)    = rugos(:,:)
+    agesno_rst(:,:)   = agesno(:,:)
+    ftsoil_rst(:,:,:) = ftsoil(:,:,:)
+
+!****************************************************************************************
+! Deallocate module variables
+!
+!****************************************************************************************
+    DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
+
+  END SUBROUTINE pbl_surface_final
+!  
+!****************************************************************************************
+!  
+
+END MODULE pbl_surface_mod
Index: /LMDZ4/trunk/libf/phylmd/surf_land_bucket_mod.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/surf_land_bucket_mod.F90	(revision 781)
+++ /LMDZ4/trunk/libf/phylmd/surf_land_bucket_mod.F90	(revision 781)
@@ -0,0 +1,317 @@
+!
+! $Header$
+!
+MODULE surf_land_bucket_mod
+!
+! Surface land bucket module
+!
+! This module is used when no external land model is choosen.
+!
+  USE fonte_neige_mod
+  USE calcul_fluxs_mod
+  USE dimphy
+  USE mod_grid_phy_lmdz
+  USE mod_phys_lmdz_para
+ 
+  IMPLICIT NONE
+
+CONTAINS
+
+  SUBROUTINE surf_land_bucket(itime, jour, knon, knindex, debut, dtime,&
+       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, &
+       spechum, petAcoef, peqAcoef, petBcoef, peqBcoef, ps, &
+       u1_lay, v1_lay, rugoro, &
+       radsol, snow, qsol, agesno, tsoil, &
+       qsurf, z0_new, alblw, evap, fluxsens, fluxlat, &
+       tsurf_new, alb_new, dflux_s, dflux_l)
+
+!****************************************************************************************
+! Bucket calculations for surface. 
+!
+    INCLUDE "clesphys.inc"
+    INCLUDE "indicesol.inc"
+    INCLUDE "dimsoil.h"
+    INCLUDE "YOMCST.inc"
+
+! Input variables  
+!****************************************************************************************
+    INTEGER, INTENT(IN)                     :: itime, jour, knon
+    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
+    LOGICAL, INTENT(IN)                     :: debut
+    REAL, INTENT(IN)                        :: dtime
+    REAL, DIMENSION(klon), INTENT(IN)       :: tsurf
+    REAL, DIMENSION(klon), INTENT(IN)       :: p1lay
+    REAL, DIMENSION(klon), INTENT(IN)       :: tq_cdrag
+    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)       :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)       :: petAcoef, peqAcoef
+    REAL, DIMENSION(klon), INTENT(IN)       :: petBcoef, peqBcoef
+    REAL, DIMENSION(klon), INTENT(IN)       :: ps
+    REAL, DIMENSION(klon), INTENT(IN)       :: u1_lay, v1_lay
+    REAL, DIMENSION(klon), INTENT(IN)       :: rugoro
+
+! In/Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
+    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+
+! Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
+    REAL, DIMENSION(klon), INTENT(OUT)       :: z0_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
+
+! Local variables
+!****************************************************************************************
+    REAL, DIMENSION(klon) :: soilcap, soilflux
+    REAL, DIMENSION(klon) :: cal, beta, dif_grnd, capsol
+    REAL, DIMENSION(klon) :: alb_neig
+    REAL, DIMENSION(klon) :: zfra
+    INTEGER               :: i
+!
+!****************************************************************************************
+
+
+!
+!* Read from limit.nc : albedo(alb_new), length of rugosity(z0_new)
+!
+    CALL interfsur_lim(itime, dtime, jour, &
+         knon, knindex, debut,  &
+         alb_new, z0_new)
+    
+!
+!* Calcultaion of fluxes 
+!
+! calculate constants
+    CALL calbeta(dtime, is_ter, knon, snow, qsol, beta, capsol, dif_grnd)
+       
+! calculate temperature, heat capacity and conduction flux in soil
+    IF (soil_model) THEN 
+       CALL soil(dtime, is_ter, knon, snow, tsurf, tsoil, soilcap, soilflux)
+       DO i=1, knon
+          cal(i) = RCPD / soilcap(i)
+          radsol(i) = radsol(i)  + soilflux(i)
+       END DO
+    ELSE 
+       cal(:) = RCPD * capsol(:)
+    ENDIF
+    
+    CALL calcul_fluxs(knon, is_ter, dtime, &
+         tsurf, p1lay, cal, beta, tq_cdrag, ps, &
+         precip_rain, precip_snow, snow, qsurf,  &
+         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+         petAcoef, peqAcoef, petBcoef, peqBcoef, &
+         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+    
+!
+!* Calculate snow height, run_off, age of snow
+!      
+    CALL fonte_neige( knon, is_ter, knindex, dtime, &
+         tsurf, precip_rain, precip_snow, &
+         snow, qsol, tsurf_new, evap)
+!
+!* Calculate the age of snow
+!
+    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
+    
+    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
+    
+    DO i=1, knon
+       zfra(i) = MAX(0.0,MIN(1.0, snow(i)/(snow(i)+10.0)))
+       alb_new(i)  = alb_neig(i) *zfra(i) + alb_new(i)*(1.0-zfra(i))
+    END DO
+
+    alblw(:) = 0.0
+    alblw(1:knon) = alb_new(1:knon)
+
+!
+!* Calculate the rugosity
+!
+    z0_new = SQRT(z0_new**2+rugoro**2)
+       
+!
+!* End
+!
+  END SUBROUTINE surf_land_bucket
+!
+!****************************************************************************************
+!
+  SUBROUTINE interfsur_lim(itime, dtime, jour, &
+       knon, knindex, debut, &
+       lmt_alb_p, lmt_rug_p)
+    
+! Cette routine sert d'interface entre le modele atmospherique et un fichier
+! de conditions aux limites
+!
+! L. Fairhead 02/2000
+!
+! input:
+!   itime        numero du pas de temps courant
+!   dtime        pas de temps de la physique (en s)
+!   jour         jour a lire dans l'annee
+!   knon         nombre de points dans le domaine a traiter
+!   knindex      index des points de la surface a traiter
+!   debut        logical: 1er appel a la physique (initialisation)
+!
+! output:
+!   lmt_alb_p      Albedo lu 
+!   lmt_rug_p      longueur de rugosite lue
+
+    INCLUDE "netcdf.inc"
+
+! Input variables
+!****************************************************************************************
+    INTEGER, INTENT(IN)                      :: itime
+    REAL   , INTENT(IN)                      :: dtime
+    INTEGER, INTENT(IN)                      :: jour
+    INTEGER, INTENT(IN)                      :: knon
+    INTEGER, DIMENSION(klon_loc), INTENT(IN) :: knindex
+    LOGICAL, INTENT(IN)                      :: debut
+
+! Output variables
+!****************************************************************************************
+    REAL, INTENT(out), DIMENSION(klon_loc)   :: lmt_alb_p
+    REAL, INTENT(out), DIMENSION(klon_loc)   :: lmt_rug_p
+
+! Local variables with attribute SAVE
+!****************************************************************************************
+    INTEGER,SAVE   :: lmt_pas     ! frequence de lecture des conditions limites 
+                                  ! (en pas de physique)
+    !$OMP THREADPRIVATE(lmt_pas) 
+    LOGICAL,SAVE   :: deja_lu_sur ! pour indiquer que le jour a lire a deja
+                                  ! lu pour une surface precedente
+    !$OMP THREADPRIVATE(deja_lu_sur) 
+    INTEGER,SAVE                           :: jour_lu_sur 
+    !$OMP THREADPRIVATE(jour_lu_sur) 
+    CHARACTER (len = 20),SAVE              :: fich ='limit.nc'
+    !$OMP THREADPRIVATE(fich) 
+    LOGICAL,SAVE                           :: check = .FALSE.
+    !$OMP THREADPRIVATE(check) 
+! Champs lus dans le fichier de CL
+    REAL, ALLOCATABLE , SAVE, DIMENSION(:) :: alb_lu_p, rug_lu_p
+    !$OMP THREADPRIVATE(alb_lu_p, rug_lu_p) 
+
+! quelques variables pour netcdf
+    INTEGER ,SAVE                          :: nid, nvarid
+    !$OMP THREADPRIVATE(nid, nvarid) 
+    INTEGER, DIMENSION(2),SAVE             :: start, epais
+    !$OMP THREADPRIVATE(start, epais) 
+
+! Other local variables
+!****************************************************************************************
+    INTEGER                                :: ii, ierr
+    CHARACTER (len = 20)                   :: modname = 'interfsur_lim'
+    CHARACTER (len = 80)                   :: abort_message
+    REAL, DIMENSION(klon_glo)              :: alb_lu
+    REAL, DIMENSION(klon_glo)              :: rug_lu
+
+!
+! End delaration
+!****************************************************************************************
+
+    IF (debut) THEN
+       lmt_pas = NINT(86400./dtime * 1.0) ! pour une lecture une fois par jour
+       jour_lu_sur = jour - 1
+       ALLOCATE(alb_lu_p(klon_loc))
+       ALLOCATE(rug_lu_p(klon_loc))
+    ENDIF
+    
+    IF ((jour - jour_lu_sur) /= 0) deja_lu_sur = .FALSE.
+  
+    IF (check) WRITE(*,*) modname,':: jour_lu_sur, deja_lu_sur', jour_lu_sur, deja_lu_sur 
+    IF (check) WRITE(*,*) modname,':: itime, lmt_pas', itime, lmt_pas
+    IF (check) CALL flush(6)
+
+!    
+! Tester d'abord si c'est le moment de lire le fichier
+!
+    IF (MOD(itime-1, lmt_pas) == 0 .AND. .NOT. deja_lu_sur) THEN
+
+!
+! Ouverture et lecture du fichier
+!
+!$OMP MASTER
+       IF (is_mpi_root) THEN
+          fich = TRIM(fich)
+          IF (check) WRITE(*,*)modname,' ouverture fichier ',fich
+          IF (check) CALL flush(6)
+          ierr = NF_OPEN (fich, NF_NOWRITE,nid)
+          IF (ierr.NE.NF_NOERR) THEN
+             abort_message = 'Pb d''ouverture du fichier de conditions aux limites'
+             CALL abort_gcm(modname,abort_message,1)
+          ENDIF
+! 
+! La tranche de donnees a lire:
+          start(1) = 1
+          start(2) = jour
+          epais(1) = klon_glo
+          epais(2) = 1
+!
+! Lecture albedo
+          ierr = NF_INQ_VARID(nid, 'ALB', nvarid)
+          IF (ierr /= NF_NOERR) THEN
+             abort_message = 'Le champ <ALB> est absent'
+             CALL abort_gcm(modname,abort_message,1)
+          ENDIF
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, alb_lu)
+#else
+          ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, alb_lu)
+#endif
+          IF (ierr /= NF_NOERR) THEN
+             abort_message = 'Lecture echouee pour <ALB>'
+             CALL abort_gcm(modname,abort_message,1)
+          ENDIF
+!
+! Lecture rugosite!
+          ierr = NF_INQ_VARID(nid, 'RUG', nvarid)
+          IF (ierr /= NF_NOERR) THEN
+             abort_message = 'Le champ <RUG> est absent'
+             CALL abort_gcm(modname,abort_message,1)
+          ENDIF
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, rug_lu)
+#else
+          ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, rug_lu)
+#endif
+          IF (ierr /= NF_NOERR) THEN
+             abort_message = 'Lecture echouee pour <RUG>'
+             CALL abort_gcm(modname,abort_message,1)
+          ENDIF
+
+!
+! Fin de lecture
+          ierr = NF_CLOSE(nid)
+
+       ENDIF ! is_mpi_root
+!$OMP END MASTER
+
+       CALL Scatter(alb_lu,alb_lu_p)
+       CALL Scatter(rug_lu,rug_lu_p)
+
+       deja_lu_sur = .TRUE.
+       jour_lu_sur = jour       
+
+    ENDIF
+  
+!
+! Recopie des variables dans les champs de sortie
+!
+    lmt_alb_p(:) = 999999.
+    lmt_rug_p(:) = 999999.
+    DO ii = 1, knon
+       lmt_alb_p(ii) = alb_lu_p(knindex(ii))
+       lmt_rug_p(ii) = rug_lu_p(knindex(ii))
+    ENDDO
+    
+
+  END SUBROUTINE interfsur_lim
+!
+!****************************************************************************************
+!
+END MODULE surf_land_bucket_mod
Index: /LMDZ4/trunk/libf/phylmd/surf_land_mod.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/surf_land_mod.F90	(revision 781)
+++ /LMDZ4/trunk/libf/phylmd/surf_land_mod.F90	(revision 781)
@@ -0,0 +1,142 @@
+!
+! $Header$
+!
+MODULE surf_land_mod
+  
+  USE surface_data, ONLY    : ok_veget
+  USE dimphy
+ 
+#ifdef CPP_VEGET
+  USE surf_land_orchidee_mod
+#endif
+  USE surf_land_bucket_mod
+
+  IMPLICIT NONE
+
+CONTAINS
+!
+!****************************************************************************************
+!  
+  SUBROUTINE surf_land(itime, dtime, date0, jour, knon, knindex, &
+       rlon, rlat, &
+       debut, lafin, zlev, epot_air, ccanopy, sollwdown, swnet, swdown, &
+       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
+       petAcoef, peqAcoef, petBcoef, peqBcoef, &
+       ps, u1_lay, v1_lay, rugoro, pctsrf, &
+       radsol, snow, qsurf, qsol, agesno, &
+       tsoil, z0_new, alblw, evap, fluxsens, fluxlat, &
+       tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_ter)
+
+    INCLUDE "indicesol.inc"
+    INCLUDE "dimsoil.h"
+
+! Input variables  
+!****************************************************************************************
+    INTEGER, INTENT(IN)                     :: itime, jour, knon
+    INTEGER, DIMENSION(klon), INTENT(IN)    :: knindex
+    REAL, INTENT(IN)                        :: date0
+    REAL, DIMENSION(klon), INTENT(IN)       :: rlon, rlat
+    LOGICAL, INTENT(IN)                     :: debut, lafin
+    REAL, INTENT(IN)                        :: dtime
+    REAL, DIMENSION(klon), INTENT(IN)       :: zlev, epot_air, ccanopy
+    REAL, DIMENSION(klon), INTENT(IN)       :: sollwdown, swnet, swdown    
+    REAL, DIMENSION(klon), INTENT(IN)       :: tsurf
+    REAL, DIMENSION(klon), INTENT(IN)       :: p1lay
+    REAL, DIMENSION(klon), INTENT(IN)       :: tq_cdrag
+    REAL, DIMENSION(klon), INTENT(IN)       :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)       :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)       :: petAcoef, peqAcoef
+    REAL, DIMENSION(klon), INTENT(IN)       :: petBcoef, peqBcoef
+    REAL, DIMENSION(klon), INTENT(IN)       :: ps
+    REAL, DIMENSION(klon), INTENT(IN)       :: u1_lay, v1_lay
+    REAL, DIMENSION(klon), INTENT(IN)       :: rugoro
+    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
+
+! In/Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
+    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+
+! Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
+    REAL, DIMENSION(klon), INTENT(OUT)       :: z0_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
+    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_ter
+
+! Local variables
+!****************************************************************************************
+    REAL, DIMENSION(klon) :: ps_tmp, p1lay_tmp
+    REAL, DIMENSION(klon) :: tsol_rad, emis_new ! output from interfsol not used
+    INTEGER               :: i
+
+
+!**************************************************************************************** 
+! Choice between call to vegetation model (ok_veget=true) or simple calculation below
+!
+!****************************************************************************************
+   IF (ok_veget) THEN
+!****************************************************************************************
+!  Call model sechiba
+!
+!****************************************************************************************
+       p1lay_tmp(:)      = 0.0
+       ps_tmp(:)         = 0.0
+       p1lay_tmp(1:knon) = p1lay(1:knon)/100.
+       ps_tmp(1:knon)    = ps(1:knon)/100.
+
+#ifdef CPP_VEGET
+       CALL surf_land_orchidee(itime, dtime, date0, knon, &
+            knindex, rlon, rlat, pctsrf, &
+            debut, lafin, &
+            zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, & 
+            tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
+            precip_rain, precip_snow, sollwdown, swnet, swdown, &
+            ps_tmp, &
+            evap, fluxsens, fluxlat, &              
+            tsol_rad, tsurf_new, alb_new, alblw, &
+            emis_new, z0_new, qsurf)       
+#endif
+
+!  
+!* Ajout de la contribution du relief
+!  
+       DO i=1,knon
+          z0_new(i) = SQRT(z0_new(i)**2 + rugoro(i)**2)
+       ENDDO
+
+    ELSE  ! not ok_veget
+!****************************************************************************************
+! No extern vegetation model choosen, call simple bucket calculations instead.
+!
+!****************************************************************************************
+       CALL surf_land_bucket(itime, jour, knon, knindex, debut, dtime,&
+            tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, &
+            spechum, petAcoef, peqAcoef, petBcoef, peqBcoef, ps, &
+            u1_lay, v1_lay, rugoro, &
+            radsol, snow, qsol, agesno, tsoil, &
+            qsurf, z0_new, alblw, evap, fluxsens, fluxlat, &
+            tsurf_new, alb_new, dflux_s, dflux_l)
+
+    ENDIF ! ok_veget
+
+!****************************************************************************************
+! Return the pourcentage of land in each grid cell, even if not changed in here!
+!
+!****************************************************************************************
+    pctsrf_ter(:) = pctsrf(:,is_ter)
+
+
+  END SUBROUTINE surf_land
+!
+!****************************************************************************************
+!  
+END MODULE surf_land_mod
+!
+!****************************************************************************************
+!  
Index: /LMDZ4/trunk/libf/phylmd/surf_land_orchidee_mod.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/surf_land_orchidee_mod.F90	(revision 781)
+++ /LMDZ4/trunk/libf/phylmd/surf_land_orchidee_mod.F90	(revision 781)
@@ -0,0 +1,651 @@
+!
+! $Header$
+!
+MODULE surf_land_orchidee_mod
+!
+! This module controles the interface towards the model ORCHIDEE
+!
+! Subroutines in this module : surf_land_orchidee
+!                              Init_orchidee_index
+!                              Get_orchidee_communicator
+!                              Init_neighbours
+#ifdef CPP_VEGET
+
+  USE dimphy
+  USE intersurf     ! module d'ORCHIDEE
+  USE cpl_mod,      ONLY : cpl_send_land_fields
+  USE surface_data, ONLY : ocean, ok_veget
+  USE comgeomphy,   ONLY : cuphy, cvphy
+  USE mod_grid_phy_lmdz
+  USE mod_phys_lmdz_para
+
+  IMPLICIT NONE
+
+  PRIVATE
+  PUBLIC  :: surf_land_orchidee
+
+CONTAINS
+!
+!****************************************************************************************
+!  
+  SUBROUTINE surf_land_orchidee(itime, dtime, date0, knon, &
+       knindex, rlon, rlat, pctsrf, &
+       debut, lafin, &
+       plev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, & 
+       tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
+       precip_rain, precip_snow, lwdown, swnet, swdown, &
+       ps, &
+       evap, fluxsens, fluxlat, &              
+       tsol_rad, tsurf_new, alb_new, alblw, &
+       emis_new, z0_new, qsurf)
+!    
+! Cette routine sert d'interface entre le modele atmospherique et le 
+! modele de sol continental. Appel a sechiba
+!
+! L. Fairhead 02/2000
+!
+! input:
+!   itime        numero du pas de temps
+!   dtime        pas de temps de la physique (en s)
+!   nisurf       index de la surface a traiter (1 = sol continental)
+!   knon         nombre de points de la surface a traiter
+!   knindex      index des points de la surface a traiter
+!   rlon         longitudes de la grille entiere
+!   rlat         latitudes de la grille entiere
+!   pctsrf       tableau des fractions de surface de chaque maille
+!   debut        logical: 1er appel a la physique (lire les restart)
+!   lafin        logical: dernier appel a la physique (ecrire les restart)
+!                     (si false calcul simplifie des fluxs sur les continents)
+!   plev         hauteur de la premiere couche (Pa)      
+!   u1_lay       vitesse u 1ere couche
+!   v1_lay       vitesse v 1ere couche
+!   temp_air     temperature de l'air 1ere couche
+!   spechum      humidite specifique 1ere couche
+!   epot_air     temp pot de l'air
+!   ccanopy      concentration CO2 canopee
+!   tq_cdrag     cdrag
+!   petAcoef     coeff. A de la resolution de la CL pour t
+!   peqAcoef     coeff. A de la resolution de la CL pour q
+!   petBcoef     coeff. B de la resolution de la CL pour t
+!   peqBcoef     coeff. B de la resolution de la CL pour q
+!   precip_rain  precipitation liquide
+!   precip_snow  precipitation solide
+!   lwdown       flux IR descendant a la surface
+!   swnet        flux solaire net
+!   swdown       flux solaire entrant a la surface
+!   ps           pression au sol
+!   radsol       rayonnement net aus sol (LW + SW)
+!   
+!
+! output:
+!   evap         evaporation totale
+!   fluxsens     flux de chaleur sensible
+!   fluxlat      flux de chaleur latente
+!   tsol_rad     
+!   tsurf_new    temperature au sol
+!   alb_new      albedo
+!   emis_new     emissivite
+!   z0_new       surface roughness
+!   qsurf        air moisture at surface
+!
+    INCLUDE "indicesol.h"
+    INCLUDE "temps.inc"
+    INCLUDE "YOMCST.inc"
+    INCLUDE "iniprint.h"
+    INCLUDE "dimensions90.h"
+  
+!
+! Parametres d'entree
+!****************************************************************************************
+    INTEGER, INTENT(IN)                       :: itime
+    REAL, INTENT(IN)                          :: dtime
+    REAL, INTENT(IN)                          :: date0
+    INTEGER, INTENT(IN)                       :: knon
+    INTEGER, DIMENSION(klon), INTENT(IN)      :: knindex
+    LOGICAL, INTENT(IN)                       :: debut, lafin
+    REAL, DIMENSION(klon,nbsrf), INTENT(IN)   :: pctsrf
+    REAL, DIMENSION(klon), INTENT(IN)         :: rlon, rlat
+    REAL, DIMENSION(klon), INTENT(IN)         :: plev
+    REAL, DIMENSION(klon), INTENT(IN)         :: u1_lay, v1_lay
+    REAL, DIMENSION(klon), INTENT(IN)         :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)         :: epot_air, ccanopy
+    REAL, DIMENSION(klon), INTENT(IN)         :: tq_cdrag
+    REAL, DIMENSION(klon), INTENT(IN)         :: petAcoef, peqAcoef
+    REAL, DIMENSION(klon), INTENT(IN)         :: petBcoef, peqBcoef
+    REAL, DIMENSION(klon), INTENT(IN)         :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)         :: lwdown, swnet, swdown, ps
+    REAL, DIMENSION(klon)                     :: swdown_vrai
+
+! Parametres de sortie
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)        :: evap, fluxsens, fluxlat, qsurf
+    REAL, DIMENSION(klon), INTENT(OUT)        :: tsol_rad, tsurf_new, alb_new, alblw
+    REAL, DIMENSION(klon), INTENT(OUT)        :: emis_new, z0_new
+
+! Local
+!****************************************************************************************
+    INTEGER                                   :: ij, jj, igrid, ireal, index
+    INTEGER                                   :: error
+    CHARACTER (len = 20)                      :: modname = 'surf_land_orchidee'
+    CHARACTER (len = 80)                      :: abort_message
+    LOGICAL,SAVE                              :: check = .FALSE.
+    !$OMP THREADPRIVATE(check)
+
+! type de couplage dans sechiba
+!  character (len=10)   :: coupling = 'implicit' 
+! drapeaux controlant les appels dans SECHIBA
+!  type(control_type), save   :: control_in
+! Preserved albedo
+    REAL, ALLOCATABLE, DIMENSION(:), SAVE     :: albedo_keep, zlev
+    !$OMP THREADPRIVATE(albedo_keep,zlev)
+! coordonnees geographiques
+    REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: lalo
+    !$OMP THREADPRIVATE(lalo)
+! pts voisins
+    INTEGER,ALLOCATABLE, DIMENSION(:,:), SAVE :: neighbours
+    !$OMP THREADPRIVATE(neighbours)
+! fractions continents
+    REAL,ALLOCATABLE, DIMENSION(:), SAVE      :: contfrac
+    !$OMP THREADPRIVATE(contfrac)
+! resolution de la grille
+    REAL, ALLOCATABLE, DIMENSION (:,:), SAVE  :: resolution
+    !$OMP THREADPRIVATE(resolution)
+
+    REAL, ALLOCATABLE, DIMENSION (:,:), SAVE  :: lon_scat, lat_scat  
+    !$OMP THREADPRIVATE(lon_scat,lat_scat)
+
+    LOGICAL, SAVE                             :: lrestart_read = .TRUE.
+    !$OMP THREADPRIVATE(lrestart_read)
+    LOGICAL, SAVE                             :: lrestart_write = .FALSE.
+    !$OMP THREADPRIVATE(lrestart_write)
+
+    REAL, DIMENSION(knon,2)                   :: albedo_out
+    !$OMP THREADPRIVATE(albedo_out)
+
+! Pb de nomenclature
+    REAL, DIMENSION(klon)                     :: petA_orc, peqA_orc
+    REAL, DIMENSION(klon)                     :: petB_orc, peqB_orc
+! Pb de correspondances de grilles
+    INTEGER, DIMENSION(:), SAVE, ALLOCATABLE  :: ig, jg
+    !$OMP THREADPRIVATE(ig,jg)
+    INTEGER :: indi, indj
+    INTEGER, SAVE, ALLOCATABLE,DIMENSION(:)   :: ktindex
+    !$OMP THREADPRIVATE(ktindex)
+
+! Essai cdrag
+    REAL, DIMENSION(klon)                     :: cdrag
+    INTEGER,SAVE                              :: offset
+    !$OMP THREADPRIVATE(offset)
+
+    REAL, DIMENSION(klon_glo)                 :: rlon_g,rlat_g
+    INTEGER, SAVE                             :: orch_comm
+    !$OMP THREADPRIVATE(orch_comm)
+
+    REAL, ALLOCATABLE, DIMENSION(:), SAVE     :: coastalflow
+    !$OMP THREADPRIVATE(coastalflow)
+    REAL, ALLOCATABLE, DIMENSION(:), SAVE     :: riverflow
+    !$OMP THREADPRIVATE(riverflow)
+!
+! Fin definition
+!****************************************************************************************
+
+    IF (check) WRITE(lunout,*)'Entree ', modname
+    IF (check) WRITE(lunout,*)'ok_veget = ',ok_veget
+  
+! Initialisation
+  
+    IF (debut) THEN
+       ALLOCATE(ktindex(knon))
+       IF ( .NOT. ALLOCATED(albedo_keep)) THEN
+          ALLOCATE(albedo_keep(klon))
+          ALLOCATE(zlev(knon))
+       ENDIF
+! Pb de correspondances de grilles
+       ALLOCATE(ig(klon))
+       ALLOCATE(jg(klon))
+       ig(1) = 1
+       jg(1) = 1
+       indi = 0
+       indj = 2
+       DO igrid = 2, klon - 1
+          indi = indi + 1
+          IF ( indi > iim) THEN
+             indi = 1
+             indj = indj + 1
+          ENDIF
+          ig(igrid) = indi
+          jg(igrid) = indj
+       ENDDO
+       ig(klon) = 1
+       jg(klon) = jjm + 1
+
+       IF ((.NOT. ALLOCATED(lalo))) THEN
+          ALLOCATE(lalo(knon,2), stat = error)
+          IF (error /= 0) THEN
+             abort_message='Pb allocation lalo'
+             CALL abort_gcm(modname,abort_message,1)
+          ENDIF
+       ENDIF
+       IF ((.NOT. ALLOCATED(lon_scat))) THEN
+          ALLOCATE(lon_scat(iim,jjm+1), stat = error)
+          IF (error /= 0) THEN
+             abort_message='Pb allocation lon_scat'
+             CALL abort_gcm(modname,abort_message,1)
+          ENDIF
+       ENDIF
+       IF ((.NOT. ALLOCATED(lat_scat))) THEN
+          ALLOCATE(lat_scat(iim,jjm+1), stat = error)
+          IF (error /= 0) THEN
+             abort_message='Pb allocation lat_scat'
+             CALL abort_gcm(modname,abort_message,1)
+          ENDIF
+       ENDIF
+       lon_scat = 0.
+       lat_scat = 0.
+       DO igrid = 1, knon
+          index = knindex(igrid)
+          lalo(igrid,2) = rlon(index)
+          lalo(igrid,1) = rlat(index)
+       ENDDO
+
+       
+       
+       CALL Gather(rlon,rlon_g)
+       CALL Gather(rlat,rlat_g)
+
+       IF (is_mpi_root) THEN
+          index = 1
+          DO jj = 2, jjm
+             DO ij = 1, iim
+                index = index + 1
+                lon_scat(ij,jj) = rlon_g(index)
+                lat_scat(ij,jj) = rlat_g(index)
+             ENDDO
+          ENDDO
+          lon_scat(:,1) = lon_scat(:,2)
+          lat_scat(:,1) = rlat_g(1)
+          lon_scat(:,jjm+1) = lon_scat(:,2)
+          lat_scat(:,jjm+1) = rlat_g(klon_glo)
+       ENDIF
+   
+
+!
+! Allouer et initialiser le tableau des voisins et des fraction de continents
+!
+       IF ( (.NOT.ALLOCATED(neighbours))) THEN
+          ALLOCATE(neighbours(knon,8), stat = error)
+          IF (error /= 0) THEN
+             abort_message='Pb allocation neighbours'
+             CALL abort_gcm(modname,abort_message,1)
+          ENDIF
+       ENDIF
+       neighbours = -1.
+       IF (( .NOT. ALLOCATED(contfrac))) THEN
+          ALLOCATE(contfrac(knon), stat = error)
+          IF (error /= 0) THEN
+             abort_message='Pb allocation contfrac'
+             CALL abort_gcm(modname,abort_message,1)
+          ENDIF
+       ENDIF
+
+       DO igrid = 1, knon
+          ireal = knindex(igrid)
+          contfrac(igrid) = pctsrf(ireal,is_ter)
+       ENDDO
+
+
+       CALL Init_neighbours(knon,neighbours,knindex,pctsrf(:,is_ter))
+
+!
+!  Allocation et calcul resolutions
+       IF ( (.NOT.ALLOCATED(resolution))) THEN
+          ALLOCATE(resolution(knon,2), stat = error)
+          IF (error /= 0) THEN
+             abort_message='Pb allocation resolution'
+             CALL abort_gcm(modname,abort_message,1)
+          ENDIF
+       ENDIF
+       DO igrid = 1, knon
+          ij = knindex(igrid)
+          resolution(igrid,1) = cuphy(ij)
+          resolution(igrid,2) = cvphy(ij)
+       ENDDO
+     
+       ALLOCATE(coastalflow(klon), stat = error)
+       IF (error /= 0) THEN
+          abort_message='Pb allocation coastalflow'
+          CALL abort_gcm(modname,abort_message,1)
+       ENDIF
+       
+       ALLOCATE(riverflow(klon), stat = error)
+       IF (error /= 0) THEN
+          abort_message='Pb allocation riverflow'
+          CALL abort_gcm(modname,abort_message,1)
+       ENDIF
+
+    ENDIF                          ! (fin debut) 
+
+! 
+! Appel a la routine sols continentaux
+!
+    IF (lafin) lrestart_write = .TRUE.
+    IF (check) WRITE(lunout,*)'lafin ',lafin,lrestart_write
+    
+    petA_orc(1:knon) = petBcoef(1:knon) * dtime
+    petB_orc(1:knon) = petAcoef(1:knon)
+    peqA_orc(1:knon) = peqBcoef(1:knon) * dtime
+    peqB_orc(1:knon) = peqAcoef(1:knon)
+
+    cdrag = 0.
+    cdrag(1:knon) = tq_cdrag(1:knon)
+
+! zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/287.05*temp_air(1:knon))*9.80665)
+    zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/RD*temp_air(1:knon))*RG)
+
+
+! PF et PASB
+!   where(cdrag > 0.01) 
+!     cdrag = 0.01
+!   endwhere
+!  write(*,*)'Cdrag = ',minval(cdrag),maxval(cdrag)
+
+!
+! Init Orchidee
+!
+!  if (pole_nord) then 
+!    offset=0
+!    ktindex(:)=ktindex(:)+iim-1
+!  else
+!    offset = klon_mpi_begin-1+iim-1
+!    ktindex(:)=ktindex(:)+MOD(offset,iim)
+!    offset=offset-MOD(offset,iim)
+!  endif
+  
+    IF (debut) THEN
+       CALL Get_orchidee_communicator(knon,orch_comm)
+       IF (knon /=0) THEN
+          CALL Init_orchidee_index(knon,orch_comm,knindex,offset,ktindex)
+          
+          CALL intersurf_main (itime+itau_phy-1, iim, jjm+1, offset, knon, ktindex, & 
+               orch_comm, dtime, lrestart_read, lrestart_write, lalo, &
+               contfrac, neighbours, resolution, date0, &
+               zlev,  u1_lay(1:knon), v1_lay(1:knon), spechum(1:knon), temp_air(1:knon), epot_air(1:knon), ccanopy(1:knon), &
+               cdrag(1:knon), petA_orc(1:knon), peqA_orc(1:knon), petB_orc(1:knon), peqB_orc(1:knon), &
+               precip_rain(1:knon), precip_snow(1:knon), lwdown(1:knon), swnet(1:knon), swdown(1:knon), ps(1:knon), &
+               evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
+               tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0_new(1:knon), &
+               lon_scat, lat_scat)
+          
+       ENDIF
+
+       albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
+
+    ENDIF
+
+!  swdown_vrai(1:knon) = swnet(1:knon)/(1. - albedo_keep(1:knon))
+    swdown_vrai(1:knon) = swdown(1:knon)
+
+    IF (knon /=0) THEN
+    
+       CALL intersurf_main (itime+itau_phy, iim, jjm+1,offset, knon, ktindex, & 
+            orch_comm,dtime, lrestart_read, lrestart_write, lalo, &
+            contfrac, neighbours, resolution, date0, &
+            zlev,  u1_lay(1:knon), v1_lay(1:knon), spechum(1:knon), temp_air(1:knon), epot_air(1:knon), ccanopy(1:knon), &
+            cdrag(1:knon), petA_orc(1:knon), peqA_orc(1:knon), petB_orc(1:knon), peqB_orc(1:knon), &
+            precip_rain(1:knon), precip_snow(1:knon), lwdown(1:knon), swnet(1:knon), swdown_vrai(1:knon), ps(1:knon), &
+            evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
+            tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0_new(1:knon), &
+            lon_scat, lat_scat)
+       
+    ENDIF
+
+    albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
+
+!* Send to coupler
+!
+    IF (ocean=='couple') THEN
+       CALL cpl_send_land_fields(itime, knon, knindex, &
+            riverflow, coastalflow)
+    ENDIF
+
+    alb_new(1:knon) = albedo_out(1:knon,1) 
+    alblw(1:knon) = albedo_out(1:knon,2)
+
+! Convention orchidee: positif vers le haut
+    fluxsens(1:knon) = -1. * fluxsens(1:knon)
+    fluxlat(1:knon)  = -1. * fluxlat(1:knon)
+    
+!  evap     = -1. * evap
+
+    IF (debut) lrestart_read = .FALSE.
+    
+  END SUBROUTINE surf_land_orchidee
+!
+!****************************************************************************************
+!
+  SUBROUTINE Init_orchidee_index(knon,orch_comm,knindex,offset,ktindex)
+    
+    INCLUDE "dimensions90.h"
+
+#ifdef CPP_PARA
+    INCLUDE 'mpif.h'
+#endif    
+
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                   :: knon
+    INTEGER, INTENT(IN)                   :: orch_comm
+    INTEGER, DIMENSION(klon), INTENT(IN)  :: knindex
+
+! Output arguments
+!****************************************************************************************
+    INTEGER, INTENT(OUT)                  :: offset
+    INTEGER, DIMENSION(knon), INTENT(OUT) :: ktindex
+
+! Local varables
+!****************************************************************************************
+#ifdef CPP_PARA
+    INTEGER, DIMENSION(MPI_STATUS_SIZE)   :: status
+#endif
+
+    INTEGER                               :: MyLastPoint
+    INTEGER                               :: LastPoint
+    INTEGER                               :: mpi_rank
+    INTEGER                               :: mpi_size
+    INTEGER                               :: ierr   
+    
+    MyLastPoint=klon_mpi_begin-1+knindex(knon)+iim-1
+    
+    IF (is_parallel) THEN
+#ifdef CPP_PARA    
+       CALL MPI_COMM_SIZE(orch_comm,mpi_size,ierr)
+       CALL MPI_COMM_RANK(orch_comm,mpi_rank,ierr)
+#endif
+    ELSE
+       mpi_rank=0
+       mpi_size=1
+    ENDIF
+    
+    IF (is_parallel) THEN
+       IF (.NOT. is_north_pole) THEN
+#ifdef CPP_PARA
+          CALL MPI_RECV(LastPoint,1,MPI_INTEGER,mpi_rank-1,1234,orch_comm,status,ierr)
+#endif
+       ENDIF
+       
+       IF (.NOT. is_south_pole) THEN
+#ifdef CPP_PARA
+          CALL MPI_SEND(MyLastPoint,1,MPI_INTEGER,mpi_rank+1,1234,orch_comm,ierr)  
+#endif
+       ENDIF
+    ENDIF
+    
+    IF (is_north_pole) THEN 
+       offset=0
+    ELSE
+       offset=LastPoint-MOD(LastPoint,iim)
+    ENDIF
+    
+    ktindex(1:knon)=knindex(1:knon)+(klon_mpi_begin+iim-1)-offset-1	
+    
+
+  END SUBROUTINE  Init_orchidee_index
+!
+!****************************************************************************************
+! 
+  SUBROUTINE Get_orchidee_communicator(knon,orch_comm)
+    
+#ifdef CPP_PARA
+    INCLUDE 'mpif.h'
+#endif    
+
+
+    INTEGER,INTENT(IN)  :: knon
+    INTEGER,INTENT(OUT) :: orch_comm
+    
+    INTEGER             :: color
+    INTEGER             :: ierr
+    
+    IF (knon==0) THEN 
+       color = 0
+    ELSE 
+       color = 1
+    ENDIF
+    
+#ifdef CPP_PARA    
+    CALL MPI_COMM_SPLIT(COMM_LMDZ_PHY,color,mpi_rank,orch_comm,ierr)
+#endif
+    
+  END SUBROUTINE Get_orchidee_communicator
+!
+!****************************************************************************************
+!  
+  SUBROUTINE Init_neighbours(knon,neighbours,ktindex,pctsrf)
+    
+    INCLUDE "indicesol.h"
+    INCLUDE "dimensions90.h"
+#ifdef CPP_PARA
+    INCLUDE 'mpif.h'
+#endif    
+
+
+    INTEGER, INTENT(IN)                     :: knon
+    INTEGER, DIMENSION(knon,8), INTENT(OUT) :: neighbours
+    INTEGER, DIMENSION(klon), INTENT(IN)    :: ktindex
+    REAL, DIMENSION(klon), INTENT(IN)       :: pctsrf
+    
+    INTEGER :: knon_nb(0:mpi_size-1)
+    INTEGER,DIMENSION(0:mpi_size-1) :: displs
+    INTEGER,ALLOCATABLE :: ktindex_g(:)
+    REAL*8  :: pctsrf_g(klon_glo)
+    INTEGER,ALLOCATABLE ::neighbours_g(:,:)
+    INTEGER :: knon_g
+    REAL*8 :: correspond(iim,jjm+1)
+    INTEGER :: i,igrid,jj,ij,iglob,ierr,ireal,index
+    INTEGER, DIMENSION(8,3) :: off_ini
+    INTEGER, DIMENSION(8)   :: offset  
+    INTEGER :: ktindex_p(knon)
+    
+    IF (is_sequential) THEN
+       knon_nb(:)=knon
+    ELSE  
+       
+#ifdef CPP_PARA  
+       CALL MPI_GATHER(knon,1,MPI_INTEGER,knon_nb,1,MPI_INTEGER,0,COMM_LMDZ_PHY,ierr)
+#endif
+       
+    ENDIF
+    
+    IF (is_mpi_root) THEN
+       knon_g=SUM(knon_nb(:))
+       ALLOCATE(ktindex_g(knon_g))
+       ALLOCATE(neighbours_g(knon_g,8))
+       neighbours_g(:,:)=-1
+       displs(0)=0
+       DO i=1,mpi_size-1
+          displs(i)=displs(i-1)+knon_nb(i-1)
+       ENDDO
+    ENDIF
+    
+    ktindex_p(1:knon)=ktindex(1:knon)+klon_mpi_begin-1+iim-1
+    
+    IF (is_sequential) THEN
+       ktindex_g(:)=ktindex_p(:)
+    ELSE
+       
+#ifdef CPP_PARA  
+       CALL MPI_GATHERV(ktindex_p,knon,MPI_INTEGER,ktindex_g,knon_nb,&
+            displs,MPI_INTEGER,0,COMM_LMDZ_PHY,ierr)
+#endif
+       
+    ENDIF
+    
+    CALL Gather(pctsrf,pctsrf_g)
+    
+    IF (is_mpi_root) THEN
+!  Initialisation des offset    
+!
+! offset bord ouest
+       off_ini(1,1) = - iim  ; off_ini(2,1) = - iim + 1; off_ini(3,1) = 1
+       off_ini(4,1) = iim + 1; off_ini(5,1) = iim      ; off_ini(6,1) = 2 * iim - 1
+       off_ini(7,1) = iim -1 ; off_ini(8,1) = - 1
+! offset point normal
+       off_ini(1,2) = - iim  ; off_ini(2,2) = - iim + 1; off_ini(3,2) = 1
+       off_ini(4,2) = iim + 1; off_ini(5,2) = iim      ; off_ini(6,2) = iim - 1
+       off_ini(7,2) = -1     ; off_ini(8,2) = - iim - 1
+! offset bord   est
+       off_ini(1,3) = - iim; off_ini(2,3) = - 2 * iim + 1; off_ini(3,3) = - iim + 1
+       off_ini(4,3) =  1   ; off_ini(5,3) = iim          ; off_ini(6,3) = iim - 1
+       off_ini(7,3) = -1   ; off_ini(8,3) = - iim - 1
+!
+!
+! Attention aux poles
+!
+       DO igrid = 1, knon_g
+          index = ktindex_g(igrid)
+          jj = INT((index - 1)/iim) + 1
+          ij = index - (jj - 1) * iim
+          correspond(ij,jj) = igrid
+       ENDDO
+       
+       DO igrid = 1, knon_g
+          iglob = ktindex_g(igrid)
+          IF (MOD(iglob, iim) == 1) THEN
+             offset = off_ini(:,1)
+          ELSE IF(MOD(iglob, iim) == 0) THEN
+             offset = off_ini(:,3)
+          ELSE
+             offset = off_ini(:,2)
+          ENDIF
+          DO i = 1, 8
+             index = iglob + offset(i)
+             ireal = (MIN(MAX(1, index - iim + 1), klon_glo))
+             IF (pctsrf_g(ireal) > EPSFRA) THEN
+                jj = INT((index - 1)/iim) + 1
+                ij = index - (jj - 1) * iim
+                neighbours_g(igrid, i) = correspond(ij, jj)
+             ENDIF
+          ENDDO
+       ENDDO
+
+    ENDIF
+    
+    DO i=1,8
+       IF (is_sequential) THEN
+          neighbours(:,i)=neighbours_g(:,i)
+       ELSE
+#ifdef CPP_PARA
+          CALL MPI_SCATTERV(neighbours_g(:,i),knon_nb,displs,MPI_INTEGER,neighbours(:,i),knon,MPI_INTEGER,0,COMM_LMDZ_PHY,ierr) 
+#endif
+       ENDIF
+    ENDDO
+    
+  END SUBROUTINE Init_neighbours
+!
+!****************************************************************************************
+!
+
+#endif
+
+END MODULE surf_land_orchidee_mod
Index: /LMDZ4/trunk/libf/phylmd/surf_landice_mod.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/surf_landice_mod.F90	(revision 781)
+++ /LMDZ4/trunk/libf/phylmd/surf_landice_mod.F90	(revision 781)
@@ -0,0 +1,170 @@
+!
+! $Header$
+!
+MODULE surf_landice_mod
+  
+  USE dimphy
+  USE surface_data,     ONLY : ocean, calice, calsno
+  USE fonte_neige_mod,  ONLY : fonte_neige, run_off_lic
+  USE cpl_mod,          ONLY : cpl_send_landice_fields
+  USE calcul_fluxs_mod, ONLY : calcul_fluxs
+
+  IMPLICIT NONE
+
+CONTAINS
+!
+!****************************************************************************************
+!
+  SUBROUTINE surf_landice(itime, dtime, knon, knindex, &
+       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
+       petAcoef, peqAcoef, petBcoef, peqBcoef, &
+       ps, u1_lay, v1_lay, rugoro, pctsrf, &
+       radsol, snow, qsurf, qsol, agesno, &
+       tsoil, z0_new, alblw, evap, fluxsens, fluxlat, &
+       tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_lic)
+
+    INCLUDE "indicesol.inc"
+    INCLUDE "dimsoil.h"
+    INCLUDE "YOMCST.inc"
+    INCLUDE "clesphys.inc"
+
+! Input variables 
+!****************************************************************************************
+    INTEGER, INTENT(IN)                           :: itime, knon
+    INTEGER, DIMENSION(klon), INTENT(in)          :: knindex
+    REAL, INTENT(in)                              :: dtime
+    REAL, DIMENSION(klon), INTENT(IN)             :: tsurf
+    REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
+    REAL, DIMENSION(klon), INTENT(IN)             :: tq_cdrag
+    REAL, DIMENSION(klon), INTENT(IN)             :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)             :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)             :: petAcoef, peqAcoef
+    REAL, DIMENSION(klon), INTENT(IN)             :: petBcoef, peqBcoef
+    REAL, DIMENSION(klon), INTENT(IN)             :: ps
+    REAL, DIMENSION(klon), INTENT(IN)             :: u1_lay, v1_lay
+    REAL, DIMENSION(klon), INTENT(IN)             :: rugoro
+    REAL, DIMENSION(klon,nbsrf), INTENT(IN)       :: pctsrf
+
+! In/Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
+    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+
+! Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
+    REAL, DIMENSION(klon), INTENT(OUT)            :: z0_new
+    REAL, DIMENSION(klon), INTENT(OUT)            :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
+    REAL, DIMENSION(klon), INTENT(OUT)            :: pctsrf_lic
+
+! Local variables
+!****************************************************************************************
+    REAL, DIMENSION(klon)    :: soilcap, soilflux
+    REAL, DIMENSION(klon)    :: cal, beta, dif_grnd
+    REAL, DIMENSION(klon)    :: zfra, alb_neig
+
+! End definition
+!****************************************************************************************
+!
+! Initialize output variables
+    alblw(:) =  999999.
+    alb_new(:) = 999999.
+
+
+!****************************************************************************************
+! Soil calculations
+! 
+!****************************************************************************************
+    IF (soil_model) THEN 
+       CALL soil(dtime, is_lic, knon, snow, tsurf, tsoil, soilcap, soilflux)
+       cal(1:knon) = RCPD / soilcap(1:knon)
+       radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
+    ELSE 
+       cal = RCPD * calice
+       WHERE (snow > 0.0) cal = RCPD * calsno
+    ENDIF
+
+
+!****************************************************************************************
+! Calulate fluxes
+!
+!****************************************************************************************
+    beta(:) = 1.0
+    dif_grnd(:) = 0.0
+
+    CALL calcul_fluxs(knon, is_lic, dtime, &
+         tsurf, p1lay, cal, beta, tq_cdrag, ps, &
+         precip_rain, precip_snow, snow, qsurf,  &
+         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+         petAcoef, peqAcoef, petBcoef, peqBcoef, &
+         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
+
+!****************************************************************************************
+! Calculate snow height, age, run-off,..
+!    
+!****************************************************************************************
+    CALL fonte_neige( knon, is_lic, knindex, dtime, &
+         tsurf, precip_rain, precip_snow, &
+         snow, qsol, tsurf_new, evap)
+
+
+!****************************************************************************************
+! Calculate albedo
+!
+!****************************************************************************************
+    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
+    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
+    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
+    alb_new(1:knon)  = alb_neig(1:knon)*zfra(1:knon) + &
+         0.6 * (1.0-zfra(1:knon))
+!
+!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
+!       alb_new(1 : knon)  = 0.6 !IM cf FH/GK 
+!       alb_new(1 : knon)  = 0.82
+!       alb_new(1 : knon)  = 0.77 !211003 Ksta0.77
+!       alb_new(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
+!IM: KstaTER0.77 & LMD_ARMIP6    
+
+! Attantion: alb_new and alblw are the same!
+    alb_new(1:knon)  = 0.77
+    alblw(1:knon)    = alb_new(1:knon)
+
+
+!****************************************************************************************
+! Rugosity
+!
+!****************************************************************************************
+    z0_new(:) = rugoro(:)
+
+
+!****************************************************************************************
+! Return the pourcentage for this sub-surface
+!
+!****************************************************************************************
+    pctsrf_lic(:) = pctsrf(:,is_lic)
+    
+
+!****************************************************************************************
+! Send run-off on land-ice to coupler if coupled ocean.
+! run_off_lic has been calculated in fonte_neige
+!
+!****************************************************************************************
+    IF (ocean=='couple') THEN
+       CALL cpl_send_landice_fields(itime, knon, knindex, run_off_lic)
+    ENDIF
+
+
+  END SUBROUTINE surf_landice
+!
+!****************************************************************************************
+!
+END MODULE surf_landice_mod
+
+
+
Index: /LMDZ4/trunk/libf/phylmd/surf_ocean_mod.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/surf_ocean_mod.F90	(revision 781)
+++ /LMDZ4/trunk/libf/phylmd/surf_ocean_mod.F90	(revision 781)
@@ -0,0 +1,164 @@
+!
+! $Header$
+!
+MODULE surf_ocean_mod
+
+  USE dimphy
+  USE surface_data, ONLY     : ocean
+  USE ocean_forced_mod, ONLY : ocean_forced_noice
+  USE ocean_slab_mod, ONLY   : ocean_slab_noice
+  USE ocean_cpl_mod, ONLY    : ocean_cpl_noice
+
+  IMPLICIT NONE
+
+CONTAINS
+!
+!****************************************************************************************
+!
+  SUBROUTINE surf_ocean(rlon, rlat, sollw, albedo, &
+       rugos, windsp, rmu0, &
+       fder, &
+       itime, dtime, jour, knon, knindex, &
+       debut, swdown, &
+       p1lay, tq_cdrag, coefm, precip_rain, precip_snow, temp_air, spechum, &
+       petAcoef, peqAcoef, petBcoef, peqBcoef, &
+       ps, u1_lay, v1_lay, rugoro, pctsrf, &
+       radsol, snow, qsurf, agesno, &
+       z0_new, alblw, evap, fluxsens, fluxlat, &
+       tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_oce)
+!
+! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force, 
+! slab or couple). The calculations of albedo and rugosity for the ocean surface are 
+! done in here because they are identical for the different modes of ocean. 
+!
+    INCLUDE "indicesol.inc"
+    INCLUDE "YOMCST.inc"
+
+! Input variables
+!****************************************************************************************
+    INTEGER, INTENT(IN)                      :: itime, jour, knon
+    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
+    REAL, INTENT(IN)                         :: dtime
+    REAL, DIMENSION(klon), INTENT(IN)        :: rlon, rlat
+    REAL, DIMENSION(klon), INTENT(IN)        :: sollw
+    REAL, DIMENSION(klon), INTENT(IN)        :: albedo
+    REAL, DIMENSION(klon), INTENT(IN)        :: rugos
+    REAL, DIMENSION(klon), INTENT(IN)        :: windsp
+    REAL, DIMENSION(klon), INTENT(IN)        :: rmu0  
+    REAL, DIMENSION(klon), INTENT(IN)        :: fder
+    REAL, DIMENSION(klon), INTENT(IN)        :: swdown    
+    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
+    REAL, DIMENSION(klon), INTENT(IN)        :: tq_cdrag
+    REAL, DIMENSION(klon), INTENT(IN)        :: coefm
+    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)        :: petAcoef, peqAcoef
+    REAL, DIMENSION(klon), INTENT(IN)        :: petBcoef, peqBcoef
+    REAL, DIMENSION(klon), INTENT(IN)        :: ps
+    REAL, DIMENSION(klon), INTENT(IN)        :: u1_lay, v1_lay
+    REAL, DIMENSION(klon), INTENT(IN)        :: rugoro
+    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
+    LOGICAL, INTENT(IN)                      :: debut
+
+! In/Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT)     :: radsol
+    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow, qsurf
+    REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno
+
+! Output variables
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: z0_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
+    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_oce
+
+
+! Local variables
+!****************************************************************************************
+    INTEGER               :: i
+    REAL, DIMENSION(klon) :: alb_eau
+
+! End definition
+!****************************************************************************************
+! 
+! Switch according to type of ocean (couple, slab or forced)
+!****************************************************************************************
+    SELECT CASE(ocean)
+    CASE('couple')
+       CALL ocean_cpl_noice(&
+            sollw, albedo, &
+            windsp, &
+            fder, & 
+            itime, dtime, knon, knindex, &
+            swdown, &
+            p1lay, tq_cdrag, precip_rain, precip_snow,temp_air,spechum,&
+            petAcoef, peqAcoef, petBcoef, peqBcoef, &
+            ps, u1_lay, v1_lay, pctsrf, &
+            radsol, snow, qsurf, agesno, &
+            evap, fluxsens, fluxlat, &
+            tsurf_new, dflux_s, dflux_l, pctsrf_oce)
+
+    CASE('slab')
+       CALL ocean_slab_noice(&
+            dtime, knon, knindex, &
+            p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum,&
+            petAcoef, peqAcoef, petBcoef, peqBcoef, &
+            ps, u1_lay, v1_lay, &
+            radsol, snow, qsurf, agesno, &
+            evap, fluxsens, fluxlat, &
+            tsurf_new, &
+            dflux_s, dflux_l, pctsrf_oce)
+       
+    CASE('force')
+       CALL ocean_forced_noice(itime, dtime, jour, knon, knindex, &
+            debut, &
+            p1lay, tq_cdrag, precip_rain, precip_snow, &
+            temp_air, spechum, &
+            petAcoef, peqAcoef, petBcoef, peqBcoef, &
+            ps, u1_lay, v1_lay, &
+            radsol, snow, qsurf, &
+            agesno, &
+            evap, fluxsens, fluxlat, &
+            tsurf_new, dflux_s, dflux_l, pctsrf_oce)
+    END SELECT
+
+!****************************************************************************************
+! Calculate albedo
+!
+!****************************************************************************************
+    IF ( MINVAL(rmu0) == MAXVAL(rmu0) .AND. MINVAL(rmu0) == -999.999 ) THEN
+       CALL alboc(FLOAT(jour),rlat,alb_eau)
+    ELSE  ! diurnal cycle
+       CALL alboc_cd(rmu0,alb_eau)
+    ENDIF
+
+    DO i =1, knon
+       alb_new(i) = alb_eau(knindex(i))
+    ENDDO
+
+!****************************************************************************************
+! Calculate the rugosity
+!
+!****************************************************************************************
+    z0_new = SQRT(rugos**2 + rugoro**2)
+
+    ! The rugosity is recalculated with another method
+    z0_new(:) = 0.0
+    DO i = 1, knon
+       z0_new(i) = 0.018*coefm(i) * (u1_lay(i)**2+v1_lay(i)**2)/RG  &
+            +  0.11*14e-6 / SQRT(coefm(i) * (u1_lay(i)**2+v1_lay(i)**2))
+       z0_new(i) = MAX(1.5e-05,z0_new(i))
+    ENDDO
+    
+    alblw(1:knon) = alb_new(1:knon)
+!
+!****************************************************************************************
+!    
+  END SUBROUTINE surf_ocean
+!
+!****************************************************************************************
+!
+END MODULE surf_ocean_mod
Index: /LMDZ4/trunk/libf/phylmd/surf_seaice_mod.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/surf_seaice_mod.F90	(revision 781)
+++ /LMDZ4/trunk/libf/phylmd/surf_seaice_mod.F90	(revision 781)
@@ -0,0 +1,133 @@
+!
+! $Header$
+!
+MODULE surf_seaice_mod
+
+  USE dimphy
+  USE surface_data, ONLY     : ocean
+  USE ocean_slab_mod, ONLY   : ocean_slab_ice
+  USE ocean_forced_mod, ONLY : ocean_forced_ice
+  USE ocean_cpl_mod, ONLY    : ocean_cpl_ice
+  IMPLICIT NONE
+
+CONTAINS
+!
+!****************************************************************************************
+!
+  SUBROUTINE surf_seaice( & 
+       rlon, rlat, sollw, albedo, &
+       fder, &
+       itime, dtime, jour, knon, knindex, &
+       debut, lafin, swdown, &
+       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
+       petAcoef, peqAcoef, petBcoef, peqBcoef, &
+       ps, u1_lay, v1_lay, rugoro, pctsrf, &
+       radsol, snow, qsurf, qsol, agesno, &
+       tsoil, z0_new, alblw, evap, fluxsens, fluxlat, &
+       tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_sic)
+!
+! This subroutine will make a call to ocean_XXX_ice according to the ocean mode (force, 
+! slab or couple). The calculation of rugosity for the sea-ice surface is also done
+! in here because it is the same calculation for the different modes of ocean.
+!
+    INCLUDE "indicesol.inc"
+    INCLUDE "dimsoil.h"
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                      :: itime, jour, knon
+    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
+    LOGICAL, INTENT(IN)                      :: debut, lafin
+    REAL, INTENT(IN)                         :: dtime
+    REAL, DIMENSION(klon), INTENT(IN)        :: rlon, rlat
+    REAL, DIMENSION(klon), INTENT(IN)        :: sollw
+    REAL, DIMENSION(klon), INTENT(IN)        :: albedo
+    REAL, DIMENSION(klon), INTENT(IN)        :: fder
+    REAL, DIMENSION(klon), INTENT(IN)        :: swdown    
+    REAL, DIMENSION(klon), INTENT(IN)        :: tsurf
+    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
+    REAL, DIMENSION(klon), INTENT(IN)        :: tq_cdrag
+    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)        :: petAcoef, peqAcoef
+    REAL, DIMENSION(klon), INTENT(IN)        :: petBcoef, peqBcoef
+    REAL, DIMENSION(klon), INTENT(IN)        :: ps
+    REAL, DIMENSION(klon), INTENT(IN)        :: u1_lay, v1_lay
+    REAL, DIMENSION(klon), INTENT(IN)        :: rugoro
+    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
+
+! In/Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsurf, qsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
+    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
+
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: z0_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
+    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_sic
+
+! End definitions
+!****************************************************************************************
+
+!****************************************************************************************
+! Switch according to type of ocean (couple, slab or forced)
+!
+!****************************************************************************************
+    SELECT CASE(ocean)
+    CASE('couple')
+       CALL ocean_cpl_ice( &
+          rlon, rlat, sollw, albedo, & 
+          fder, & 
+          itime, dtime, knon, knindex, &
+          lafin,&
+          swdown, &
+          p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum,&
+          petAcoef, peqAcoef, petBcoef, peqBcoef, &
+          ps, u1_lay, v1_lay, pctsrf, &
+          radsol, snow, qsurf, &
+          alblw, evap, fluxsens, fluxlat, &
+          tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_sic)
+
+    CASE('slab')
+       CALL ocean_slab_ice( & 
+          itime, dtime, jour, knon, knindex, &
+          debut, &
+          tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum,&
+          petAcoef, peqAcoef, petBcoef, peqBcoef, &
+          ps, u1_lay, v1_lay, &
+          radsol, snow, qsurf, qsol, agesno, &
+          tsoil, &
+          alblw, evap, fluxsens, fluxlat, &
+          tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_sic)
+    
+    CASE('force')
+       CALL ocean_forced_ice(itime, dtime, jour, knon, knindex, &
+          debut, &
+          tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum,&
+          petAcoef, peqAcoef, petBcoef, peqBcoef, &
+          ps, u1_lay, v1_lay, &
+          radsol, snow, qsurf, qsol, agesno, &
+          tsoil, alblw, evap, fluxsens, fluxlat, &
+          tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_sic)
+    END SELECT
+
+!****************************************************************************************
+! Calculate rugosity
+!
+!****************************************************************************************
+    z0_new = 0.002
+    z0_new = SQRT(z0_new**2+rugoro**2)
+
+
+  END SUBROUTINE surf_seaice
+!
+!****************************************************************************************
+!
+END MODULE surf_seaice_mod
+
Index: /LMDZ4/trunk/libf/phylmd/surface_data.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/surface_data.F90	(revision 781)
+++ /LMDZ4/trunk/libf/phylmd/surface_data.F90	(revision 781)
@@ -0,0 +1,22 @@
+!
+! $Header$
+!
+MODULE surface_data
+
+  USE dimphy, ONLY : klon
+
+  REAL, PARAMETER        :: calice=1.0/(5.1444e+06*0.15)
+  REAL, PARAMETER        :: tau_gl=86400.*5.
+  REAL, PARAMETER        :: calsno=1./(2.3867e+06*.15)
+  
+  LOGICAL, SAVE          :: debug_surf=.FALSE.
+  !$OMP THREADPRIVATE(debug_surf)
+  LOGICAL, SAVE          :: newwind=.FALSE.
+  !$OMP THREADPRIVATE(newwind)
+  
+  LOGICAL, SAVE          :: ok_veget       ! true for use of vegetation model ORCHIDEE
+  !$OMP THREADPRIVATE(ok_veget)
+  CHARACTER(len=6), SAVE :: ocean          ! force/slab/couple
+  !$OMP THREADPRIVATE(ocean)
+
+END MODULE surface_data
