Index: LMDZ4/trunk/libf/phylmd/ocean_cpl_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/ocean_cpl_mod.F90	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/ocean_cpl_mod.F90	(revision 888)
@@ -65,14 +65,13 @@
 !
   SUBROUTINE ocean_cpl_noice( &
-       sollw, albedo, &
+       swnet, lwnet, alb1, &
        windsp, &
        fder_old, &
        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_in, &
-       radsol, snow, qsurf, agesno, &
-       evap, fluxsens, fluxlat, &
+       radsol, snow, agesno, &
+       qsurf, evap, fluxsens, fluxlat, &
        tsurf_new, dflux_s, dflux_l, pctsrf_oce)
 !
@@ -89,9 +88,9 @@
     INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
     REAL, INTENT(IN)                         :: dtime
-    REAL, DIMENSION(klon), INTENT(IN)        :: sollw
-    REAL, DIMENSION(klon), INTENT(IN)        :: albedo
+    REAL, DIMENSION(klon), INTENT(IN)        :: swnet
+    REAL, DIMENSION(klon), INTENT(IN)        :: lwnet
+    REAL, DIMENSION(klon), INTENT(IN)        :: alb1   ! albedo in visible SW interval
     REAL, DIMENSION(klon), INTENT(IN)        :: windsp
     REAL, DIMENSION(klon), INTENT(IN)        :: fder_old
-    REAL, DIMENSION(klon), INTENT(IN)        :: swdown    
     REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
     REAL, DIMENSION(klon), INTENT(IN)        :: tq_cdrag
@@ -199,6 +198,6 @@
 
     CALL cpl_send_ocean_fields(itime, knon, knindex, &
-         swdown, sollw, fluxlat, fluxsens, &
-         precip_rain, precip_snow, evap, tsurf_new, fder_new, albedo, taux, tauy, windsp)
+         swnet, lwnet, fluxlat, fluxsens, &
+         precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, taux, tauy, windsp)
     
 
@@ -208,15 +207,14 @@
 !
   SUBROUTINE ocean_cpl_ice( &
-       rlon, rlat, sollw, albedo, &
+       rlon, rlat, swnet, lwnet, alb1, &
        fder_old, &
        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_in, &
        radsol, snow, qsurf, &
-       alblw, evap, fluxsens, fluxlat, &
-       tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_sic)
+       alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+       tsurf_new, dflux_s, dflux_l, pctsrf_sic)
 !
 ! This subroutine treats the ocean where there is ice. The subroutine first receives 
@@ -234,8 +232,8 @@
     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)        :: swnet
+    REAL, DIMENSION(klon), INTENT(IN)        :: lwnet
+    REAL, DIMENSION(klon), INTENT(IN)        :: alb1   ! albedo in visible SW interval
     REAL, DIMENSION(klon), INTENT(IN)        :: fder_old
-    REAL, DIMENSION(klon), INTENT(IN)        :: swdown
     REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
     REAL, DIMENSION(klon), INTENT(IN)        :: tq_cdrag
@@ -256,7 +254,7 @@
 !****************************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
-    REAL, DIMENSION(klon), INTENT(OUT)       :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new, alb2_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
-    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
     REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_sic
@@ -271,4 +269,5 @@
     REAL, DIMENSION(klon)   :: tsurf_cpl, fder_new
     REAL, DIMENSION(klon)   :: taux, tauy
+    REAL, DIMENSION(klon)   :: alb_cpl
 
 ! End definitions
@@ -278,5 +277,5 @@
 
 !****************************************************************************************
-! Receive ocean temperature(tsurf_cpl), albedo(alb_new) and new fraction of 
+! Receive ocean temperature(tsurf_cpl), albedo(alb_cpl) and new fraction of 
 ! seaice(pctsrf_sic) from coupler
 !
@@ -284,5 +283,9 @@
 
     CALL cpl_receive_seaice_fields(knon, knindex, &
-         tsurf_cpl, alb_new, pctsrf_sic)
+         tsurf_cpl, alb_cpl, pctsrf_sic)
+
+    alb1_new(1:knon) = alb_cpl(1:knon)
+    alb2_new(1:knon) = alb_cpl(1:knon)    
+
     
 !****************************************************************************************
@@ -342,9 +345,7 @@
     CALL cpl_send_seaice_fields(itime, dtime, knon, knindex, &
        pctsrf_in, lafin, rlon, rlat, &
-       swdown, sollw, fluxlat, fluxsens, &
-       precip_rain, precip_snow, evap, tsurf_new, fder_new, albedo, taux, tauy)
+       swnet, lwnet, fluxlat, fluxsens, &
+       precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, taux, tauy)
  
-
-    alblw(1:knon) = alb_new(1:knon)    
 
   END SUBROUTINE ocean_cpl_ice
Index: LMDZ4/trunk/libf/phylmd/ocean_forced_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/ocean_forced_mod.F90	(revision 888)
@@ -69,7 +69,6 @@
        petAcoef, peqAcoef, petBcoef, peqBcoef, &
        ps, u1_lay, v1_lay, &
-       radsol, snow, qsurf, &
-       agesno, &
-       evap, fluxsens, fluxlat, &
+       radsol, snow, agesno, & 
+       qsurf, evap, fluxsens, fluxlat, &
        tsurf_new, dflux_s, dflux_l, pctsrf_oce)
 !
@@ -191,7 +190,7 @@
        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)
+       radsol, snow, qsol, agesno, tsoil, &
+       qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+       tsurf_new, dflux_s, dflux_l, pctsrf_sic)
 !
 ! This subroutine treats the ocean where there is ice. 
@@ -232,7 +231,8 @@
 !****************************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
-    REAL, DIMENSION(klon), INTENT(OUT)            :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
+    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
     REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
-    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
     REAL, DIMENSION(klon), INTENT(OUT)            :: pctsrf_sic
@@ -286,5 +286,4 @@
     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)
@@ -320,10 +319,11 @@
     WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.
 
-    alb_new(:) = 0.0
+    alb1_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)
+       alb1_new(i) = alb_neig(i) * zfra +  0.6 * (1.0-zfra)
     ENDDO
-!!      alb_new(1 : knon) = 0.6
+
+    alb2_new(:) = alb1_new(:)
 
 !****************************************************************************************
@@ -344,7 +344,4 @@
     ENDDO
 
-!!$    z0_new = 0.002
-!!$    z0_new = SQRT(z0_new**2+rugoro**2)
-    alblw(1:knon) = alb_new(1:knon)
 
 !****************************************************************************************
@@ -355,5 +352,4 @@
 
     pctsrf_sic(:) = pctsrf_lim(:,is_sic)
-
 
   END SUBROUTINE ocean_forced_ice
Index: LMDZ4/trunk/libf/phylmd/ocean_slab_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 888)
@@ -139,8 +139,7 @@
        petAcoef, peqAcoef, petBcoef, peqBcoef, &
        ps, u1_lay, v1_lay, &
-       radsol, snow, qsurf, agesno, &
-       evap, fluxsens, fluxlat, &
-       tsurf_new, &
-       dflux_s, dflux_l, pctsrf_oce)
+       radsol, snow, agesno, &
+       qsurf, evap, fluxsens, fluxlat, &
+       tsurf_new, dflux_s, dflux_l, pctsrf_oce)
 
     INCLUDE "indicesol.h"
@@ -179,5 +178,5 @@
     INTEGER                :: i
     REAL, DIMENSION(klon)  :: cal, beta, dif_grnd
-    REAL, DIMENSION(klon)  :: alb_neig, tsurf_temp
+    REAL, DIMENSION(klon)  :: tsurf_temp
 
 !****************************************************************************************
@@ -191,5 +190,4 @@
     beta = 1.
     dif_grnd = 0.
-    alb_neig(:) = 0.
     agesno(:) = 0.
     
@@ -221,8 +219,7 @@
        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)
+       radsol, snow, qsurf, qsol, agesno, tsoil, &
+       alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+       tsurf_new, dflux_s, dflux_l, pctsrf_sic)
 
     INCLUDE "indicesol.h"
@@ -258,7 +255,8 @@
 !****************************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
-    REAL, DIMENSION(klon), INTENT(OUT)            :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
+    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
     REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
-    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
     REAL, DIMENSION(klon), INTENT(OUT)            :: pctsrf_sic
@@ -282,5 +280,5 @@
 
 ! Initialization of output variables
-    alb_new(:) = 0.0
+    alb1_new(:) = 0.0
 
 !****************************************************************************************
@@ -347,7 +345,8 @@
     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) + & 
+    alb1_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + & 
          0.6 * (1.0-zfra(1:knon))
-    !!      alb_new(1 : knon) = 0.6
+    
+    alb2_new(:) = alb1_new(:)
 
 !
@@ -397,7 +396,4 @@
     ENDDO
      
-!!$    z0_new = 0.002
-!!$    z0_new = SQRT(z0_new**2+rugoro**2)
-    alblw(1:knon) = alb_new(1:knon)
         
 !****************************************************************************************
Index: LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 888)
@@ -26,19 +26,19 @@
 
 ! Declaration of variables saved in restart file
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: qsol
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: qsol   ! water height in the soil (mm)
   !$OMP THREADPRIVATE(qsol)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder   ! flux drift
   !$OMP THREADPRIVATE(fder)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: snow
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: snow   ! snow at surface
   !$OMP THREADPRIVATE(snow)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf  ! humidity at surface
   !$OMP THREADPRIVATE(qsurf)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: evap
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: evap   ! evaporation at surface
   !$OMP THREADPRIVATE(evap)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: rugos
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: rugos  ! rugosity at surface (m)
   !$OMP THREADPRIVATE(rugos)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: agesno
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: agesno ! age of snow at surface
   !$OMP THREADPRIVATE(agesno)
-  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: ftsoil  
+  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: ftsoil ! soil temperature
   !$OMP THREADPRIVATE(ftsoil)
 
@@ -182,7 +182,7 @@
        t,         q,         u,        v,             &
        pplay,     paprs,     pctsrf,                  &
-       ts,        albe,      alblw,    u10m,   v10m,  &
-       sollwdown, cdragh,    cdragm,   zu1,    zv1,   &
-       albsol,    albsollw,  zxsens,   zxevap,        &
+       ts,        alb1,      alb2,     u10m,   v10m,  &
+       lwdown_m,  cdragh,    cdragm,   zu1,    zv1,   &
+       alb1_m,    alb2_m,    zxsens,   zxevap,        &
        zxtsol,    zxfluxlat, zt2m,     qsat2m,        &
        d_t,       d_q,       d_u,      d_v,           & 
@@ -270,83 +270,101 @@
 ! 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
+    REAL,                         INTENT(IN)        :: dtime   ! time interval (s)
+    REAL,                         INTENT(IN)        :: date0   ! initial day
+    INTEGER,                      INTENT(IN)        :: itap    ! time step
+    INTEGER,                      INTENT(IN)        :: jour    ! current day of the year
+    LOGICAL,                      INTENT(IN)        :: debut   ! true if first run step
+    LOGICAL,                      INTENT(IN)        :: lafin   ! true if last run step
+    REAL, DIMENSION(klon),        INTENT(IN)        :: rlon    ! longitudes in degrees
+    REAL, DIMENSION(klon),        INTENT(IN)        :: rlat    ! latitudes in degrees
+    REAL, DIMENSION(klon),        INTENT(IN)        :: rugoro  ! rugosity length
+    REAL, DIMENSION(klon),        INTENT(IN)        :: rmu0    ! cosine of solar zenith angle
+    REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f  ! rain fall
+    REAL, DIMENSION(klon),        INTENT(IN)        :: snow_f  ! snow fall
+    REAL, DIMENSION(klon),        INTENT(IN)        :: solsw_m ! net shortwave radiation at mean surface
+    REAL, DIMENSION(klon),        INTENT(IN)        :: sollw_m ! net longwave radiation at mean surface
+    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t       ! temperature (K)
+    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q       ! water vapour (kg/kg)
+    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: u       ! u speed
+    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: v       ! v speed
+    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pplay   ! mid-layer pression (Pa)
+    REAL, DIMENSION(klon,klev+1), INTENT(IN)        :: paprs   ! pression between layers (Pa) 
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN)        :: pctsrf  ! sub-surface fraction
 
 ! 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
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts      ! temperature at surface (K)
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: alb1    ! albedo in visible SW interval
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: alb2    ! albedo in near infra-red SW interval
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: u10m    ! u speed at 10m
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: v10m    ! v speed at 10m
 
 ! 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)       :: lwdown_m   ! Downcoming longwave radiation
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh     ! drag coefficient for T and Q
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm     ! drag coefficient for wind
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu1        ! u wind speed in first layer
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv1        ! v wind speed in first layer
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: alb1_m     ! mean albedo in visible SW interval
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: alb2_m     ! mean albedo in near IR SW interval
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens     ! sensible heat flux at surface with inversed sign 
+                                                                  ! (=> positive sign upwards)
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxevap     ! water vapour flux at surface, positiv upwards
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxtsol     ! temperature at surface, mean for each grid point
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat  ! latent flux, mean for each grid point
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zt2m       ! temperature at 2m, mean for each grid point
     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
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t        ! change in temperature 
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_q        ! change in water vapour
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_u        ! change in u speed
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_v        ! change in v speed
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zcoefh     ! coef for turbulent diffusion of T and Q, mean for each grid point
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: pctsrf_new ! new sub-surface fraction
 
 ! 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
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsol_d     ! water height in the soil (mm)
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zq2m       ! water vapour at 2m, mean for each grid point
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh     ! height of the planetary boundary layer(HPBL)
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl     ! condensation level
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_capCL    ! CAPE of PBL
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_oliqCL   ! liquid water intergral of PBL
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_cteiCL   ! cloud top instab. crit. of PBL
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblT     ! temperature at PBLH
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_therm    ! thermal virtual temperature excess
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb1    ! deep cape, mean for each grid point
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb2    ! inhibition, mean for each grid point
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb3    ! point Omega, mean for each grid point
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxrugs     ! rugosity at surface (m), mean for each grid point
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu10m      ! u speed at 10m, mean for each grid point
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv10m      ! v speed at 10m, mean for each grid point
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxqsurf    ! humidity at surface, mean for each grid point
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: rh2m       ! relative humidity at 2m
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxu    ! u wind tension, mean for each grid point
+    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxv    ! v wind tension, mean for each grid point
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: rugos_d    ! rugosity length (m)
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: agesno_d   ! age of snow at surface
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: solsw      ! net shortwave radiation at surface 
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: sollw      ! net longwave radiation at surface
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: d_ts       ! change in temperature at surface
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: evap_d     ! evaporation at surface
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: fluxlat    ! latent flux
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m        ! temperature at 2 meter height
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfbils     ! heat balance at surface
+    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfbilo     ! water balance at surface
+    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t     ! sensible heat flux (CpT) J/m**2/s (W/m**2)
+                                                                  ! positve orientation downwards
+    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u     ! u wind tension (kg m/s)/(m**2 s) or Pascal
+    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v     ! v wind tension (kg m/s)/(m**2 s) or Pascal
 
 ! 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
+    REAL, DIMENSION(klon),       INTENT(OUT)        :: dflux_t    ! change of sensible heat flux 
+    REAL, DIMENSION(klon),       INTENT(OUT)        :: dflux_q    ! change of water vapour flux
+    REAL, DIMENSION(klon),       INTENT(OUT)        :: zxsnow     ! snow at surface, mean for each grid point
+    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxt    ! sensible heat flux, mean for each grid point
+    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxq    ! water vapour flux, mean for each grid point
+    REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m        ! water vapour at 2 meter height
+    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q     ! water vapour flux(latent flux) (kg/m**2/s)
 
 ! Input/output
@@ -356,5 +374,5 @@
 ! Local variables with attribute SAVE
 !****************************************************************************************
-    INTEGER, SAVE                            :: nhoridbg, nidbg
+    INTEGER, SAVE                            :: nhoridbg, nidbg   ! variables for IOIPSL
 !$OMP THREADPRIVATE(nhoridbg, nidbg)
     LOGICAL, SAVE                            :: debugindex=.FALSE.
@@ -373,14 +391,14 @@
     REAL                               :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
     REAL                               :: amn, amx
+    REAL                               :: f1 ! fraction de longeurs visibles parmi tout SW intervalle
     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)              :: yalb, yalb1, yalb2
     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)              :: ysolsw, ysollw
     REAL, DIMENSION(klon)              :: yfder
-    REAL, DIMENSION(klon)              :: yrads,yrugoro
+    REAL, DIMENSION(klon)              :: yrugoro
     REAL, DIMENSION(klon)              :: yfluxlat
     REAL, DIMENSION(klon)              :: y_d_ts
@@ -409,10 +427,11 @@
     REAL, DIMENSION(klon)              :: qairsol, zgeo1
     REAL, DIMENSION(klon)              :: rugo1
-    REAL, DIMENSION(klon)              :: yfluxsens, swdown
+    REAL, DIMENSION(klon)              :: yfluxsens
     REAL, DIMENSION(klon)              :: petAcoef, peqAcoef, petBcoef, peqBcoef
-    REAL, DIMENSION(klon)              :: ypsref, epot_air
-    REAL, DIMENSION(klon)              :: yevap, ytsurf_new, yalb_new
+    REAL, DIMENSION(klon)              :: ypsref
+    REAL, DIMENSION(klon)              :: yevap, ytsurf_new, yalb1_new, yalb2_new
     REAL, DIMENSION(klon)              :: pctsrf_nsrf
     REAL, DIMENSION(klon)              :: ztsol
+    REAL, DIMENSION(klon)              :: alb_m  ! mean albedo for whole SW interval
     REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q
     REAL, DIMENSION(klon,klev)         :: y_d_u, y_d_v
@@ -441,6 +460,6 @@
 
 
-    REAL, DIMENSION(klon,nbsrf)        :: pblh
-    REAL, DIMENSION(klon,nbsrf)        :: plcl
+    REAL, DIMENSION(klon,nbsrf)        :: pblh         ! height of the planetary boundary layer
+    REAL, DIMENSION(klon,nbsrf)        :: plcl         ! condensation level
     REAL, DIMENSION(klon,nbsrf)        :: capCL
     REAL, DIMENSION(klon,nbsrf)        :: oliqCL
@@ -448,10 +467,12 @@
     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)        :: trmb1        ! deep cape
+    REAL, DIMENSION(klon,nbsrf)        :: trmb2        ! inhibition
+    REAL, DIMENSION(klon,nbsrf)        :: trmb3        ! point Omega
     REAL, DIMENSION(klon,nbsrf)        :: zx_rh2m, zx_qsat2m
     REAL, DIMENSION(klon,nbsrf)        :: zx_qs1, zx_t1
     REAL, DIMENSION(klon,nbsrf)        :: zdelta1, zcor1
+    REAL, DIMENSION(klon, nbsrf)       :: alb          ! mean albedo for whole SW interval
+    REAL, DIMENSION(klon)              :: ylwdown      ! jg : temporary (ysollwdown)
 
 
@@ -521,8 +542,8 @@
     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    
+    zv1 = 0.0     ; yqsurf = 0.0     ; yalb1 = 0.0     ; yalb2 = 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
+    ysollw = 0.0  ; yrugos = 0.0     ; yu1 = 0.0    
+    yv1 = 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 
@@ -538,5 +559,6 @@
 ! 3) - Calculate pressure thickness of each layer
 !    - Calculate the wind at first layer
-!
+!    - Mean calculations of albedo
+!    - Calculate net radiance at sub-surface
 !****************************************************************************************
     DO k = 1, klev
@@ -555,5 +577,5 @@
 !****************************************************************************************
 ! Test for rugos........ from physiq.. A la fin plutot???
-! Calcul de l'abedo moyen par maille
+!
 !****************************************************************************************
 
@@ -566,17 +588,41 @@
     ENDDO
 
-! Calcul de l'abedo moyen par maille
-    albsol(:)   = 0.0
-    albsollw(:) = 0.0
+! Mean calculations of albedo
+!
+! Albedo at sub-surface
+! * alb1 : albedo in visible SW interval
+! * alb2 : albedo in near infrared SW interval
+! * alb  : mean albedo for whole SW interval
+!
+! Mean albedo for grid point
+! * alb1_m : albedo in visible SW interval
+! * alb2_m : albedo in near infrared SW interval
+! * alb_m  : mean albedo at whole SW interval
+
+    alb1_m(:) = 0.0
+    alb2_m(:) = 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)
+          alb1_m(i) = alb1_m(i) + alb1(i,nsrf) * pctsrf(i,nsrf)
+          alb2_m(i) = alb2_m(i) + alb2(i,nsrf) * pctsrf(i,nsrf)
        ENDDO
     ENDDO
 
-
-
-! Calcule de ztsol (aussi fait dans physiq.F, pourrait etre un argument)
+! We here suppose the fraction f1 of incoming radiance of visible radiance 
+! as a fraction of all shortwave radiance 
+!    f1 = 0.5 
+    f1 = 1    ! put f1=1 to recreate old calculations
+
+    DO nsrf = 1, nbsrf
+       DO i = 1, klon
+          alb(i,nsrf) = f1*alb1(i,nsrf) + (1-f1)*alb2(i,nsrf)
+       ENDDO
+    ENDDO
+
+    DO i = 1, klon
+       alb_m(i) = f1*alb1_m(i) + (1-f1)*alb2_m(i)
+    END DO
+
+! Calculation of mean temperature at surface grid points
     ztsol(:) = 0.0
     DO nsrf = 1, nbsrf
@@ -586,16 +632,17 @@
     ENDDO
 
-
-! Repartition du longwave par sous-surface linearisee
+! Linear distrubution on sub-surface of long- and shortwave net radiance
     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))
+          solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
        ENDDO
     ENDDO
 
 
+! Downwelling longwave radiation at mean surface
+    lwdown_m(:) = 0.0
     DO i = 1, klon
-       sollwdown(i) = sollw_m(i) + RSIGMA*ztsol(i)**4
+       lwdown_m(i) = sollw_m(i) + RSIGMA*ztsol(i)**4
     ENDDO
 
@@ -644,22 +691,21 @@
        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)
+          ypct(j)    = pctsrf(i,nsrf)
+          yts(j)     = ts(i,nsrf)
+          ysnow(j)   = snow(i,nsrf)
+          yqsurf(j)  = qsurf(i,nsrf)
+          yalb(j)    = alb(i,nsrf)
+          yalb1(j)   = alb1(i,nsrf)
+          yalb2(j)   = alb2(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)
+          yfder(j)   = fder(i)
+          ysolsw(j)  = solsw(i,nsrf)
+          ysollw(j)  = sollw(i,nsrf)
+          yrugos(j)  = rugos(i,nsrf)
           yrugoro(j) = rugoro(i)
-          yu1(j) = u1lay(i)
-          yv1(j) = v1lay(i)
-          yrads(j) =  ysolsw(j)+ ysollw(j)
+          yu1(j)     = u1lay(i)
+          yv1(j)     = v1lay(i)
           ypaprs(j,klev+1) = paprs(i,klev+1)
           yu10mx(j) = u10m(i,nsrf)
@@ -730,17 +776,9 @@
 !
 !****************************************************************************************
-          
+
+! - Reference pressure is given the values at surface level          
        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
+
+! - Constant CO2 is copied to global grid
        r_co2_ppm(:) = co2_ppm
 
@@ -755,48 +793,53 @@
      
        CASE(is_ter)
+          ! ylwdown : to be removed, calculation is now done at land surface in surf_land
+          ylwdown(:)=0.0
+          DO i=1,knon
+             ylwdown(i)=lwdown_m(ni(i))
+          END DO
           CALL surf_land(itap, dtime, date0, jour, knon, ni,&
                rlon, rlat, &
-               debut, lafin, ydelp(:,1), epot_air, r_co2_ppm, ysollwdown, ysolsw, swdown, &
+               debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
                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)
+               ysnow, yqsol, yagesno, ytsoil, &
+               yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
+               yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, pctsrf_nsrf, &
+               ylwdown)
      
        CASE(is_lic)
           CALL surf_landice(itap, dtime, knon, ni, &
+               ysolsw, ysollw, yts, ypplay(:,1), &
+               ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+               petAcoef, peqAcoef, petBcoef, peqBcoef, &
+               ypsref, yu1, yv1, yrugoro, pctsrf, &
+               ysnow, yqsurf, yqsol, yagesno, &
+               ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
+               ytsurf_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
+          
+       CASE(is_oce)
+          CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb1, &
+               yrugos, ywindsp, rmu0, yfder, &
+               itap, dtime, jour, knon, ni, &
+               debut, &
+               ypplay(:,1), ycoefh(:,1), ycoefm(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+               petAcoef, peqAcoef, petBcoef, peqBcoef, &
+               ypsref, yu1, yv1, yrugoro, pctsrf, &
+               ysnow, yqsurf, yagesno, &
+               yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
+               ytsurf_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
+          
+       CASE(is_sic)
+          CALL surf_seaice( &
+               rlon, rlat, ysolsw, ysollw, yalb1, yfder, &
+               itap, dtime, jour, knon, ni, &
+               debut, lafin, &
                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)
+               ysnow, yqsurf, yqsol, yagesno, ytsoil, &
+               yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
+               ytsurf_new, y_dflux_t, y_dflux_q, pctsrf_nsrf)
           
 
@@ -815,11 +858,8 @@
 !****************************************************************************************
 ! 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)
-
 !****************************************************************************************
 !
@@ -888,6 +928,6 @@
        evap(:,nsrf) = - flux_q(:,1,nsrf)
        
-       albe(:, nsrf) = 0.
-       alblw(:, nsrf) = 0.
+       alb1(:, nsrf) = 0.
+       alb2(:, nsrf) = 0.
        snow(:, nsrf) = 0.
        qsurf(:, nsrf) = 0.
@@ -897,6 +937,6 @@
           i = ni(j)
           d_ts(i,nsrf) = y_d_ts(j)
-          albe(i,nsrf) = yalb(j)  
-          alblw(i,nsrf) = yalblw(j)
+          alb1(i,nsrf) = yalb1_new(j)  
+          alb2(i,nsrf) = yalb2_new(j)
           snow(i,nsrf) = ysnow(j)  
           qsurf(i,nsrf) = yqsurf(j)
Index: LMDZ4/trunk/libf/phylmd/phyetat0.F
===================================================================
--- LMDZ4/trunk/libf/phylmd/phyetat0.F	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/phyetat0.F	(revision 888)
@@ -7,5 +7,5 @@
      .           rlat_p, rlon_p, pctsrf_p, tsol_p,
      .           ocean_in, ok_veget_in,
-     .           albe_p, alblw_p,
+     .           alb1_p, alb2_p,
      .           rain_fall_p, snow_fall_p,solsw_p, sollw_p,
      .           radsol_p,clesphy0,
@@ -53,8 +53,6 @@
       REAL qsol_p(klon)
       REAL snow_p(klon,nbsrf)
-      REAL albe_p(klon,nbsrf)
-cIM BEG alblw
-      REAL alblw_p(klon,nbsrf)
-cIM END alblw
+      REAL alb1_p(klon,nbsrf)   ! albedo in visible SW interval
+      REAL alb2_p(klon,nbsrf)   ! albedo in near IR interval
       REAL evap_p(klon,nbsrf)
       REAL radsol_p(klon)
@@ -93,6 +91,6 @@
       REAL qsol(klon_glo)
       REAL snow(klon_glo,nbsrf)
-      REAL albe(klon_glo,nbsrf)
-      REAL alblw(klon_glo,nbsrf)
+      REAL alb1(klon_glo,nbsrf)
+      REAL alb2(klon_glo,nbsrf)
       REAL evap(klon_glo,nbsrf)
       REAL radsol(klon_glo)
@@ -728,5 +726,5 @@
       ENDIF
 c
-c Lecture de albedo au sol:
+c Lecture de albedo de l'interval visible au sol:
 c
       ierr = NF_INQ_VARID (nid, "ALBE", nvarid)
@@ -746,7 +744,7 @@
            ENDIF
 #ifdef NC_DOUBLE
-           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, albe(1,nsrf))
-#else
-           ierr = NF_GET_VAR_REAL(nid, nvarid, albe(1,nsrf))
+           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, alb1(1,nsrf))
+#else
+           ierr = NF_GET_VAR_REAL(nid, nvarid, alb1(1,nsrf))
 #endif
            IF (ierr.NE.NF_NOERR) THEN
@@ -757,6 +755,6 @@
            xmax = -1.0E+20
            DO i = 1, klon_glo
-              xmin = MIN(albe(i,nsrf),xmin)
-              xmax = MAX(albe(i,nsrf),xmax)
+              xmin = MIN(alb1(i,nsrf),xmin)
+              xmax = MAX(alb1(i,nsrf),xmax)
            ENDDO
            PRINT*,'Albedo du sol ALBE**:', nsrf, xmin, xmax
@@ -766,7 +764,7 @@
          PRINT*, '          J ignore donc les autres ALBE**'
 #ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, albe(1,1))
-#else
-         ierr = NF_GET_VAR_REAL(nid, nvarid, albe(1,1))
+         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, alb1(1,1))
+#else
+         ierr = NF_GET_VAR_REAL(nid, nvarid, alb1(1,1))
 #endif
          IF (ierr.NE.NF_NOERR) THEN
@@ -777,17 +775,17 @@
          xmax = -1.0E+20
          DO i = 1, klon_glo
-            xmin = MIN(albe(i,1),xmin)
-            xmax = MAX(albe(i,1),xmax)
+            xmin = MIN(alb1(i,1),xmin)
+            xmax = MAX(alb1(i,1),xmax)
          ENDDO
          PRINT*,'Neige du sol <ALBE>', xmin, xmax
          DO nsrf = 2, nbsrf
          DO i = 1, klon_glo
-            albe(i,nsrf) = albe(i,1)
-         ENDDO
-         ENDDO
-      ENDIF
-
-c
-c Lecture de albedo au sol LW:
+            alb1(i,nsrf) = alb1(i,1)
+         ENDDO
+         ENDDO
+      ENDIF
+
+c
+c Lecture de albedo au sol dans l'interval proche infra-rouge:
 c
       ierr = NF_INQ_VARID (nid, "ALBLW", nvarid)
@@ -798,5 +796,5 @@
          DO nsrf = 1, nbsrf
            DO i = 1, klon_glo
-             alblw(i,nsrf) = albe(i,nsrf)
+             alb2(i,nsrf) = alb1(i,nsrf)
            ENDDO
          ENDDO
@@ -805,7 +803,7 @@
          PRINT*, '          J ignore donc les autres ALBLW**'
 #ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, alblw(1,1))
-#else
-         ierr = NF_GET_VAR_REAL(nid, nvarid, alblw(1,1))
+         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, alb2(1,1))
+#else
+         ierr = NF_GET_VAR_REAL(nid, nvarid, alb2(1,1))
 #endif
          IF (ierr.NE.NF_NOERR) THEN
@@ -816,11 +814,11 @@
          xmax = -1.0E+20
          DO i = 1, klon_glo
-            xmin = MIN(alblw(i,1),xmin)
-            xmax = MAX(alblw(i,1),xmax)
+            xmin = MIN(alb2(i,1),xmin)
+            xmax = MAX(alb2(i,1),xmax)
          ENDDO
          PRINT*,'Neige du sol <ALBLW>', xmin, xmax
          DO nsrf = 2, nbsrf
          DO i = 1, klon_glo
-            alblw(i,nsrf) = alblw(i,1)
+            alb2(i,nsrf) = alb2(i,1)
          ENDDO
          ENDDO
@@ -1583,6 +1581,6 @@
       call Scatter( qsol,qsol_p)
       call Scatter( snow,snow_p)
-      call Scatter( albe,albe_p)
-      call Scatter( alblw,alblw_p)
+      call Scatter( alb1,alb1_p)
+      call Scatter( alb2,alb2_p)
       call Scatter( evap,evap_p)
       call Scatter( radsol,radsol_p)
Index: LMDZ4/trunk/libf/phylmd/phyredem.F
===================================================================
--- LMDZ4/trunk/libf/phylmd/phyredem.F	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/phyredem.F	(revision 888)
@@ -5,5 +5,5 @@
       SUBROUTINE phyredem (fichnom,dtime,radpas,ocean,
      .           rlat_p,rlon_p, pctsrf_p,tsol_p,
-     .           albedo_p, alblw_p,
+     .           alb1_p, alb2_p,
      .           rain_fall_p, snow_fall_p,solsw_p, sollw_p,
      .           radsol_p,zmea_p,zstd_p,zsig_p,
@@ -46,8 +46,6 @@
       REAL qsol_p(klon)
       REAL snow_p(klon,nbsrf)
-      REAL albedo_p(klon,nbsrf)
-cIM BEG
-      REAL alblw_p(klon,nbsrf)
-cIM END
+      REAL alb1_p(klon,nbsrf)   ! albedo in visible SW interval  
+      REAL alb2_p(klon,nbsrf)   ! albedo in near IR interval
       REAL evap_p(klon,nbsrf)
       REAL rain_fall_p(klon)
@@ -80,6 +78,6 @@
       REAL qsol(klon_glo)
       REAL snow(klon_glo,nbsrf)
-      REAL albedo(klon_glo,nbsrf)
-      REAL alblw(klon_glo,nbsrf)
+      REAL alb1(klon_glo,nbsrf)
+      REAL alb2(klon_glo,nbsrf)
       REAL evap(klon_glo,nbsrf)
       REAL rain_fall(klon_glo)
@@ -147,6 +145,6 @@
       call Gather( qsol_p,qsol)
       call Gather( snow_p,snow)
-      call Gather( albedo_p,albedo)
-      call Gather( alblw_p,alblw)
+      call Gather( alb1_p,alb1)
+      call Gather( alb2_p,alb2)
       call Gather( evap_p,evap)
       call Gather( radsol_p,radsol)
@@ -481,11 +479,10 @@
         ENDIF
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,albedo(1,nsrf))
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,albedo(1,nsrf))
-#endif
-      ENDDO
-
-cIM BEG albedo LW
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,alb1(1,nsrf))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,alb1(1,nsrf))
+#endif
+      ENDDO
+
         DO nsrf = 1, nbsrf
         IF (nsrf.LE.99) THEN
@@ -505,10 +502,10 @@
         ENDIF
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,alblw(1,nsrf))
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,alblw(1,nsrf))
-#endif
-      ENDDO
-cIM END albedo LW
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,alb2(1,nsrf))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,alb2(1,nsrf))
+#endif
+      ENDDO
+c
 c
       DO nsrf = 1, nbsrf
Index: LMDZ4/trunk/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/trunk/libf/phylmd/physiq.F	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/physiq.F	(revision 888)
@@ -690,9 +690,9 @@
       REAL,save ::  solarlong0
 c
-      REAL,allocatable,save :: falbe(:,:)
-c$OMP THREADPRIVATE(falbe)       ! albedo par type de surface
-c
-      REAL,allocatable,save :: falblw(:,:)
-c$OMP THREADPRIVATE(falblw)       ! albedo par type de surface
+      REAL,allocatable,save :: falb1(:,:)
+c$OMP THREADPRIVATE(falb1)       ! albedo par type de surface pour SW visible 
+c
+      REAL,allocatable,save :: falb2(:,:)
+c$OMP THREADPRIVATE(falb2)       ! albedo par type de surface pour SW proche IR
 
 c
@@ -969,10 +969,9 @@
 cIM
 cym      SAVE pctsrf                 ! sous-fraction du sol
-      REAL,allocatable,save :: albsol(:)
-c$OMP THREADPRIVATE(albsol)
-cym      SAVE albsol                 ! albedo du sol total
-      REAL,allocatable,save :: albsollw(:)
-c$OMP THREADPRIVATE(albsollw)      
-cym      SAVE albsollw                 ! albedo du sol total
+
+      REAL,allocatable,save :: albsol1(:) ! albedo du sol total pour SW visible
+c$OMP THREADPRIVATE(albsol1)
+      REAL,allocatable,save :: albsol2(:) ! albedo du sol total pour SW proche IR
+c$OMP THREADPRIVATE(albsol2)      
 
       REAL,allocatable,save :: wo(:,:)
@@ -1527,6 +1526,6 @@
       allocate( ftsol(klon,nbsrf))
       allocate( deltat(klon))
-      allocate( falbe(klon,nbsrf))
-      allocate( falblw(klon,nbsrf))
+      allocate( falb1(klon,nbsrf))
+      allocate( falb2(klon,nbsrf))
       allocate( zmea(klon))
       allocate( zstd(klon))
@@ -1573,6 +1572,6 @@
       allocate( total_rain(klon), nday_rain(klon))
       allocate( pctsrf(klon,nbsrf))
-      allocate( albsol(klon))
-      allocate( albsollw(klon))
+      allocate( albsol1(klon))
+      allocate( albsol2(klon))
       allocate( wo(klon,klev))
       allocate( clwcon(klon,klev),rnebcon(klon,klev))
@@ -1731,4 +1730,5 @@
 c        nhistoW(:,:,:,:) = 0.0
 c        histoW(:,:,:,:) = 0.0
+! fin anne
 
 cIM      
@@ -1773,5 +1773,5 @@
      .       rlat,rlon,pctsrf, ftsol,
      .       ocean, ok_veget,
-     .       falbe, falblw, rain_fall,snow_fall,
+     .       falb1, falb2, rain_fall,snow_fall,
      .       solsw, sollw,
      .       radsol,clesphy0,
@@ -2318,7 +2318,7 @@
      e     t_seri,    q_seri,    u_seri,  v_seri,   
      e     pplay,     paprs,     pctsrf,            
-     +     ftsol,     falbe,     falblw,  u10m,   v10m,
+     +     ftsol,     falb1,     falb2,   u10m,   v10m,
      s     sollwdown, cdragh,    cdragm,  yu1,    yv1,
-     s     albsol,    albsollw,  sens,    evap,  
+     s     albsol1,   albsol2,   sens,    evap,  
      s     zxtsol,    zxfluxlat, zt2m,    qsat2m, 
      s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
@@ -3233,5 +3233,5 @@
      $                          zxsnow_dummy,
      $                          solsw,
-     $                          albsol,
+     $                          albsol1,
      $                          rain_fall,
      $                          snow_fall,
@@ -3294,12 +3294,12 @@
 
       DO i = 1, klon
-         albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
-     .             + falbe(i,is_lic) * pctsrf(i,is_lic)
-     .             + falbe(i,is_ter) * pctsrf(i,is_ter)
-     .             + falbe(i,is_sic) * pctsrf(i,is_sic)
-         albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
-     .               + falblw(i,is_lic) * pctsrf(i,is_lic)
-     .               + falblw(i,is_ter) * pctsrf(i,is_ter)
-     .               + falblw(i,is_sic) * pctsrf(i,is_sic)
+         albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
+     .             + falb1(i,is_lic) * pctsrf(i,is_lic)
+     .             + falb1(i,is_ter) * pctsrf(i,is_ter)
+     .             + falb1(i,is_sic) * pctsrf(i,is_sic)
+         albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
+     .               + falb2(i,is_lic) * pctsrf(i,is_lic)
+     .               + falb2(i,is_ter) * pctsrf(i,is_ter)
+     .               + falb2(i,is_sic) * pctsrf(i,is_sic)
       ENDDO
 
@@ -3335,5 +3335,5 @@
       CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
      e            (dist, rmu0, fract, 
-     e             paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
+     e             paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
      e             wo,
      e             cldfra, cldemi, cldtau,
@@ -3574,5 +3574,5 @@
      I                   pphis,
      I                   pphi,
-     I                   albsol,
+     I                   albsol1,
      I                   qx(1,1,1), 
      I                   rhcl,
@@ -3827,16 +3827,16 @@
 c pour ecriture dans histxxx.nc 
 
-      ! Get some variables from module mod_fonte_neige
+      ! Get some variables from module fonte_neige_mod
       CALL fonte_neige_get_vars(pctsrf, 
      .     zxfqcalving, zxfqfonte, zxffonte)
 
       IF (ocean == 'slab') THEN
-         ! Get some variables from module oceanslab
+         ! Get some variables from module ocean_slab_mod
          CALL ocean_slab_get_vars(tslab, seaice, fluxo, fluxg)
       ELSEIF (ocean == 'couple') THEN
-         ! Get some variables from module oceancpl
+         ! Get some variables from module ocean_cpl_mod
          CALL ocean_cpl_get_vars(fluxo, fluxg)
       ELSE 
-         ! Get some variables from module oceanforced
+         ! Get some variables from module ocean_forced_mod
          CALL ocean_forced_get_vars(fluxo, fluxg)
       ENDIF
@@ -3884,5 +3884,5 @@
          CALL phyredem ("restartphy.nc",dtime,radpas,ocean,
      .      rlat, rlon, pctsrf, ftsol,
-     .      falbe,falblw, rain_fall, 
+     .      falb1, falb2, rain_fall, 
      .      snow_fall,
      .      solsw, sollw,
Index: LMDZ4/trunk/libf/phylmd/radlwsw.F
===================================================================
--- LMDZ4/trunk/libf/phylmd/radlwsw.F	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/radlwsw.F	(revision 888)
@@ -3,5 +3,5 @@
 !
       SUBROUTINE radlwsw(dist, rmu0, fract, 
-     .                  paprs, pplay,tsol,albedo, alblw, t,q,wo,
+     .                  paprs, pplay,tsol,alb1, alb2, t,q,wo,
      .                  cldfra, cldemi, cldtaupd,
      .                  heat,heat0,cool,cool0,radsol,albpla,
@@ -30,5 +30,6 @@
 c pplay----input-R- pression au milieu de couche (Pa)
 c tsol-----input-R- temperature du sol (en K)
-c albedo---input-R- albedo du sol (entre 0 et 1)
+c alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible
+c alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge
 c t--------input-R- temperature (K)
 c q--------input-R- vapeur d'eau (en kg/kg)
@@ -81,5 +82,5 @@
 c
       real paprs(klon,klev+1), pplay(klon,klev)
-      real albedo(klon), alblw(klon), tsol(klon)
+      real alb1(klon), alb2(klon), tsol(klon)
       real t(klon,klev), q(klon,klev), wo(klon,klev)
       real cldfra(klon,klev), cldemi(klon,klev), cldtaupd(klon,klev)
@@ -209,10 +210,10 @@
          zfract(i) = fract(iof+i)
          zrmu0(i) = rmu0(iof+i)
-         PALBD(i,1) = albedo(iof+i)
-!         PALBD(i,2) = albedo(iof+i)
-         PALBD(i,2) = alblw(iof+i)
-         PALBP(i,1) = albedo(iof+i)
-!         PALBP(i,2) = albedo(iof+i)
-         PALBP(i,2) = alblw(iof+i)
+         PALBD(i,1) = alb1(iof+i)
+!         PALBD(i,2) = alb1(iof+i)
+         PALBD(i,2) = alb2(iof+i)
+         PALBP(i,1) = alb1(iof+i)
+!         PALBP(i,2) = alb1(iof+i)
+         PALBP(i,2) = alb2(iof+i)
 cIM cf. JLD pour etre en accord avec ORCHIDEE il faut mettre PEMIS(i) = 0.96
          PEMIS(i) = 1.0 
Index: LMDZ4/trunk/libf/phylmd/surf_land_bucket_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_land_bucket_mod.F90	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/surf_land_bucket_mod.F90	(revision 888)
@@ -20,9 +20,9 @@
   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)
+       spechum, petAcoef, peqAcoef, petBcoef, peqBcoef, pref, &
+       u1_lay, v1_lay, rugoro, swnet, lwnet, &
+       snow, qsol, agesno, tsoil, &
+       qsurf, z0_new, alb1_new, alb2_new, evap, &
+       fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l)
 
 !****************************************************************************************
@@ -47,11 +47,11 @@
     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)       :: pref
     REAL, DIMENSION(klon), INTENT(IN)       :: u1_lay, v1_lay
     REAL, DIMENSION(klon), INTENT(IN)       :: rugoro
+    REAL, DIMENSION(klon), INTENT(IN)       :: swnet, lwnet
 
 ! In/Output variables
 !****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
     REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
@@ -62,7 +62,7 @@
     REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
     REAL, DIMENSION(klon), INTENT(OUT)       :: z0_new
-    REAL, DIMENSION(klon), INTENT(OUT)       :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new, alb2_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
-    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
 
@@ -71,6 +71,7 @@
     REAL, DIMENSION(klon) :: soilcap, soilflux
     REAL, DIMENSION(klon) :: cal, beta, dif_grnd, capsol
-    REAL, DIMENSION(klon) :: alb_neig
+    REAL, DIMENSION(klon) :: alb_neig, alb_lim
     REAL, DIMENSION(klon) :: zfra
+    REAL, DIMENSION(klon) :: radsol       ! total net radiance at surface
     INTEGER               :: i
 !
@@ -79,13 +80,18 @@
 
 !
-!* Read from limit.nc : albedo(alb_new), length of rugosity(z0_new)
+!* Read from limit.nc : albedo(alb_lim), length of rugosity(z0_new)
 !
     CALL interfsur_lim(itime, dtime, jour, &
          knon, knindex, debut,  &
-         alb_new, z0_new)
+         alb_lim, z0_new)
     
 !
 !* Calcultaion of fluxes 
 !
+
+! calculate total absorbed radiance at surface
+       radsol(:) = 0.0
+       radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
+
 ! calculate constants
     CALL calbeta(dtime, is_ter, knon, snow, qsol, beta, capsol, dif_grnd)
@@ -103,5 +109,5 @@
     
     CALL calcul_fluxs(knon, is_ter, dtime, &
-         tsurf, p1lay, cal, beta, tq_cdrag, ps, &
+         tsurf, p1lay, cal, beta, tq_cdrag, pref, &
          precip_rain, precip_snow, snow, qsurf,  &
          radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
@@ -124,10 +130,16 @@
     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))
+       alb_lim(i)  = alb_neig(i) *zfra(i) + alb_lim(i)*(1.0-zfra(i))
     END DO
 
-    alblw(:) = 0.0
-    alblw(1:knon) = alb_new(1:knon)
-
+!
+!* Return albedo : 
+!    alb1_new and alb2_new are here given the same values
+!
+    alb1_new(:) = 0.0
+    alb2_new(:) = 0.0
+    alb1_new(1:knon) = alb_lim(1:knon)
+    alb2_new(1:knon) = alb_lim(1:knon)
+       
 !
 !* Calculate the rugosity
Index: LMDZ4/trunk/libf/phylmd/surf_land_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_land_mod.F90	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/surf_land_mod.F90	(revision 888)
@@ -20,14 +20,16 @@
   SUBROUTINE surf_land(itime, dtime, date0, jour, knon, knindex, &
        rlon, rlat, &
-       debut, lafin, zlev, epot_air, ccanopy, sollwdown, swnet, swdown, &
+       debut, lafin, zlev, ccanopy, swnet, lwnet, albedo, &
        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)
+       pref, u1_lay, v1_lay, rugoro, pctsrf, &
+       snow, qsol, agesno, tsoil, &
+       z0_new, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+       qsurf, tsurf_new, dflux_s, dflux_l, pctsrf_ter, &
+       lwdown_m)
 
     INCLUDE "indicesol.h"
     INCLUDE "dimsoil.h"
+    INCLUDE "YOMCST.h"
 
 ! Input variables  
@@ -39,6 +41,7 @@
     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)       :: zlev, ccanopy
+    REAL, DIMENSION(klon), INTENT(IN)       :: swnet, lwnet
+    REAL, DIMENSION(klon), INTENT(IN)       :: albedo  ! albedo for whole short-wave interval
     REAL, DIMENSION(klon), INTENT(IN)       :: tsurf
     REAL, DIMENSION(klon), INTENT(IN)       :: p1lay
@@ -48,12 +51,14 @@
     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)       :: pref   ! pressure reference
     REAL, DIMENSION(klon), INTENT(IN)       :: u1_lay, v1_lay
     REAL, DIMENSION(klon), INTENT(IN)       :: rugoro
     REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
 
+    REAL, DIMENSION(klon), INTENT(IN)       :: lwdown_m  ! downwelling longwave radiation at mean surface
+                                                         ! corresponds to previous sollwdown
+
 ! In/Output variables
 !****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
     REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
@@ -62,9 +67,11 @@
 ! Output variables
 !****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: z0_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new ! albdeo for shortwave interval 1(visible)
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alb2_new ! albedo for shortwave interval 2(near infrared)
+    REAL, DIMENSION(klon), INTENT(OUT)       :: evap
+    REAL, DIMENSION(klon), INTENT(OUT)       :: fluxsens, fluxlat
     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)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
     REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_ter
@@ -72,5 +79,9 @@
 ! Local variables
 !****************************************************************************************
-    REAL, DIMENSION(klon) :: ps_tmp, p1lay_tmp
+    REAL, DIMENSION(klon) :: p1lay_tmp
+    REAL, DIMENSION(klon) :: pref_tmp
+    REAL, DIMENSION(klon) :: swdown     ! downwelling shortwave radiation at land surface
+    REAL, DIMENSION(klon) :: lwdown     ! downwelling longwave radiation at land surface
+    REAL, DIMENSION(klon) :: epot_air           ! potential air temperature
     REAL, DIMENSION(klon) :: tsol_rad, emis_new ! output from interfsol not used
     INTEGER               :: i
@@ -83,13 +94,30 @@
    IF (ok_veget) THEN
 !****************************************************************************************
-!  Call model sechiba
+!  Call model sechiba in model ORCHIDEE
 !
 !****************************************************************************************
        p1lay_tmp(:)      = 0.0
-       ps_tmp(:)         = 0.0
+       pref_tmp(:)       = 0.0
        p1lay_tmp(1:knon) = p1lay(1:knon)/100.
-       ps_tmp(1:knon)    = ps(1:knon)/100.
+       pref_tmp(1:knon)  = pref(1:knon)/100.
+! 
+!* Calculate incoming flux for SW and LW interval: swdown, lwdown
+!
+       swdown(:) = 0.0
+       lwdown(:) = 0.0
+       DO i = 1, knon
+          swdown(i) = swnet(i)/(1-albedo(i))
+          lwdown(i) = lwnet(i) + RSIGMA*tsurf(i)**4
+       END DO
+!
+!* Calculate potential air temperature
+!
+       epot_air(:) = 0.0
+       DO i = 1, knon
+          epot_air(i) = RCPD*temp_air(i)*(pref(i)/p1lay(i))**RKAPPA
+       END DO
 
 #ifdef CPP_VEGET
+       ! temporary for keeping same results using lwdown_m instead of lwdown
        CALL surf_land_orchidee(itime, dtime, date0, knon, &
             knindex, rlon, rlat, pctsrf, &
@@ -97,13 +125,13 @@
             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, &
+            precip_rain, precip_snow, lwdown_m, swnet, swdown, &
+            pref_tmp, &
             evap, fluxsens, fluxlat, &              
-            tsol_rad, tsurf_new, alb_new, alblw, &
+            tsol_rad, tsurf_new, alb1_new, alb2_new, &
             emis_new, z0_new, qsurf)       
 #endif
 
 !  
-!* Ajout de la contribution du relief
+!* Add contribution of relief to surface roughness
 !  
        DO i=1,knon
@@ -118,9 +146,9 @@
        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)
+            spechum, petAcoef, peqAcoef, petBcoef, peqBcoef, pref, &
+            u1_lay, v1_lay, rugoro, swnet, lwnet, &
+            snow, qsol, agesno, tsoil, &
+            qsurf, z0_new, alb1_new, alb2_new, evap, &
+            fluxsens, fluxlat, tsurf_new, dflux_s, dflux_l)
 
     ENDIF ! ok_veget
Index: LMDZ4/trunk/libf/phylmd/surf_land_orchidee_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_land_orchidee_mod.F90	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/surf_land_orchidee_mod.F90	(revision 888)
@@ -37,5 +37,5 @@
        ps, &
        evap, fluxsens, fluxlat, &              
-       tsol_rad, tsurf_new, alb_new, alblw, &
+       tsol_rad, tsurf_new, alb1_new, alb2_new, &
        emis_new, z0_new, qsurf)
 !    
@@ -84,5 +84,6 @@
 !   tsol_rad     
 !   tsurf_new    temperature au sol
-!   alb_new      albedo
+!   alb1_new     albedo in visible SW interval
+!   alb2_new     albedo in near IR interval
 !   emis_new     emissivite
 !   z0_new       surface roughness
@@ -120,5 +121,6 @@
 !****************************************************************************************
     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)        :: tsol_rad, tsurf_new
+    REAL, DIMENSION(klon), INTENT(OUT)        :: alb1_new, alb2_new
     REAL, DIMENSION(klon), INTENT(OUT)        :: emis_new, z0_new
 
@@ -442,6 +444,6 @@
     ENDIF
 
-    alb_new(1:knon) = albedo_out(1:knon,1) 
-    alblw(1:knon) = albedo_out(1:knon,2)
+    alb1_new(1:knon) = albedo_out(1:knon,1) 
+    alb2_new(1:knon) = albedo_out(1:knon,2)
 
 ! Convention orchidee: positif vers le haut
Index: LMDZ4/trunk/libf/phylmd/surf_landice_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_landice_mod.F90	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/surf_landice_mod.F90	(revision 888)
@@ -17,10 +17,11 @@
 !
   SUBROUTINE surf_landice(itime, dtime, knon, knindex, &
-       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
+       swnet, lwnet, 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)
+       snow, qsurf, qsol, agesno, &
+       tsoil, z0_new, alb1, alb2, evap, fluxsens, fluxlat, &
+       tsurf_new, dflux_s, dflux_l, pctsrf_lic)
 
     INCLUDE "indicesol.h"
@@ -34,4 +35,6 @@
     INTEGER, DIMENSION(klon), INTENT(in)          :: knindex
     REAL, INTENT(in)                              :: dtime
+    REAL, DIMENSION(klon), INTENT(IN)             :: swnet ! net shortwave radiance
+    REAL, DIMENSION(klon), INTENT(IN)             :: lwnet ! net longwave radiance
     REAL, DIMENSION(klon), INTENT(IN)             :: tsurf
     REAL, DIMENSION(klon), INTENT(IN)             :: p1lay
@@ -48,5 +51,4 @@
 ! In/Output variables
 !****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
     REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
@@ -57,7 +59,8 @@
     REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
     REAL, DIMENSION(klon), INTENT(OUT)            :: z0_new
-    REAL, DIMENSION(klon), INTENT(OUT)            :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1  ! new albedo in visible SW interval
+    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2  ! new albedo in near IR interval
     REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
-    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
     REAL, DIMENSION(klon), INTENT(OUT)            :: pctsrf_lic
@@ -68,4 +71,5 @@
     REAL, DIMENSION(klon)    :: cal, beta, dif_grnd
     REAL, DIMENSION(klon)    :: zfra, alb_neig
+    REAL, DIMENSION(klon)    :: radsol
 
 ! End definition
@@ -73,7 +77,13 @@
 !
 ! Initialize output variables
-    alblw(:) =  999999.
-    alb_new(:) = 999999.
+    alb2(:) = 999999.
+    alb1(:) = 999999.
 
+!****************************************************************************************
+! Calculate total absorbed radiance at surface
+!
+!****************************************************************************************
+    radsol(:) = 0.0
+    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
 
 !****************************************************************************************
@@ -122,17 +132,17 @@
     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) + &
+    alb1(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
+!       alb1(1 : knon)  = 0.6 !IM cf FH/GK 
+!       alb1(1 : knon)  = 0.82
+!       alb1(1 : knon)  = 0.77 !211003 Ksta0.77
+!       alb1(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)
+! Attantion: alb1 and alb2 are the same!
+    alb1(1:knon)  = 0.77
+    alb2(1:knon)  = alb1(1:knon)
 
 
Index: LMDZ4/trunk/libf/phylmd/surf_ocean_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_ocean_mod.F90	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/surf_ocean_mod.F90	(revision 888)
@@ -16,15 +16,14 @@
 !****************************************************************************************
 !
-  SUBROUTINE surf_ocean(rlon, rlat, sollw, albedo, &
-       rugos, windsp, rmu0, &
-       fder, &
+  SUBROUTINE surf_ocean(rlon, rlat, swnet, lwnet, alb1, &
+       rugos, windsp, rmu0, fder, &
        itime, dtime, jour, knon, knindex, &
-       debut, swdown, &
+       debut, &
        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)
+       snow, qsurf, agesno, &
+       z0_new, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+       tsurf_new, dflux_s, dflux_l, pctsrf_oce)
 !
 ! This subroutine will make a call to ocean_XXX_noice according to the ocean mode (force, 
@@ -41,11 +40,11 @@
     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)        :: swnet  ! net shortwave radiation at surface  
+    REAL, DIMENSION(klon), INTENT(IN)        :: lwnet  ! net longwave radiation at surface  
+    REAL, DIMENSION(klon), INTENT(IN)        :: alb1   ! albedo in visible SW interval
     REAL, DIMENSION(klon), INTENT(IN)        :: 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
@@ -63,6 +62,6 @@
 ! In/Output variables
 !****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(INOUT)     :: radsol
-    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow, qsurf
+    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
+    REAL, DIMENSION(klon), INTENT(INOUT)     :: qsurf
     REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno
 
@@ -70,7 +69,8 @@
 !****************************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT)       :: z0_new
-    REAL, DIMENSION(klon), INTENT(OUT)       :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new  ! new albedo in visible SW interval
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alb2_new  ! new albedo in near IR interval
     REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
-    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
     REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_oce
@@ -81,38 +81,47 @@
     INTEGER               :: i
     REAL, DIMENSION(klon) :: alb_eau
+    REAL, DIMENSION(klon) :: radsol
 
 ! End definition
 !****************************************************************************************
-! 
+
+
+!****************************************************************************************
+! Calculate total net radiance at surface
+!
+!****************************************************************************************
+    radsol(:) = 0.0
+    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
+
+!****************************************************************************************
 ! Switch according to type of ocean (couple, slab or forced)
 !****************************************************************************************
     SELECT CASE(ocean)
     CASE('couple')
-       CALL ocean_cpl_noice(&
-            sollw, albedo, &
+       CALL ocean_cpl_noice( &
+            swnet, lwnet, alb1, &
             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, &
+            radsol, snow, agesno, &
+            qsurf, evap, fluxsens, fluxlat, &
             tsurf_new, dflux_s, dflux_l, pctsrf_oce)
 
     CASE('slab')
-       CALL ocean_slab_noice(&
+       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)
+            radsol, snow, agesno, &
+            qsurf, evap, fluxsens, fluxlat, &
+            tsurf_new, dflux_s, dflux_l, pctsrf_oce)
        
     CASE('force')
-       CALL ocean_forced_noice(itime, dtime, jour, knon, knindex, &
+       CALL ocean_forced_noice( &
+            itime, dtime, jour, knon, knindex, &
             debut, &
             p1lay, tq_cdrag, precip_rain, precip_snow, &
@@ -120,7 +129,6 @@
             petAcoef, peqAcoef, petBcoef, peqBcoef, &
             ps, u1_lay, v1_lay, &
-            radsol, snow, qsurf, &
-            agesno, &
-            evap, fluxsens, fluxlat, &
+            radsol, snow, agesno, &
+            qsurf, evap, fluxsens, fluxlat, &
             tsurf_new, dflux_s, dflux_l, pctsrf_oce)
     END SELECT
@@ -137,6 +145,7 @@
 
     DO i =1, knon
-       alb_new(i) = alb_eau(knindex(i))
+       alb1_new(i) = alb_eau(knindex(i))
     ENDDO
+    alb2_new(1:knon) = alb1_new(1:knon)
 
 !****************************************************************************************
@@ -154,5 +163,4 @@
     ENDDO
     
-    alblw(1:knon) = alb_new(1:knon)
 !
 !****************************************************************************************
Index: LMDZ4/trunk/libf/phylmd/surf_seaice_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/surf_seaice_mod.F90	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/surf_seaice_mod.F90	(revision 888)
@@ -16,14 +16,13 @@
 !
   SUBROUTINE surf_seaice( & 
-       rlon, rlat, sollw, albedo, &
-       fder, &
+       rlon, rlat, swnet, lwnet, alb1, fder, &
        itime, dtime, jour, knon, knindex, &
-       debut, lafin, swdown, &
+       debut, lafin, &
        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)
+       snow, qsurf, qsol, agesno, tsoil, &
+       z0_new, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+       tsurf_new, dflux_s, dflux_l, pctsrf_sic)
 !
 ! This subroutine will make a call to ocean_XXX_ice according to the ocean mode (force, 
@@ -41,8 +40,8 @@
     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)        :: swnet  ! net shortwave radiation at surface  
+    REAL, DIMENSION(klon), INTENT(IN)        :: lwnet  ! net longwave radiation at surface  
+    REAL, DIMENSION(klon), INTENT(IN)        :: alb1   ! albedo in visible SW interval
     REAL, DIMENSION(klon), INTENT(IN)        :: fder
-    REAL, DIMENSION(klon), INTENT(IN)        :: swdown    
     REAL, DIMENSION(klon), INTENT(IN)        :: tsurf
     REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
@@ -59,5 +58,4 @@
 ! In/Output arguments
 !****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
     REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsurf, qsol
     REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
@@ -67,12 +65,25 @@
 !****************************************************************************************
     REAL, DIMENSION(klon), INTENT(OUT)       :: z0_new
-    REAL, DIMENSION(klon), INTENT(OUT)       :: alblw
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alb1_new  ! new albedo in visible SW interval
+    REAL, DIMENSION(klon), INTENT(OUT)       :: alb2_new  ! new albedo in near IR interval
     REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
-    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new, alb_new
+    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
     REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_sic
 
+! Local arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon)  :: radsol
+!
 ! End definitions
 !****************************************************************************************
+
+
+!****************************************************************************************
+! Calculate total net radiance at surface
+!
+!****************************************************************************************
+    radsol(:) = 0.0
+    radsol(1:knon) = swnet(1:knon) + lwnet(1:knon)
 
 !****************************************************************************************
@@ -83,15 +94,14 @@
     CASE('couple')
        CALL ocean_cpl_ice( &
-          rlon, rlat, sollw, albedo, & 
+          rlon, rlat, swnet, lwnet, alb1, & 
           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)
+          alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+          tsurf_new, dflux_s, dflux_l, pctsrf_sic)
 
     CASE('slab')
@@ -102,8 +112,7 @@
           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)
+          radsol, snow, qsurf, qsol, agesno, tsoil, &
+          alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+          tsurf_new, dflux_s, dflux_l, pctsrf_sic)
     
     CASE('force')
@@ -113,7 +122,7 @@
           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)
+          radsol, snow, qsol, agesno, tsoil, &
+          qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, &
+          tsurf_new, dflux_s, dflux_l, pctsrf_sic)
     END SELECT
 
Index: LMDZ4/trunk/libf/phylmd/write_histday.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/write_histday.h	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/write_histday.h	(revision 888)
@@ -515,5 +515,5 @@
      $                     zx_tmp_fi2d)
 C
-        zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)
+        zx_tmp_fi2d(1 : klon) = falb1( 1 : klon, nsrf)
 cym        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
         CALL histwrite_phy(nid_day,"albe_"//clnsurf(nsrf),itau_w,
Index: LMDZ4/trunk/libf/phylmd/write_histhf.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/write_histhf.h	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/write_histhf.h	(revision 888)
@@ -100,5 +100,5 @@
       CALL histwrite_phy(nid_hf,"SWnetOR",itau_w, zx_tmp_fi2d)
 c
-      zx_tmp_fi2d(1:klon) = solsw(1:klon)/(1.-albsol(1:klon))
+      zx_tmp_fi2d(1:klon) = solsw(1:klon)/(1.-albsol1(1:klon))
 cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d, zx_tmp_2d)
       CALL histwrite_phy(nid_hf,"SWdownOR",itau_w, zx_tmp_fi2d)
Index: LMDZ4/trunk/libf/phylmd/write_histins.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/write_histins.h	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/write_histins.h	(revision 888)
@@ -145,5 +145,5 @@
      $      zx_tmp_fi2d) 
 C
-        zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)
+        zx_tmp_fi2d(1 : klon) = falb1( 1 : klon, nsrf)
 cym        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
         CALL histwrite_phy(nid_ins,"albe_"//clnsurf(nsrf),itau_w,
@@ -151,9 +151,9 @@
 C 
       END DO  
-cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol,zx_tmp_2d)
-      CALL histwrite_phy(nid_ins,"albs",itau_w,albsol)
+cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol1,zx_tmp_2d)
+      CALL histwrite_phy(nid_ins,"albs",itau_w,albsol1)
       
-cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsollw,zx_tmp_2d)
-      CALL histwrite_phy(nid_ins,"albslw",itau_w,albsollw)
+cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol2,zx_tmp_2d)
+      CALL histwrite_phy(nid_ins,"albslw",itau_w,albsol2)
 
 c
Index: LMDZ4/trunk/libf/phylmd/write_histmth.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/write_histmth.h	(revision 887)
+++ LMDZ4/trunk/libf/phylmd/write_histmth.h	(revision 888)
@@ -671,5 +671,5 @@
        DO nsrf=1, nbsrf
 c
-        zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)
+        zx_tmp_fi2d(1 : klon) = falb1( 1 : klon, nsrf)
 cym        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
         CALL histwrite_phy(nid_mth,"albe_"//clnsurf(nsrf),itau_w,
@@ -688,9 +688,9 @@
       ENDDO !nsrf=1, nbsrf
 c
-cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol,zx_tmp_2d)
-      CALL histwrite_phy(nid_mth,"albs",itau_w,albsol)
-c
-cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsollw,zx_tmp_2d)
-      CALL histwrite_phy(nid_mth,"albslw",itau_w,albsollw)
+cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol1,zx_tmp_2d)
+      CALL histwrite_phy(nid_mth,"albs",itau_w,albsol1)
+c
+cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol2,zx_tmp_2d)
+      CALL histwrite_phy(nid_mth,"albslw",itau_w,albsol2)
 c
       ENDIF !lev_histmth.GE.3
@@ -1703,5 +1703,5 @@
        DO nsrf=1, nbsrf
 c
-        zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)
+        zx_tmp_fi2d(1 : klon) = falb1( 1 : klon, nsrf)
 cym        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
         CALL histwrite_phy(nid_mth,"albe_"//clnsurf(nsrf),itau_w,
@@ -1720,9 +1720,9 @@
       ENDDO !nsrf=1, nbsrf
 c
-cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol,zx_tmp_2d)
-      CALL histwrite_phy(nid_mth,"albs",itau_w,albsol)
-c
-cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsollw,zx_tmp_2d)
-      CALL histwrite_phy(nid_mth,"albslw",itau_w,albsollw)
+cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol1,zx_tmp_2d)
+      CALL histwrite_phy(nid_mth,"albs",itau_w,albsol1)
+c
+cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol2,zx_tmp_2d)
+      CALL histwrite_phy(nid_mth,"albslw",itau_w,albsol2)
 c
       ENDIF !lev_histmth.GE.3
