Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/YOMCST.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/YOMCST.h	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/YOMCST.h	(revision 98)
@@ -1,23 +1,23 @@
-C A1.0 Fundamental constants
+! A1.0 Fundamental constants
       REAL RPI,RCLUM,RHPLA,RKBOL,RNAVO
-C A1.1 Astronomical constants
+! A1.1 Astronomical constants
       REAL RDAY,REA,REPSM,RSIYEA,RSIDAY,ROMEGA
-c A1.1.bis Constantes concernant l'orbite de la Terre:
+! A1.1.bis Constantes concernant l'orbite de la Terre:
       REAL R_ecc, R_peri, R_incl
-C A1.2 Geoide
+! A1.2 Geoide
       REAL RA,RG,R1SA
-C A1.3 Radiation
+! A1.3 Radiation
       REAL RSIGMA,RI0
-C A1.4 Thermodynamic gas phase
+! A1.4 Thermodynamic gas phase
       REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV
       REAL RKAPPA,RETV
-C A1.5,6 Thermodynamic liquid,solid phases
+! A1.5,6 Thermodynamic liquid,solid phases
       REAL RCW,RCS
-C A1.7 Thermodynamic transition of phase
+! A1.7 Thermodynamic transition of phase
       REAL RLVTT,RLSTT,RLMLT,RTT,RATM
-C A1.8 Curve of saturation
+! A1.8 Curve of saturation
       REAL RESTT,RALPW,RBETW,RGAMW,RALPS,RBETS,RGAMS
       REAL RALPD,RBETD,RGAMD
-C
+!
       COMMON/YOMCST/RPI   ,RCLUM ,RHPLA ,RKBOL ,RNAVO
      S      ,RDAY  ,REA   ,REPSM ,RSIYEA,RSIDAY,ROMEGA
@@ -31,3 +31,3 @@
      S      ,RESTT ,RALPW ,RBETW ,RGAMW ,RALPS ,RBETS ,RGAMS
      S      ,RALPD ,RBETD ,RGAMD
-C    ------------------------------------------------------------------
+!    ------------------------------------------------------------------
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/clim.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/clim.h	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/clim.h	(revision 98)
@@ -42,5 +42,5 @@
       PARAMETER ( CLIM_Void    = 0  )
       PARAMETER ( CLIM_MaxMod  = 8 )
-      PARAMETER ( CLIM_MaxPort = 16 )
+      PARAMETER ( CLIM_MaxPort = 40 )
       PARAMETER ( CLIM_MaxSegments = 160 )
       PARAMETER ( CLIM_MaxLink = CLIM_MaxMod * CLIM_MaxPort )
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 98)
@@ -1,6 +1,8 @@
       SUBROUTINE clmain(dtime,pctsrf,t,q,u,v,
-     .                  soil_model,ts,soilcap,soilflux,
-     .                  paprs,pplay,radsol,snow,qsol,
-     .                  xlat, rugos,
+     .                  ok_veget,ts,
+     .                  paprs,pplay,radsol,snow,qsol,evap,albe,
+     .                  rain_f, snow_f, solsw, sollw,
+     .                  rlon, rlat, rugos,
+     .                  debut, lafin,
      .                  d_t,d_q,d_u,d_v,d_ts,
      .                  flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
@@ -41,5 +43,5 @@
 c beta-----input-R- coefficient de l'evaporation reelle (0 a 1)
 c dif_grnd-input-R- coeff. de diffusion (chaleur) vers le sol profond
-c xlat-----input-R- latitude en degree
+c rlat-----input-R- latitude en degree
 c rugos----input-R- longeur de rugosite (en m)
 c
@@ -75,12 +77,13 @@
       REAL u(klon,klev), v(klon,klev)
       REAL paprs(klon,klev+1), pplay(klon,klev), radsol(klon)
-      REAL xlat(klon)
+      REAL rlon(klon), rlat(klon)
       REAL d_t(klon, klev), d_q(klon, klev)
       REAL d_u(klon, klev), d_v(klon, klev)
-      REAL flux_t(klon,klev), flux_q(klon,klev)
+      REAL flux_t(klon,klev, nbsrf), flux_q(klon,klev, nbsrf)
       REAL dflux_t(klon), dflux_q(klon)
-      REAL flux_u(klon,klev), flux_v(klon,klev)
+      REAL flux_u(klon,klev, nbsrf), flux_v(klon,klev, nbsrf)
       REAL rugmer(klon)
       REAL cdragh(klon), cdragm(klon)
+      LOGICAL debut, lafin, ok_veget
 cAA      INTEGER itr
 cAA      REAL tr(klon,klev,nbtr)
@@ -93,4 +96,8 @@
       REAL snow(klon,nbsrf)
       REAL qsol(klon,nbsrf)
+      REAL evap(klon,nbsrf)
+      REAL albe(klon,nbsrf)
+      real rain_f(klon), snow_f(klon)
+      REAL sollw(klon), solsw(klon)
       REAL rugos(klon,nbsrf)
 cAA
@@ -103,6 +110,9 @@
 c======================================================================
       REAL yts(klon), yrugos(klon), ypct(klon)
-      REAL ycal(klon), ybeta(klon), ydif(klon)
+      REAL ycal(klon), ybeta(klon), ydif(klon), yalb(klon),yevap(klon)
       REAL yu1(klon), yv1(klon)
+      real ysnow(klon), yqsol(klon)
+      real yrain_f(klon), ysnow_f(klon)
+      real ysollw(klon), ysolsw(klon), ysolswnet(klon)
       REAL yrugm(klon), yrads(klon)
       REAL y_d_ts(klon)
@@ -144,6 +154,4 @@
 c======================================================================
 
-      write(*,*)'CLMAIN.NEW'
-
       DO k = 1, klev   ! epaisseur de couche
       DO i = 1, klon
@@ -174,15 +182,20 @@
          d_ts(i,nsrf) = 0.0
       ENDDO
-      ENDDO
+      END DO
+C§§§ PB
+      flux_t = 0.
+      flux_q = 0.
+      flux_u = 0.
+      flux_v = 0.
       DO k = 1, klev
       DO i = 1, klon
          d_t(i,k) = 0.0
          d_q(i,k) = 0.0
-         flux_t(i,k) = 0.0
-         flux_q(i,k) = 0.0
+c$$$         flux_t(i,k) = 0.0
+c$$$         flux_q(i,k) = 0.0
          d_u(i,k) = 0.0
          d_v(i,k) = 0.0
-         flux_u(i,k) = 0.0
-         flux_v(i,k) = 0.0
+c$$$         flux_u(i,k) = 0.0
+c$$$         flux_v(i,k) = 0.0
          zcoefh(i,k) = 0.0
       ENDDO
@@ -203,21 +216,23 @@
 c
 c prescrire les proprietes du sol:
-      CALL calbeta(dtime,nsrf,snow,qsol, beta,capsol,dif_grnd)
-      IF (.NOT.soil_model) THEN
-         DO i = 1, klon
-            cal(i) = RCPD * capsol(i)
-            totalflu(i) = radsol(i)
-         ENDDO
-      ELSE
-         DO i = 1, klon
-            totalflu(i) = soilflux(i,nsrf) + radsol(i)
-            IF (nsrf.EQ.is_oce) THEN
-               cal(i) = 0.0
-            ELSE
-               cal(i) = RCPD / soilcap(i,nsrf)
-            ENDIF
-         ENDDO
-      ENDIF
-c
+c      CALL calbeta(dtime,nsrf,snow,qsol, beta,capsol,dif_grnd)
+c      IF (.NOT.soil_model) THEN
+c         DO i = 1, klon
+c            cal(i) = RCPD * capsol(i)
+c            totalflu(i) = radsol(i)
+c         ENDDO
+c      ELSE
+c         DO i = 1, klon
+c            totalflu(i) = soilflux(i,nsrf) + radsol(i)
+c            IF (nsrf.EQ.is_oce) THEN
+c               cal(i) = 0.0
+c            ELSE
+c               cal(i) = RCPD / soilcap(i,nsrf)
+c            ENDIF
+c         ENDDO
+c      ENDIF
+c
+      totalflu = radsol
+
 c chercher les indices:
       DO j = 1, klon
@@ -237,11 +252,19 @@
         ypct(j) = pctsrf(i,nsrf)
         yts(j) = ts(i,nsrf)
+        ysnow(j) = snow(i,nsrf)
+        yevap(j) = evap(i,nsrf)
+        yqsol(j) = qsol(i,nsrf)
+        yalb(j) = albe(i,nsrf)
+        yrain_f(j) = rain_f(i)
+        ysnow_f(j) = snow_f(i)
+        ysolsw(j) = solsw(i)
+        ysollw(j) = sollw(i)
         yrugos(j) = rugos(i,nsrf)
         yu1(j) = u1lay(i)
         yv1(j) = v1lay(i)
         yrads(j) = totalflu(i)
-        ycal(j) = cal(i)
-        ybeta(j) = beta(i)
-        ydif(j) = dif_grnd(i)
+c        ycal(j) = cal(i)
+c        ybeta(j) = beta(i)
+c        ydif(j) = dif_grnd(i)
         ypaprs(j,klev+1) = paprs(i,klev+1)
       ENDDO
@@ -259,18 +282,4 @@
       ENDDO
 c
-cAA      IF (itr.GE.1) THEN
-cAA      DO it = 1, itr
-cAA        DO k = 1, klev
-cAA        DO j = 1, knon
-cAA        i = ni(j)
-cAA          ytr(j,k,it) = tr(i,k,it)
-cAA        ENDDO
-cAA        ENDDO
-cAA        DO j = 1, knon
-cAA        i = ni(j)
-cAA          yflxsrf(j,it) = flux_surf(i,it)
-cAA        ENDDO
-cAA      ENDDO
-cAA      ENDIF
 c
 c calculer Cdrag et les coefficients d'echange
@@ -287,45 +296,11 @@
       ENDDO
 c
-c parametrisation non-locale:
-      IF (ok_nonloc) THEN
-         DO i = 1, knon
-            y_cd_h(i) = ycoefh(i,1)
-            y_cd_m(i) = ycoefm(i,1)
-         ENDDO
-         CALL nonlocal(knon, ypaprs, ypplay,
-     .        yts,ybeta,yu,yv,yt,yq,
-     .        y_cd_h, y_cd_m, ycoefm0, ycoefh0, ygamt, ygamq)
-         DO k = 1, klev
-         DO i = 1, knon
-            ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
-            ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
-         ENDDO
-         ENDDO
-      ELSE
-         IF (.NOT.contreg) THEN
-            DO k = 2, klev
-            DO i = 1, knon
-               ygamq(i,k) = 0.0
-               ygamt(i,k) = 0.0
-            ENDDO
-            ENDDO
-         ELSE
-         DO k = 3, klev
-         DO i = 1, knon
-            ygamq(i,k) = 0.0
-            ygamt(i,k) = -1.0E-03
-         ENDDO
-         ENDDO
-         DO i = 1, knon
-            ygamq(i,2) = 0.0
-            ygamt(i,2) = -2.5E-03
-         ENDDO
-         ENDIF
-      ENDIF
 c
 c calculer la diffusion de "q" et de "h"
-      CALL clqh(knon, dtime, nsrf,yu1, yv1,
+      CALL clqh(knon, dtime, nsrf, ni, pctsrf, rlon, rlat,
+     e          yu1, yv1,
      e          ycoefh,yt,yq,yts,ypaprs,ypplay,ydelp,yrads,
-     e          ycal,ybeta,ydif,ygamt,ygamq,
+     e          yevap,yalb, ysnow, yqsol, yrain_f, ysnow_f,
+     e          ysollw, ysolsw,
      s          y_d_t, y_d_q, y_d_ts, 
      s          y_flux_t, y_flux_q, y_dflux_t, y_dflux_q)
@@ -363,21 +338,32 @@
 c
       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)
+C§§§ PB
+          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)
+c$$$ PB        y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
+c$$$ PB        y_flux_q(j,k) = y_flux_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)
+c$$$ PB        y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
+c$$$ PB        y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
+        ENDDO
+      ENDDO
+
+      evap(:,nsrf) = - flux_q(:,1,nsrf)
+c
       DO j = 1, knon
-         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_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
-         y_flux_q(j,k) = y_flux_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)
-         y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
-         y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
-      ENDDO
-      ENDDO
-c
-      DO j = 1, knon
-      i = ni(j)
+         i = ni(j)
          d_ts(i,nsrf) = y_d_ts(j)
+         albe(i,nsrf) = yalb(j)
+         snow(i,nsrf) = ysnow(j)
+         qsol(i,nsrf) = yqsol(j)
          rugmer(i) = yrugm(j)
          cdragh(i) = cdragh(i) + ycoefh(j,1)
@@ -400,10 +386,10 @@
          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)
-         flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
-         flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
+c$$$ PB        flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
+c$$$         flux_q(i,k) = flux_q(i,k) + y_flux_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)
-         flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
-         flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
+c$$$  PB       flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
+c$$$         flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
          zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
       ENDDO
@@ -430,9 +416,13 @@
 99999 CONTINUE
 c
+
       RETURN
       END
-      SUBROUTINE clqh(knon,dtime,nisurf,u1lay,v1lay,coef,
+      SUBROUTINE clqh(knon,dtime,nisurf,knindex,pctsrf, rlon, rlat,
+     e                u1lay,v1lay,coef,
      e                t,q,ts,paprs,pplay,
-     e                delp,radsol,cal,beta,dif_grnd, gamt,gamq,
+     e                delp,radsol,evap,albedo,snow,qsol, 
+     e                precip_rain, precip_snow,
+     e                lwdown, swdown,
      s                d_t, d_q, d_ts, flux_t, flux_q,dflux_s,dflux_l)
 
@@ -446,4 +436,8 @@
 #include "dimensions.h"
 #include "dimphy.h"
+#include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
+#include "indicesol.h"
 c Arguments:
       INTEGER knon
@@ -462,8 +456,16 @@
       REAL q(klon,klev)       ! humidite specifique (kg/kg)
       REAL ts(klon)           ! temperature du sol (K)
+      REAL evap(klon)         ! evaporation au sol
       REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
       REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
       REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)
       REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
+      REAL albedo(klon)       ! albedo de la surface
+      REAL snow(klon)         ! hauteur de neige
+      REAL qsol(klon)         ! humidite de la surface
+      real precip_rain(klon), precip_snow(klon)
+      integer knindex(klon)
+      real pctsrf(klon,nbsrf)
+      real rlon(klon), rlat(klon)
 c
       REAL d_t(klon,klev)     ! incrementation de "t"
@@ -521,8 +523,6 @@
 c======================================================================
       REAL zcor, zdelta, zcvm5
-#include "YOMCST.h"
-#include "YOETHF.h"
-#include "FCTTRE.h"
-#include "indicesol.h"
+      logical contreg
+      parameter (contreg=.true.)
 c======================================================================
 c Rajout pour l'interface
@@ -530,13 +530,10 @@
       integer jour
       integer nisurf
-      integer knindex(klon)
       logical debut, lafin, ok_veget
-      real rlon(klon), rlat(klon)
-      real zlev(klon), zlflu(klon)
+      real zlev1(klon)
       real temp_air(klon), spechum(klon)
       real hum_air(klon), ccanopy(klon)
       real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
       real petBcoef(klon), peqBcoef(klon)
-      real precip_rain(klon), precip_snow(klon)
       real lwdown(klon), swnet(klon), swdown(klon), ps(klon)
       real p1lay(klon)
@@ -545,10 +542,29 @@
 
 ! Parametres de sortie
-      real evap(klon), fluxsens(klon), fluxlat(klon)
+      real fluxsens(klon), fluxlat(klon)
       real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
       real emis_new(klon), z0_new(klon)
-      real dflux_l(klon), dflux_s(klon)
       real pctsrf_new(klon,nbsrf)
 c
+
+      if (.not. contreg) then
+        do k = 2, klev
+          do i = 1, knon
+            gamq(i,k) = 0.0
+            gamt(i,k) = 0.0
+          enddo
+        enddo
+      else
+        do k = 3, klev
+          do i = 1, knon
+            gamq(i,k)= 0.0
+            gamt(i,k)=  -1.0e-03
+          enddo
+        enddo
+        do i = 1, knon
+          gamq(i,2) = 0.0
+          gamt(i,2) = -2.5e-03
+        enddo
+      endif
 
       DO i = 1, knon
@@ -650,19 +666,28 @@
       spechum=q(:,1)
       p1lay = pplay(:,1)
+      zlev1 = delp(:,1)
+      swnet = swdown * (1 - albedo)
+c En attendant mieux
+      hum_air = 0.
+      ccanopy = 0.
+      tq_cdrag = 0.
 
       CALL interfsurf(itime, dtime, jour,
-     . klon, nisurf, knon, knindex, rlon, rlat,
+     . klon, iim, jjm, nisurf, knon, knindex, pctsrf, rlon, rlat,
      . debut, lafin, ok_veget,
-     . zlev, zlflu, u1lay, v1lay, temp_air, spechum, hum_air, ccanopy, 
+     . zlev1,  u1lay, v1lay, temp_air, spechum, hum_air, ccanopy, 
      . tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,
      . precip_rain, precip_snow, lwdown, swnet, swdown,
-     . ts, p1lay, cal, beta, coef1lay, psref, radsol, dif_grnd,
+     . albedo, snow, qsol,
+     . ts, p1lay, coef1lay, psref, radsol,
      . ocean,
      . evap, fluxsens, fluxlat, dflux_l, dflux_s,              
-     . tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new)
+     . tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new, 
+     . zmasq)
 
       flux_t(:,1) = fluxsens
-      flux_q(:,1) = fluxlat
+      flux_q(:,1) = - evap
       d_ts = tsurf_new - ts
+      albedo = alb_new
 
 c==== une fois on a zx_h_ts, on peut faire l'iteration ========
@@ -1228,5 +1253,5 @@
       RETURN
       END
-      SUBROUTINE calbeta(dtime,indice,snow,qsol,
+      SUBROUTINE calbeta(dtime,indice,knon,snow,qsol,
      .                    vbeta,vcal,vdif)
       IMPLICIT none
@@ -1238,6 +1263,6 @@
 c Calculer quelques parametres pour appliquer la couche limite
 c ------------------------------------------------------------
-#include "dimensions.h"
-#include "dimphy.h"
+!#include "dimensions.h"
+!#include "dimphy.h"
 #include "YOMCST.h"
 #include "indicesol.h"
@@ -1256,13 +1281,13 @@
 c
       REAL dtime
-      REAL snow(klon,nbsrf), qsol(klon,nbsrf)
-      INTEGER indice
-C
-      REAL vbeta(klon)
-      REAL vcal(klon)
-      REAL vdif(klon)
+      REAL snow(knon,nbsrf), qsol(knon,nbsrf)
+      INTEGER indice, knon
+C
+      REAL vbeta(knon)
+      REAL vcal(knon)
+      REAL vdif(knon)
 C
       IF (indice.EQ.is_oce) THEN
-      DO i = 1, klon
+      DO i = 1, knon
           vcal(i)   = 0.0
           vbeta(i)  = 1.0
@@ -1272,5 +1297,5 @@
 c
       IF (indice.EQ.is_sic) THEN
-      DO i = 1, klon
+      DO i = 1, knon
           vcal(i) = calice
           IF (snow(i,is_sic) .GT. 0.0) vcal(i) = calsno
@@ -1282,5 +1307,5 @@
 c
       IF (indice.EQ.is_ter) THEN
-      DO i = 1, klon
+      DO i = 1, knon
           vcal(i) = calsol
           IF (snow(i,is_ter) .GT. 0.0) vcal(i) = calsno
@@ -1291,5 +1316,5 @@
 c
       IF (indice.EQ.is_lic) THEN
-      DO i = 1, klon
+      DO i = 1, knon
           vcal(i) = calice
           IF (snow(i,is_lic) .GT. 0.0) vcal(i) = calsno
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/condsurf.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/condsurf.F	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/condsurf.F	(revision 98)
@@ -31,5 +31,5 @@
 
       LOGICAL newlmt
-      PARAMETER (newlmt=.FALSE.)
+      PARAMETER (newlmt=.TRUE.)
 
       INTEGER     nannemax
@@ -114,21 +114,23 @@
       IF (newlmt) THEN
 c
-c Fraction "ocean":
-      ierr = NF_INQ_VARID (nid, "FOCE", nvarid)
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "condsurf: Le champ <FOCE> est absent"
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_oce))
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_oce))
-#endif
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "condsurf: Lecture echouee pour <FOCE>"
-         CALL abort
-      ENDIF
+c$$$c Fraction "ocean":
+c$$$      ierr = NF_INQ_VARID (nid, "FOCE", nvarid)
+c$$$      IF (ierr .NE. NF_NOERR) THEN
+c$$$         PRINT*, "condsurf: Le champ <FOCE> est absent"
+c$$$         CALL abort
+c$$$      ENDIF
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_oce))
+c$$$#else
+c$$$      ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_oce))
+c$$$#endif
+c$$$      IF (ierr .NE. NF_NOERR) THEN
+c$$$         PRINT*, "condsurf: Lecture echouee pour <FOCE>"
+c$$$         CALL abort
+c$$$      ENDIF
 c
 c Fraction "glace de mer":
+c
+c
       ierr = NF_INQ_VARID (nid, "FSIC", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
@@ -145,36 +147,54 @@
          CALL abort
       ENDIF
-c
-c Fraction "terre":
-      ierr = NF_INQ_VARID (nid, "FTER", nvarid)
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "condsurf: Le champ <FTER> est absent"
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_ter))
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_ter))
-#endif
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "condsurf: Lecture echouee pour <FTER>"
-         CALL abort
-      ENDIF
-c
-c Fraction "glacier terre":
-      ierr = NF_INQ_VARID (nid, "FLIC", nvarid)
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "condsurf: Le champ <FLIC> est absent"
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_lic))
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_lic))
-#endif
-      IF (ierr .NE. 0) THEN
-         PRINT*, "condsurf: Lecture echouee pour <FLIC>"
-         CALL abort
-      ENDIF
+C 
+C positionnement % ocean libre et verification qu'il y a compatibilite des soussurfaces
+      pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon)) 
+     $    - pctsrf(1 : klon, is_sic)
+      DO i = 1, klon
+        IF ( pctsrf(i, is_sic) .GT. (1. - zmasq(i)) ) THEN
+            WRITE(*,*) 'condsurf : sea-ice et masque pb en ', i,
+     $     pctsrf(i, is_sic), (1. - zmasq(i))     
+            pctsrf(i, is_sic) = (1. - zmasq(i))
+            pctsrf(i, is_oce) = 0.
+        ENDIF          
+        IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) + 
+     $       pctsrf(i, is_oce) + pctsrf(i, is_sic)  - 1.) .GT. EPSFRA) 
+     $       THEN 
+             WRITE(*,*) 'physiq : pb sous surface au point ', i, 
+     $           pctsrf(i, 1 : nbsrf)
+         ENDIF 
+      END DO 
+c
+c$$$c Fraction "terre":
+c$$$      ierr = NF_INQ_VARID (nid, "FTER", nvarid)
+c$$$      IF (ierr .NE. NF_NOERR) THEN
+c$$$         PRINT*, "condsurf: Le champ <FTER> est absent"
+c$$$         CALL abort
+c$$$      ENDIF
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_ter))
+c$$$#else
+c$$$      ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_ter))
+c$$$#endif
+c$$$      IF (ierr .NE. NF_NOERR) THEN
+c$$$         PRINT*, "condsurf: Lecture echouee pour <FTER>"
+c$$$         CALL abort
+c$$$      ENDIF
+c$$$c
+c$$$c Fraction "glacier terre":
+c$$$      ierr = NF_INQ_VARID (nid, "FLIC", nvarid)
+c$$$      IF (ierr .NE. NF_NOERR) THEN
+c$$$         PRINT*, "condsurf: Le champ <FLIC> est absent"
+c$$$         CALL abort
+c$$$      ENDIF
+c$$$#ifdef NC_DOUBLE
+c$$$      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,debut,epais,pctsrf(1,is_lic))
+c$$$#else
+c$$$      ierr = NF_GET_VARA_REAL(nid,nvarid,debut,epais,pctsrf(1,is_lic))
+c$$$#endif
+c$$$      IF (ierr .NE. 0) THEN
+c$$$         PRINT*, "condsurf: Lecture echouee pour <FLIC>"
+c$$$         CALL abort
+c$$$      ENDIF
 c
       ELSE ! test sur newlmt
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/dimphy.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/dimphy.h	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/dimphy.h	(revision 98)
@@ -7,2 +7,4 @@
       PARAMETER (nbtr=nqmx-2+1/(nqmx-1))
 c-----------------------------------------------------------------------
+      REAL zmasq(KLON)
+      COMMON/terreoce/zmasq
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/inc_cpl.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/inc_cpl.h	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/inc_cpl.h	(revision 98)
@@ -1,16 +1,19 @@
 C
-C -- inc_cpl.h  
+C -- inc_cpl.h   1998-04
 C    **********
 C@
-C@  Contents : variables describing field restart file names 
+C@  Contents : variables describing pipe and field names
 C@  --------
 C@
-C@ -- cl_write/cl_read  : for fields to write/READ
-C@ -- cl_f_write/cl_f_read  : for fields to write/read
+C@ -- cl_write  : for fields to write
+C@
+C@ -- cl_read  : for fields to read
 C@
 C     -------------------------------------------------------------------
 C
+      INTEGER jpread, jpwrit
+      PARAMETER (jpread=0, jpwrit=1)
       CHARACTER*8 cl_writ(jpmaxfld), cl_read(jpmaxfld)
-      CHARACTER*8 cl_f_writ(jpmaxfld), cl_f_read(jpmaxfld)
+      CHARACTER*6 cl_f_writ(jpmaxfld), cl_f_read(jpmaxfld)
       COMMON / comcpl / cl_writ, cl_read, cl_f_writ, cl_f_read
 C     -------------------------------------------------------------------
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/inc_sipc.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/inc_sipc.h	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/inc_sipc.h	(revision 98)
@@ -18,6 +18,6 @@
       INTEGER  mpoolinitr
       INTEGER  mpoolinitw
-      INTEGER  mpoolwrit(jpflda2o)
-      INTEGER  mpoolread(jpfldo2a)
+      INTEGER  mpoolwrit(jpmaxfld)
+      INTEGER  mpoolread(jpmaxfld)
       COMMON / compool / mpoolinitr, mpoolinitw, mpoolwrit, mpoolread
 C     -------------------------------------------------------------------
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/indicesol.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/indicesol.h	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/indicesol.h	(revision 98)
@@ -13,2 +13,5 @@
       REAL epsfra
       PARAMETER (epsfra=1.0E-05)
+!
+      CHARACTER *3 clnsurf(nbsrf)
+      DATA clnsurf/'ter', 'lic', 'oce', 'sic'/
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 98)
@@ -45,13 +45,14 @@
 !
   SUBROUTINE interfsurf_hq(itime, dtime, jour, &
-      & klon, nisurf, knon, knindex, rlon, rlat, &
+      & klon, iim, jjm, nisurf, knon, knindex, pctsrf, rlon, rlat, &
       & debut, lafin, ok_veget, &
-      & zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, & 
+      & zlev,  u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, & 
       & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
       & precip_rain, precip_snow, lwdown, swnet, swdown, &
-      & tsurf, p1lay, cal, beta, coef1lay, ps, radsol, dif_grnd, &
+      & albedo, snow, qsol, &
+      & tsurf, p1lay, coef1lay, ps, radsol, &
       & ocean, &
       & evap, fluxsens, fluxlat, dflux_l, dflux_s, &              
-      & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new)
+      & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new, zmasq)
 
 
@@ -67,4 +68,5 @@
 !   itime        numero du pas de temps
 !   klon         nombre total de points de grille
+!   iim, jjm     nbres de pts de grille
 !   dtime        pas de temps de la physique (en s)
 !   jour         jour dans l'annee en cours
@@ -72,4 +74,5 @@
 !   knon         nombre de points de la surface a traiter
 !   knindex      index des points de la surface a traiter
+!   pctsrf       tableau des pourcentages de surface de chaque maille
 !   rlon         longitudes
 !   rlat         latitudes
@@ -79,5 +82,312 @@
 !                     (si false calcul simplifie des fluxs sur les continents)
 !   zlev         hauteur de la premiere couche
-!   zlflu        epaisseur de la premier couche
+!   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
+!   hum_air      humidite 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 entrant a la surface
+!   swnet        flux solaire net
+!   swdown       flux solaire entrant a la surface
+!   albedo       albedo de la surface
+!   tsurf        temperature de surface
+!   p1lay        pression 1er niveau (milieu de couche)
+!   coef1lay     coefficient d'echange
+!   ps           pression au sol
+!   radsol       rayonnement net aus sol (LW + SW)
+!   ocean        type d'ocean utilise (force, slab, couple)
+!
+! 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
+!   pctsrf_new   nouvelle repartition des surfaces
+!   zmasq        masque terre/ocean
+
+  include 'indicesol.h'
+
+! Parametres d'entree
+  integer, intent(IN) :: itime
+  integer, intent(IN) :: iim, jjm
+  integer, intent(IN) :: klon
+  real, intent(IN) :: dtime
+  integer, intent(IN) :: jour
+  integer, intent(IN) :: nisurf
+  integer, intent(IN) :: knon
+  integer, dimension(knon), intent(in) :: knindex
+  real, dimension(klon,nbsrf), intent(IN) :: pctsrf
+  logical, intent(IN) :: debut, lafin, ok_veget
+  real, dimension(klon), intent(IN) :: rlon, rlat
+  real, dimension(knon), intent(IN) :: zlev
+  real, dimension(knon), intent(IN) :: u1_lay, v1_lay
+  real, dimension(knon), intent(IN) :: temp_air, spechum
+  real, dimension(knon), intent(IN) :: hum_air, ccanopy
+  real, dimension(knon), intent(IN) :: tq_cdrag, petAcoef, peqAcoef
+  real, dimension(knon), intent(IN) :: petBcoef, peqBcoef
+  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
+  real, dimension(knon), intent(IN) :: lwdown, swnet, swdown, ps, albedo
+  real, dimension(knon), intent(IN) :: tsurf, p1lay, coef1lay
+  real, dimension(knon), intent(IN) :: radsol
+  real, dimension(klon), intent(IN) :: zmasq
+  character (len = 6)  :: ocean
+  real, dimension(knon), intent(INOUT) :: evap, snow, qsol
+
+! Parametres de sortie
+  real, dimension(knon), intent(OUT):: fluxsens, fluxlat
+  real, dimension(knon), intent(OUT):: tsol_rad, tsurf_new, alb_new
+  real, dimension(knon), intent(OUT):: emis_new, z0_new
+  real, dimension(knon), intent(OUT):: dflux_l, dflux_s
+  real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
+
+! Local
+  character (len = 20) :: modname = 'interfsurf_hq'
+  character (len = 80) :: abort_message 
+  logical, save        :: first_call = .true.
+  integer              :: error
+  logical              :: check = .true.
+  real, dimension(knon):: cal, beta, dif_grnd, capsol
+  real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=1./86400.*5.
+  real, parameter      :: calsno=1./(2.3867e+06*.15)
+
+#include "YOMCST.inc"
+
+  if (check) write(*,*) 'Entree ', modname
+!
+! On doit commencer par appeler les schemas de surfaces continentales
+! car l'ocean a besoin du ruissellement qui est y calcule
+!
+  if (first_call) then
+    if (nisurf /= is_ter .and. klon > 1) then 
+      write(*,*)' *** Warning ***'
+      write(*,*)' nisurf = ',nisurf,' /= is_ter = ',is_ter
+      write(*,*)'or on doit commencer par les surfaces continentales'
+      abort_message='voir ci-dessus'
+      call abort_gcm(modname,abort_message,1)
+    endif
+    if (ocean /= 'slab  ' .and. ocean /= 'force ' .and. ocean /= 'couple') then
+      write(*,*)' *** Warning ***'
+      write(*,*)'Option couplage pour l''ocean = ', ocean
+      abort_message='option pour l''ocean non valable'
+      call abort_gcm(modname,abort_message,1)
+    endif
+  endif
+  first_call = .false.
+!
+! Calcul age de la neige
+!
+  
+! Aiguillage vers les differents schemas de surface
+
+  if (nisurf == is_ter) then
+!
+! Surface "terre" appel a l'interface avec les sols continentaux
+!
+! allocation du run-off
+    if (.not. allocated(run_off)) then
+      allocate(run_off(knon), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation run_off'
+        call abort_gcm(modname,abort_message,1)
+      endif
+    else if (size(run_off) /= knon) then
+      write(*,*)'Bizarre, le nombre de points continentaux'
+      write(*,*)'a change entre deux appels. Je continue ...'
+      deallocate(run_off, stat = error)
+      allocate(run_off(knon), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation run_off'
+        call abort_gcm(modname,abort_message,1)
+      endif
+    endif
+!
+    if (.not. ok_veget) then
+!
+! calcul snow et qsol, hydrol adapté
+!
+      call calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
+      cal = RCPD * capsol
+      call calcul_fluxs( knon, nisurf, dtime, &
+     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
+     &   precip_rain, precip_snow, snow, qsol,  &
+     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
+     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
+!
+! calcul albedo: lecture albedo fichier CL puis ajout albedo neige 
+! 
+    else
+!
+!  appel a sechiba
+!
+      call interfsol(itime, klon, dtime, nisurf, knon, &
+     &  knindex, rlon, rlat, &
+     &  debut, lafin, ok_veget, &
+     &  zlev,  u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, & 
+     &  tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
+     &  precip_rain, precip_snow, lwdown, swnet, swdown, &
+     &  tsurf, p1lay, coef1lay, ps, radsol, &
+     &  evap, fluxsens, fluxlat, &              
+     &  tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
+    endif    
+!
+  else if (nisurf == is_oce) then
+
+    if (check) write(*,*)'ocean, nisurf = ',nisurf 
+
+
+!
+! Surface "ocean" appel a l'interface avec l'ocean
+!
+!    if (ocean == 'couple') then
+!      call interfoce(nisurf, ocean)
+!    else if (ocean == 'slab  ') then
+!      call interfoce(nisurf)
+!    else                              ! lecture conditions limites
+!      call interfoce(itime, dtime, jour, & 
+!     &  klon, nisurf, knon, knindex, &
+!     &  debut, &
+!     &  tsurf_new, alb_new, z0_new, pctsrf_new)
+!    endif
+!
+    cal = 0.
+    beta = 1.
+    dif_grnd = 0.
+
+    call calcul_fluxs( knon, nisurf, dtime, &
+     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
+     &   precip_rain, precip_snow, snow, qsol,  &
+     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
+     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
+
+!
+! calcul albedo
+!
+
+!
+  else if (nisurf == is_sic) then
+
+    if (check) write(*,*)'sea ice, nisurf = ',nisurf 
+
+!
+! Surface "glace de mer" appel a l'interface avec l'ocean
+!
+!    call interfoce(nisurf, ocean)
+!
+
+    cal = calice
+    where (snow > 0.0) cal = calsno
+    beta = 1.0
+    dif_grnd = 1.0 / tau_gl
+
+    call calcul_fluxs( knon, nisurf, dtime, &
+     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
+     &   precip_rain, precip_snow, snow, qsol,  &
+     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
+     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
+
+  else if (nisurf == is_lic) then
+
+    if (check) write(*,*)'glacier, nisurf = ',nisurf 
+
+!
+! Surface "glacier continentaux" appel a l'interface avec le sol
+!
+!    call interfsol(nisurf)
+
+    cal = calice
+    where (snow > 0.0) cal = calsno
+    beta = 1.0
+    dif_grnd = 0.0
+
+    call calcul_fluxs( knon, nisurf, dtime, &
+     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
+     &   precip_rain, precip_snow, snow, qsol,  &
+     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
+     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
+  else
+    write(*,*)'Index surface = ',nisurf
+    abort_message = 'Index surface non valable'
+    call abort_gcm(modname,abort_message,1)
+  endif
+
+  END SUBROUTINE interfsurf_hq
+
+!
+!#########################################################################
+!
+  SUBROUTINE interfsurf_vent(nisurf, knon   &         
+  &                     )
+!
+! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general 
+! (sols continentaux, oceans, glaces) pour les tensions de vents.
+! En pratique l'interface se fait entre la couche limite du modele 
+! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
+!
+! 
+! L.Fairhead 02/2000
+!
+! input:
+!   nisurf       index de la surface a traiter (1 = sol continental)
+!   knon         nombre de points de la surface a traiter
+
+! Parametres d'entree
+  integer, intent(IN) :: nisurf
+  integer, intent(IN) :: knon
+
+
+  return
+  END SUBROUTINE interfsurf_vent
+!
+!#########################################################################
+!
+  SUBROUTINE interfsol(itime, klon, dtime, nisurf, knon, &
+     & knindex, rlon, rlat, &
+     & debut, lafin, ok_veget, &
+     & zlev,  u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, & 
+     & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
+     & precip_rain, precip_snow, lwdown, swnet, swdown, &
+     & tsurf, p1lay, coef1lay, ps, radsol, &
+     & evap, fluxsens, fluxlat, &              
+     & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
+
+! 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
+!   klon         nombre total de points de grille
+!   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
+!   debut        logical: 1er appel a la physique (lire les restart)
+!   lafin        logical: dernier appel a la physique (ecrire les restart)
+!   ok_veget     logical: appel ou non au schema de surface continental
+!                     (si false calcul simplifie des fluxs sur les continents)
+!   zlev         hauteur de la premiere couche       
 !   u1_lay       vitesse u 1ere couche
 !   v1_lay       vitesse v 1ere couche
@@ -98,11 +408,11 @@
 !   tsurf        temperature de surface
 !   p1lay        pression 1er niveau (milieu de couche)
-!   cal          capacite calorifique du sol
-!   beta         evap reelle
 !   coef1lay     coefficient d'echange
 !   ps           pression au sol
 !   radsol       rayonnement net aus sol (LW + SW)
-!   dif_grnd     coeff. diffusion vers le sol profond
-!   ocean        type d'ocean utilise (force, slab, couple)
+!   
+!
+! input/output
+!   run_off      ruissellement total
 !
 ! output:
@@ -116,255 +426,4 @@
 !   z0_new       surface roughness
 
-  include 'indicesol.h'
-
-! Parametres d'entree
-  integer, intent(IN) :: itime
-  integer, intent(IN) :: klon
-  real, intent(IN) :: dtime
-  integer, intent(IN) :: jour
-  integer, intent(IN) :: nisurf
-  integer, intent(IN) :: knon
-  integer, dimension(knon), intent(in) :: knindex
-  logical, intent(IN) :: debut, lafin, ok_veget
-  real, dimension(klon), intent(IN) :: rlon, rlat
-  real, dimension(knon), intent(IN) :: zlev, zlflu
-  real, dimension(knon), intent(IN) :: u1_lay, v1_lay
-  real, dimension(knon), intent(IN) :: temp_air, spechum
-  real, dimension(knon), intent(IN) :: hum_air, ccanopy
-  real, dimension(knon), intent(IN) :: tq_cdrag, petAcoef, peqAcoef
-  real, dimension(knon), intent(IN) :: petBcoef, peqBcoef
-  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
-  real, dimension(knon), intent(IN) :: lwdown, swnet, swdown, ps
-  real, dimension(knon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
-  real, dimension(knon), intent(IN) :: radsol, dif_grnd
-  character (len = 6)  :: ocean
-
-! Parametres de sortie
-  real, dimension(knon), intent(OUT):: evap, fluxsens, fluxlat
-  real, dimension(knon), intent(OUT):: tsol_rad, tsurf_new, alb_new
-  real, dimension(knon), intent(OUT):: emis_new, z0_new
-  real, dimension(knon), intent(OUT):: dflux_l, dflux_s
-  real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
-
-! Local
-  character (len = 20) :: modname = 'interfsurf_hq'
-  character (len = 80) :: abort_message 
-  logical, save        :: first_call = .true.
-  integer              :: error
-  logical              :: check = .true.
-
-  if (check) write(*,*) 'Entree ', modname
-!
-! On doit commencer par appeler les schemas de surfaces continentales
-! car l'ocean a besoin du ruissellement qui est y calcule
-!
-  if (first_call) then
-    if (nisurf /= is_ter .and. klon > 1) then 
-      write(*,*)' *** Warning ***'
-      write(*,*)' nisurf = ',nisurf,' /= is_ter = ',is_ter
-      write(*,*)'or on doit commencer par les surfaces continentales'
-      abort_message='voir ci-dessus'
-      call abort_gcm(modname,abort_message,1)
-    endif
-    if (ocean /= 'slab  ' .and. ocean /= 'force ' .and. ocean /= 'couple') then
-      write(*,*)' *** Warning ***'
-      write(*,*)'Option couplage pour l''ocean = ', ocean
-      abort_message='option pour l''ocean non valable'
-      call abort_gcm(modname,abort_message,1)
-    endif
-  endif
-  first_call = .false.
-
-! Aiguillage vers les differents schemas de surface
-
-  if (nisurf == is_ter) then
-!
-! Surface "terre" appel a l'interface avec les sols continentaux
-!
-! allocation du run-off
-    if (.not. allocated(run_off)) then
-      allocate(run_off(knon), stat = error)
-      if (error /= 0) then
-        abort_message='Pb allocation run_off'
-        call abort_gcm(modname,abort_message,1)
-      endif
-    else if (size(run_off) /= knon) then
-      write(*,*)'Bizarre, le nombre de points continentaux'
-      write(*,*)'a change entre deux appels. Je continue ...'
-      deallocate(run_off, stat = error)
-      allocate(run_off(knon), stat = error)
-      if (error /= 0) then
-        abort_message='Pb allocation run_off'
-        call abort_gcm(modname,abort_message,1)
-      endif
-    endif
-!
-    call interfsol(itime, klon, dtime, nisurf, knon, &
-     &  knindex, rlon, rlat, &
-     &  debut, lafin, ok_veget, &
-     &  zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, & 
-     &  tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
-     &  precip_rain, precip_snow, lwdown, swnet, swdown, &
-     &  tsurf, p1lay, cal, beta, coef1lay, ps, radsol, dif_grnd, &
-     &  evap, fluxsens, fluxlat, &              
-     &  tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
-      
-!
-  else if (nisurf == is_oce) then
-
-    if (check) write(*,*)'ocean, nisurf = ',nisurf 
-!
-! Surface "ocean" appel a l'interface avec l'ocean
-!
-!    if (ocean == 'couple') then
-!      call interfoce(nisurf, ocean)
-!    else if (ocean == 'slab  ') then
-!      call interfoce(nisurf)
-!    else                              ! lecture conditions limites
-!      call interfoce(itime, dtime, jour, & 
-!     &  klon, nisurf, knon, knindex, &
-!     &  debut, &
-!     &  tsurf_new, alb_new, z0_new, pctsrf_new)
-!    endif
-!
-    call calcul_fluxs( knon, dtime, &
-     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
-     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
-     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
-     &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
-!
-  else if (nisurf == is_sic) then
-
-    if (check) write(*,*)'sea ice, nisurf = ',nisurf 
-
-!
-! Surface "glace de mer" appel a l'interface avec l'ocean
-!
-!    call interfoce(nisurf, ocean)
-!
-    call calcul_fluxs( knon, dtime, &
-     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
-     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
-     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
-     &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
-
-  else if (nisurf == is_lic) then
-
-    if (check) write(*,*)'glacier, nisurf = ',nisurf 
-
-!
-! Surface "glacier continentaux" appel a l'interface avec le sol
-!
-!    call interfsol(nisurf)
-    call calcul_fluxs( knon, dtime, &
-     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
-     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
-     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
-     &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
-
-  else
-    write(*,*)'Index surface = ',nisurf
-    abort_message = 'Index surface non valable'
-    call abort_gcm(modname,abort_message,1)
-  endif
-
-  END SUBROUTINE interfsurf_hq
-
-!
-!#########################################################################
-!
-  SUBROUTINE interfsurf_vent(nisurf, knon   &         
-  &                     )
-!
-! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general 
-! (sols continentaux, oceans, glaces) pour les tensions de vents.
-! En pratique l'interface se fait entre la couche limite du modele 
-! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
-!
-! 
-! L.Fairhead 02/2000
-!
-! input:
-!   nisurf       index de la surface a traiter (1 = sol continental)
-!   knon         nombre de points de la surface a traiter
-
-! Parametres d'entree
-  integer, intent(IN) :: nisurf
-  integer, intent(IN) :: knon
-
-
-  return
-  END SUBROUTINE interfsurf_vent
-!
-!#########################################################################
-!
-  SUBROUTINE interfsol(itime, klon, dtime, nisurf, knon, &
-     & knindex, rlon, rlat, &
-     & debut, lafin, ok_veget, &
-     & zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, & 
-     & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
-     & precip_rain, precip_snow, lwdown, swnet, swdown, &
-     & tsurf, p1lay, cal, beta, coef1lay, ps, radsol, dif_grnd, &
-     & evap, fluxsens, fluxlat, &              
-     & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
-
-! 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
-!   klon         nombre total de points de grille
-!   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
-!   debut        logical: 1er appel a la physique (lire les restart)
-!   lafin        logical: dernier appel a la physique (ecrire les restart)
-!   ok_veget     logical: appel ou non au schema de surface continental
-!                     (si false calcul simplifie des fluxs sur les continents)
-!   zlev         hauteur de la premiere couche
-!   zlflu        
-!   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
-!   hum_air      humidite 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 entrant a la surface
-!   swnet        flux solaire net
-!   swdown       flux solaire entrant a la surface
-!   tsurf        temperature de surface
-!   p1lay        pression 1er niveau (milieu de couche)
-!   cal          capacite calorifique du sol
-!   beta         evap reelle
-!   coef1lay     coefficient d'echange
-!   ps           pression au sol
-!   radsol       rayonnement net aus sol (LW + SW)
-!   dif_grnd     coeff. diffusion vers le sol profond      
-!   
-!
-! input/output
-!   run_off      ruissellement total
-!
-! 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
 
 ! Parametres d'entree
@@ -377,5 +436,5 @@
   logical, intent(IN) :: debut, lafin, ok_veget
   real, dimension(klon), intent(IN) :: rlon, rlat
-  real, dimension(knon), intent(IN) :: zlev, zlflu
+  real, dimension(knon), intent(IN) :: zlev
   real, dimension(knon), intent(IN) :: u1_lay, v1_lay
   real, dimension(knon), intent(IN) :: temp_air, spechum
@@ -385,6 +444,6 @@
   real, dimension(knon), intent(IN) :: precip_rain, precip_snow
   real, dimension(knon), intent(IN) :: lwdown, swnet, swdown, ps
-  real, dimension(knon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
-  real, dimension(knon), intent(IN) :: radsol, dif_grnd
+  real, dimension(knon), intent(IN) :: tsurf, p1lay, coef1lay
+  real, dimension(knon), intent(IN) :: radsol
 ! Parametres de sortie
   real, dimension(knon), intent(OUT):: evap, fluxsens, fluxlat
@@ -400,4 +459,5 @@
   character (len = 80) :: abort_message
   logical              :: check = .true.
+  real, dimension(knon) :: cal, beta, dif_grnd, capsol
 ! type de couplage dans sechiba
 !  character (len=10)   :: coupling = 'implicit' 
@@ -414,15 +474,10 @@
   integer, save          :: rest_id_stom, hist_id_stom
 
+  real, dimension(knon):: snow, qsol
+
   if (check) write(*,*)'Entree ', modname
   if (check) write(*,*)'ok_veget = ',ok_veget
 
 ! initialisation
-  if (.not. ok_veget) then
-    call calcul_fluxs( knon, dtime, &
-     &   tsurf, p1lay, cal, beta, coef1lay, ps, &
-     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
-     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
-     &   tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
-  else
 !    if (debut) then
 !      !
@@ -476,5 +531,5 @@
 !     & debut, lafin, coupling, control_in, &
 !     & lalo, neighbours, resolution,&
-!     & zlev, zlflu, u1_lay, v1_lay, spechum, temp_air,hum_air , ccanopy, &
+!     & zlev,  u1_lay, v1_lay, spechum, temp_air,hum_air , ccanopy, &
 !     & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
 !     & precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
@@ -486,5 +541,4 @@
 ! Sauvegarde dans fichiers histoire
 !
-  endif  ! fin ok_veget
 
   END SUBROUTINE interfsol
@@ -492,5 +546,12 @@
 !#########################################################################
 !
-  SUBROUTINE interfoce_cpl(nisurf, ocean)
+  SUBROUTINE interfoce_cpl(itime, dtime, &
+      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
+      & ocean, nexca, debut, lafin, &
+      & swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
+      & fder, albsol, taux, tauy, &
+      & tsurf_new, alb_new, alb_ice, pctsrf_new)
+
+
 
 ! Cette routine sert d'interface entre le modele atmospherique et un 
@@ -500,28 +561,227 @@
 !
 ! input:
+!   itime        numero du pas de temps
+!   iim, jjm     nbres de pts de grille
+!   dtime        pas de tempsde la physique
+!   klon         nombre total de points de grille
 !   nisurf       index de la surface a traiter (1 = sol continental)
+!   pctsrf       tableau des fractions de surface de chaque maille
+!   knon         nombre de points de la surface a traiter
+!   knindex      index des points de la surface a traiter
+!   rlon         longitudes
+!   rlat         latitudes
+!   debut        logical: 1er appel a la physique
+!   lafin        logical: dernier appel a la physique
 !   ocean        type d'ocean
+!   nexca        frequence de couplage
+!   swdown       flux solaire entrant a la surface
+!   lwdown       flux IR entrant a la surface
+!   precip_rain  precipitation liquide
+!   precip_snow  precipitation solide
+!   evap         evaporation
+!   tsurf        temperature de surface
+!   fder         derivee dF/dT
+!   albsol       albedo du sol (coherent avec swdown)
+!   taux         tension de vent en x
+!   tauy         tension de vent en y
+!   nexca        frequence de couplage
+!
 !
 ! output:
-!
+!   tsurf_new    temperature au sol
+!   alb_new      albedo
+!   pctsrf_new   nouvelle repartition des surfaces
+!   alb_ice      albedo de la glace
+!
+
+#include 'indicesol.h'
 
 ! Parametres d'entree
+  integer, intent(IN) :: itime
+  integer, intent(IN) :: iim, jjm
+  real, intent(IN) :: dtime
+  integer, intent(IN) :: klon
   integer, intent(IN) :: nisurf
+  integer, intent(IN) :: knon
+  real, dimension(klon,nbsrf), intent(IN) :: pctsrf
+  integer, dimension(knon), intent(in) :: knindex
+  logical, intent(IN) :: debut, lafin
+  real, dimension(klon), intent(IN) :: rlon, rlat
   character (len = 6)  :: ocean
+  real, dimension(knon), intent(IN) :: lwdown, swdown
+  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
+  real, dimension(knon), intent(IN) :: tsurf, fder, albsol, taux, tauy
+  integer              :: nexca
+
+  real, dimension(knon), intent(INOUT) :: evap
 
 ! Parametres de sortie
+  real, dimension(knon), intent(OUT):: tsurf_new, alb_new, alb_ice
+  real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
 
 ! Variables locales
-
-
+  integer                    :: j, error, sum_error, ig
+  integer                    :: npas
+  character (len = 20) :: modname = 'interfoce_cpl'
+  character (len = 80) :: abort_message
+  logical              :: check = .true.
+! variables pour moyenner les variables de couplage
+  real, allocatable, dimension(:),save :: cpl_sols, cpl_nsol, cpl_rain
+  real, allocatable, dimension(:),save :: cpl_snow, cpl_evap, cpl_tsol
+  real, allocatable, dimension(:),save :: cpl_fder, cpl_albe, cpl_taux
+  real, allocatable, dimension(:),save :: cpl_tauy, cpl_ruis
+! variables a passer au coupleur
+  real, dimension(iim, jjm+1) :: wri_sols, wri_nsol, wri_rain
+  real, dimension(iim, jjm+1) :: wri_snow, wri_evap, wri_tsol
+  real, dimension(iim, jjm+1) :: wri_fder, wri_albe, wri_taux
+  real, dimension(iim, jjm+1) :: wri_tauy, wri_ruis
+! variables relues par le coupleur
+  real, dimension(iim, jjm+1) :: read_sst, read_sic
+  real, dimension(iim, jjm+1) :: read_alb_sst, read_alb_sic
+! variable tampon
+  real, dimension(klon)       :: tamp
+  real, dimension(knon)       :: tamp_sic
+
+
+!
 ! Initialisation
+!
+  if (debut) then
+    sum_error = 0
+    allocate(cpl_sols(knon), stat = error)
+    sum_error = sum_error + error
+    allocate(cpl_nsol(knon), stat = error)
+    sum_error = sum_error + error
+    allocate(cpl_rain(knon), stat = error)
+    sum_error = sum_error + error
+    allocate(cpl_snow(knon), stat = error)
+    sum_error = sum_error + error
+    allocate(cpl_evap(knon), stat = error)
+    sum_error = sum_error + error
+    allocate(cpl_tsol(knon), stat = error)
+    sum_error = sum_error + error
+    allocate(cpl_fder(knon), stat = error)
+    sum_error = sum_error + error
+    allocate(cpl_albe(knon), stat = error)
+    sum_error = sum_error + error
+    allocate(cpl_taux(knon), stat = error)
+    sum_error = sum_error + error
+    allocate(cpl_tauy(knon), stat = error)
+    sum_error = sum_error + error
+    allocate(cpl_ruis(knon), stat = error)
+    sum_error = sum_error + error
+    if (sum_error /= 0) then
+      abort_message='Pb allocation variables couplees'
+      call abort_gcm(modname,abort_message,1)
+    endif
+    cpl_sols = 0.
+    cpl_nsol = 0.
+    cpl_rain = 0.
+    cpl_snow = 0.
+    cpl_evap = 0.
+    cpl_tsol = 0.
+    cpl_fder = 0.
+    cpl_albe = 0.
+    cpl_taux = 0.
+    cpl_tauy = 0.
+    cpl_ruis = 0.
+!
+! initialisation couplage
+!
+    call inicma(npas, nexca, dtime)
+!
+! 1ere lecture champs ocean
+!
+    call fromcpl(itime,(jjm+1)*iim,                                          &
+     &        read_sst, read_sic, read_alb_sst, read_alb_sic)
+    call cpl2gath(read_sst, tsurf_new, klon, knon,iim,jjm, knindex)
+    call cpl2gath(read_sic, tamp_sic , klon, knon,iim,jjm, knindex)
+    call cpl2gath(read_alb_sst, alb_new, klon, knon,iim,jjm, knindex)
+    call cpl2gath(read_alb_sic, alb_ice, klon, knon,iim,jjm, knindex)
+!
+! transformer tamp_sic en pctsrf_new
+!
+    do ig = 1, klon
+      IF (pctsrf(ig,is_oce) > epsfra .OR.            &
+     &             pctsrf(ig,is_sic) > epsfra) THEN
+            pctsrf_new(ig,is_oce) = pctsrf(ig,is_oce)    &
+     &                        - (tamp_sic(ig)-pctsrf(ig,is_sic))
+            pctsrf_new(ig,is_sic) = tamp_sic(ig)
+      endif
+    enddo
+
+  endif ! fin if (debut)
+
 ! fichier restart et fichiers histoires
 
 ! calcul des fluxs a passer
 
+  cpl_sols = cpl_sols + swdown      / FLOAT(nexca)
+  cpl_nsol = cpl_nsol + lwdown      / FLOAT(nexca)
+  cpl_rain = cpl_rain + precip_rain / FLOAT(nexca)
+  cpl_snow = cpl_snow + precip_snow / FLOAT(nexca)
+  cpl_evap = cpl_evap + evap        / FLOAT(nexca)
+  cpl_tsol = cpl_tsol + tsurf       / FLOAT(nexca)
+  cpl_fder = cpl_fder + fder        / FLOAT(nexca)
+  cpl_albe = cpl_albe + albsol      / FLOAT(nexca)
+  cpl_taux = cpl_taux + taux        / FLOAT(nexca)
+  cpl_tauy = cpl_tauy + tauy        / FLOAT(nexca)
+  cpl_ruis = cpl_ruis + run_off     / FLOAT(nexca)/dtime
+
+  if (mod(itime, nexca) == 0) then
+!
+! Mise sur la bonne grille des champs a passer au coupleur
+    call gath2cpl(cpl_sols, wri_sols, klon, knon,iim,jjm, knindex)
+    call gath2cpl(cpl_nsol, wri_nsol, klon, knon,iim,jjm, knindex)
+    call gath2cpl(cpl_rain, wri_rain, klon, knon,iim,jjm, knindex)
+    call gath2cpl(cpl_snow, wri_snow, klon, knon,iim,jjm, knindex)
+    call gath2cpl(cpl_evap, wri_evap, klon, knon,iim,jjm, knindex)
+    call gath2cpl(cpl_tsol, wri_tsol, klon, knon,iim,jjm, knindex)
+    call gath2cpl(cpl_fder, wri_fder, klon, knon,iim,jjm, knindex)
+    call gath2cpl(cpl_albe, wri_albe, klon, knon,iim,jjm, knindex)
+    call gath2cpl(cpl_taux, wri_taux, klon, knon,iim,jjm, knindex)
+    call gath2cpl(cpl_tauy, wri_tauy, klon, knon,iim,jjm, knindex)
+    call gath2cpl(cpl_ruis, wri_ruis, klon, knon,iim,jjm, knindex)
+!
+! Passage des champs au coupleur
+!
+    call intocpl(itime, iim, jjm , wri_sols, wri_nsol, wri_rain, wri_snow, &
+      & wri_evap, wri_tsol, wri_fder, wri_albe, wri_taux, wri_tauy,       &
+      & wri_ruis ) 
+    cpl_sols = 0.
+    cpl_nsol = 0.
+    cpl_rain = 0.
+    cpl_snow = 0.
+    cpl_evap = 0.
+    cpl_tsol = 0.
+    cpl_fder = 0.
+    cpl_albe = 0.
+    cpl_taux = 0.
+    cpl_tauy = 0.
+    cpl_ruis = 0.
+
+    call fromcpl(itime,(jjm+1)*iim,                                          &
+     &        read_sst, read_sic, read_alb_sst, read_alb_sic)
+    call cpl2gath(read_sst, tsurf_new, klon, knon,iim,jjm, knindex)
+    call cpl2gath(read_sic, tamp_sic , klon, knon,iim,jjm, knindex)
+    call cpl2gath(read_alb_sst, alb_new, klon, knon,iim,jjm, knindex)
+    call cpl2gath(read_alb_sic, alb_ice, klon, knon,iim,jjm, knindex)
+! transformer tamp_sic en pctsrf_new
+
+    do ig = 1, klon
+      IF (pctsrf(ig,is_oce) > epsfra .OR.            &
+     &             pctsrf(ig,is_sic) > epsfra) THEN
+            pctsrf_new(ig,is_oce) = pctsrf(ig,is_oce)    &
+     &                        - (tamp_sic(ig)-pctsrf(ig,is_sic))
+            pctsrf_new(ig,is_sic) = tamp_sic(ig)
+      endif
+    enddo
+  endif
+
   END SUBROUTINE interfoce_cpl
 !
 !#########################################################################
 !
+
   SUBROUTINE interfoce_slab(nisurf)
 
@@ -825,9 +1085,10 @@
 !
 
-  SUBROUTINE calcul_fluxs( knon, dtime, &
+  SUBROUTINE calcul_fluxs( knon, nisurf, dtime, &
      & tsurf, p1lay, cal, beta, coef1lay, ps, &
+     & precip_rain, precip_snow, snow, qsol, &
      & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
      & petAcoef, peqAcoef, petBcoef, peqBcoef, &
-     & tsurf_new, fluxlat, fluxsens, dflux_s, dflux_l)
+     & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
 
 ! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
@@ -838,4 +1099,5 @@
 ! input:
 !   knon         nombre de points a traiter
+!   nisurf       surface a traiter
 !   tsurf        temperature de surface
 !   p1lay        pression 1er niveau (milieu de couche)
@@ -844,4 +1106,9 @@
 !   coef1lay     coefficient d'echange
 !   ps           pression au sol
+!   precip_rain  precipitations liquides
+!   precip_snow  precipitations solides
+!   snow         champs hauteur de neige
+!   qsol         humidite du sol
+!   runoff       runoff en cas de trop plein
 !   petAcoef     coeff. A de la resolution de la CL pour t
 !   peqAcoef     coeff. A de la resolution de la CL pour q
@@ -862,7 +1129,8 @@
 #include "YOETHF.inc"
 #include "FCTTRE.inc"
+#include 'indicesol.h'
 
 ! Parametres d'entree
-  integer, intent(IN) :: knon
+  integer, intent(IN) :: knon, nisurf
   real   , intent(IN) :: dtime
   real, dimension(knon), intent(IN) :: petAcoef, peqAcoef
@@ -870,9 +1138,11 @@
   real, dimension(knon), intent(IN) :: ps, q1lay
   real, dimension(knon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
+  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
   real, dimension(knon), intent(IN) :: radsol, dif_grnd
   real, dimension(knon), intent(IN) :: t1lay, u1lay, v1lay
+  real, dimension(knon), intent(INOUT) :: snow, qsol
 
 ! Parametres sorties
-  real, dimension(knon), intent(OUT):: tsurf_new, fluxlat, fluxsens
+  real, dimension(knon), intent(OUT):: tsurf_new, evap, fluxsens, fluxlat
   real, dimension(knon), intent(OUT):: dflux_s, dflux_l
 
@@ -885,9 +1155,32 @@
   real, dimension(knon) :: zx_h_ts, zx_q_0 , d_ts
   real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
+  real                  :: bilan_f, fq_fonte
   real, parameter :: t_grnd = 271.35, t_coup = 273.15
   logical         :: check = .true.
   character (len = 20)  :: modname = 'calcul_fluxs'
+  logical         :: fonte_neige = .false.
+  real            :: max_eau_sol = 150.0
+  character (len = 80) :: abort_message 
 
   if (check) write(*,*)'Entree ', modname
+
+  if (size(run_off) /= knon .AND. nisurf == is_ter) then
+    write(*,*)'Bizarre, le nombre de points continentaux'
+    write(*,*)'a change entre deux appels. J''arrete ...'
+    abort_message='Pb run_off'
+    call abort_gcm(modname,abort_message,1)
+  endif
+!
+! Traitement neige et humidite du sol
+!
+    if (nisurf == is_oce) then
+      snow = 0.
+      qsol = max_eau_sol
+    else
+      snow = max(0.0, snow + (precip_snow - evap) * dtime)
+      qsol = qsol + (precip_rain - evap) * dtime
+    endif 
+
+
 ! 
 ! Initialisation
@@ -961,4 +1254,12 @@
      &                       zx_nh(i) + zx_sl(i) * zx_nq(i)) &  
      &                     + dtime * dif_grnd(i))
+
+!
+! Y'a-t-il fonte de neige?
+!
+    fonte_neige = (nisurf /= is_oce) .AND. &
+     & (snow(i) > 0. .OR. nisurf == is_sic .OR. nisurf == is_lic) &
+     & .AND. (tsurf_new(i) >= RTT)
+    if (fonte_neige) tsurf_new(i) = RTT  
     zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
     d_ts(i) = tsurf_new(i) - tsurf(i)
@@ -966,9 +1267,25 @@
 !== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
 !== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
-    fluxlat(i) = zx_mq(i) + zx_nq(i) * tsurf_new(i) 
+    evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i) 
+    fluxlat(i) = - evap(i) * zx_sl(i)
     fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
 ! Derives des flux dF/dTs (W m-2 K-1):
     dflux_s(i) = zx_nh(i)
     dflux_l(i) = (zx_sl(i) * zx_nq(i))
+!
+! en cas de fonte de neige
+!
+    if (fonte_neige) then
+      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
+     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
+     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
+      bilan_f = max(0., bilan_f)
+      fq_fonte = bilan_f / zx_sl(i)
+      snow(i) = max(0., snow(i) - fq_fonte * dtime)
+      qsol(i) = qsol(i) + (fq_fonte * dtime)
+    endif
+    if (nisurf == is_ter)  &
+     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
+    qsol(i) = min(qsol(i), max_eau_sol) 
   ENDDO
 
@@ -978,3 +1295,201 @@
 !
 
+  SUBROUTINE sol_dem_write(itime, klon, rlon, rlat, &
+  &                     pctsrf_new,tsurf_new,alb_new)
+
+! Routine d'ecriture de l'etat de redemarrage pour le sol
+!
+! L.Fairhead
+!
+! input:
+!   itime        numero du pas de temps
+!   klon         nombre total de points de grille
+!   rlon         longitudes
+!   rlat         latitudes
+!   tsurf_new    temperature au sol
+!   alb_new      albedo
+!   pctsrf_new   repartition des surfaces
+
+  include 'indicesol.h'
+#include 'temps.inc'
+  include 'netcdf.inc'
+
+! Parametres d'entree
+  integer, intent(IN) :: itime
+  integer, intent(IN) :: klon
+  real, dimension(klon), intent(IN) :: rlon, rlat
+  real, dimension(klon,nbsrf), intent(IN)  :: tsurf_new, alb_new
+  real, dimension(klon,nbsrf), intent(IN) :: pctsrf_new
+
+! Variables locales
+  integer             :: ierr, nid
+  integer             :: idim1, idim2, idim3
+  integer,parameter   :: length = 100
+  character (len = 20) :: modname = 'sol_dem_write'
+  character (len = 80) :: abort_message 
+  real, dimension(length) :: tab_cntrl = 0.
+  integer                 :: nvarid
+
+  ierr = NF_CREATE('restartsol', NF_CLOBBER, nid)
+  IF (ierr.NE.NF_NOERR) THEN
+    abort_message=' Pb d''ouverture du fichier restartsol'
+    CALL abort_gcm(modname,abort_message,ierr)
+  ENDIF
+
+  ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 23,    &
+     &                    "Fichier redemmarage sol")
+  ierr = NF_DEF_DIM (nid, "index", length, idim1)
+  ierr = NF_DEF_DIM (nid, "points_physiques", klon, idim2)
+  ierr = NF_DEF_DIM (nid, "nombre_surfaces", nbsrf, idim3)
+  ierr = NF_ENDDEF(nid)
+
+  tab_cntrl(13) = day_end
+  tab_cntrl(14) = anne_ini
+
+  ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+  ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
+#else
+  ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
+#endif
+  ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 22,  &
+     &                        "Parametres de controle")
+  ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
+#else
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
+#endif
+
+  ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+  ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
+#else
+  ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
+#endif
+  ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 32,  &
+     &               "Longitudes de la grille physique")
+  ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlon)
+#else
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,rlon)
+#endif
+!
+  ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+  ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
+#else
+  ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
+#endif
+  ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 31,    &
+     &                        "Latitudes de la grille physique")
+  ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlat)
+#else
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,rlat)
+#endif
+        ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+        ierr = NF_DEF_VAR (nid, "TS", NF_DOUBLE, 1, idim2,nvarid)
+#else
+        ierr = NF_DEF_VAR (nid, "TS", NF_FLOAT, 1, idim2,nvarid)
+#endif
+        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 22, &
+     &                        "Temperature de surface")
+        ierr = NF_ENDDEF(nid)
+
+
+
+
+  END SUBROUTINE sol_dem_write
+!
+!#########################################################################
+!
+  SUBROUTINE gath2cpl(champ_in, champ_out, klon, knon, iim, jjm, knindex)
+
+! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
+! au coupleur.
+!
+! 
+! input:         
+!   champ_in     champ sur la grille gathere        
+!   knon         nombre de points dans le domaine a traiter
+!   knindex      index des points de la surface a traiter
+!   klon         taille de la grille
+!   iim,jjm      dimension de la grille 2D
+!
+! output:
+!   champ_out    champ sur la grille 2D
+!
+! input
+  integer                   :: klon, knon, iim, jjm
+  real, dimension(knon)     :: champ_in
+  integer, dimension(knon)  :: knindex
+! output
+  real, dimension(iim,jjm+1)  :: champ_out
+! local
+  integer                   :: i, ig, j
+  real, dimension(klon)     :: tamp
+
+  do i = 1, knon
+    ig = knindex(i)
+    tamp(ig) = champ_in(i)
+  enddo    
+  champ_out(:,1) = tamp(1)
+  do j = 2, jjm
+    do i = 1, iim
+      champ_out(i,j) = tamp((j-2)*jjm + i + 1)
+    enddo
+  enddo
+  champ_out(:,jjm+1) = tamp(klon)
+
+  END SUBROUTINE gath2cpl
+!
+!#########################################################################
+!
+  SUBROUTINE cpl2gath(champ_in, champ_out, klon, knon, iim, jjm, knindex)
+
+! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
+! au coupleur.
+!
+! 
+! input:         
+!   champ_in     champ sur la grille gathere        
+!   knon         nombre de points dans le domaine a traiter
+!   knindex      index des points de la surface a traiter
+!   klon         taille de la grille
+!   iim,jjm      dimension de la grille 2D
+!
+! output:
+!   champ_out    champ sur la grille 2D
+!
+! input
+  integer                   :: klon, knon, iim, jjm
+  real, dimension(iim,jjm+1)     :: champ_in
+  integer, dimension(knon)  :: knindex
+! output
+  real, dimension(knon)  :: champ_out
+! local
+  integer                   :: i, ig, j
+  real, dimension(klon)     :: tamp
+
+  tamp(1) = champ_in(1,1)
+  do j = 2, jjm
+    do i = 1, iim
+      tamp((j-2)*jjm + i + 1) = champ_in(i,j)
+    enddo
+  enddo
+  tamp(klon) = champ_in(1,jjm+1)
+
+  do i = 1, knon
+    ig = knindex(i)
+    champ_out(i) = tamp(ig)
+  enddo    
+
+  END SUBROUTINE cpl2gath
+!
+!#########################################################################
+!
   END MODULE interface_surf
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F	(revision 98)
@@ -386,5 +386,5 @@
       END
 
-      SUBROUTINE fromcpl(kt, imjm, sst, gla)
+      SUBROUTINE fromcpl(kt, imjm, sst,sic, alb_sst, alb_sic )
       IMPLICIT none
 c
@@ -399,5 +399,8 @@
       INTEGER imjm, kt
       REAL sst(imjm)          ! sea-surface-temperature
-      REAL gla(imjm)          ! sea-ice
+      REAL alb_sst(imjm)  ! open sea albedo
+      REAL sic(imjm)      ! sea ice cover
+      REAL alb_sic(imjm)  ! sea ice albedo
+
 c
       INTEGER nuout             ! listing output unit
@@ -447,5 +450,11 @@
      $          nuout)
             IF (jf.eq.2)
-     $          CALL locread(cl_read(jf), gla, imjm, nuread, iflag,
+     $          CALL locread(cl_read(jf), sic, imjm, nuread, iflag,
+     $          nuout)
+            IF (jf.eq.3)
+     $          CALL locread(cl_read(jf), alb_sst, imjm, nuread, iflag,
+     $          nuout)
+            IF (jf.eq.4)
+     $          CALL locread(cl_read(jf), alb_sic, imjm, nuread, iflag,
      $          nuout)
             CLOSE (nuread)
@@ -474,5 +483,15 @@
 c
           CALL SIPC_Read_Model(index, imjm, cmodinf, 
-     $              cljobnam_r,infos, gla)
+     $              cljobnam_r,infos, sic)
+c         Index of open sea albedo in total number of fields jpfldo2a: 
+          index = 3
+c
+          CALL SIPC_Read_Model(index, imjm, cmodinf, 
+     $              cljobnam_r,infos, alb_sst)
+c         Index of sea-ice albedo in total number of fields jpfldo2a: 
+          index = 4
+c
+          CALL SIPC_Read_Model(index, imjm, cmodinf, 
+     $              cljobnam_r,infos, alb_sic)
 c
 c
@@ -482,8 +501,10 @@
 c exchanges from ocean=CPL to atmosphere
 c
-          DO jf=1,jpfldo2a
-            IF (jf.eq.1) CALL CLIM_Import (cl_read(jf) , kt, sst, info)
-            IF (jf.eq.2) CALL CLIM_Import (cl_read(jf) , kt, gla, info)
-            IF ( info .NE. CLIM_Ok) THEN
+        DO jf=1,jpfldo2a
+          IF (jf.eq.1) CALL CLIM_Import (cl_read(jf) , kt, sst, info)
+          IF (jf.eq.2) CALL CLIM_Import (cl_read(jf) , kt, sic, info)
+         IF (jf.eq.3) CALL CLIM_Import (cl_read(jf) , kt, alb_sst, info)
+         IF (jf.eq.4) CALL CLIM_Import (cl_read(jf) , kt, alb_sic, info)
+          IF ( info .NE. CLIM_Ok) THEN
                 WRITE(nuout,*)'Pb in reading ', cl_read(jf), jf
                 WRITE(nuout,*)'Couplage kt is = ',kt
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.h	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.h	(revision 98)
@@ -1,9 +1,14 @@
-      LOGICAL ok_oasis
-      PARAMETER (ok_oasis = .FALSE.)
-c
-      CHARACTER*4 cchan
-      PARAMETER (cchan="PIPE")
-c      PARAMETER (cchan="CLIM")
+C
+C -- oasis.h   
+C    ******
+C@
+C@  Contents : choice for the OASIS version: clim or pipe
+C@  --------
 
-      INTEGER jpmaxfld
-      PARAMETER(jpmaxfld = 20)
+      logical ok_oasis
+      parameter(ok_oasis = .false.)
+
+      CHARACTER*8 cchan
+      PARAMETER ( cchan='PIPE' )
+C
+C     --- end of oasis.h
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/param_cou.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/param_cou.h	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/param_cou.h	(revision 98)
@@ -2,10 +2,13 @@
 C -- param_cou.h
 C
+	INTEGER jpmaxfld
+	PARAMETER(jpmaxfld = 100)        ! Number of maximum fields
+                                         ! exchange betwwen ocean and atmosphere
 	INTEGER jpflda2o
-	PARAMETER(jpflda2o = 1)          ! Number of fields exchanged from
+	PARAMETER(jpflda2o = 8)          ! Number of fields exchanged from
                                          ! atmosphere to ocean
 C
 	INTEGER jpfldo2a
-	PARAMETER(jpfldo2a = 1)          ! Number of fields exchanged from
+	PARAMETER(jpfldo2a = 2)          ! Number of fields exchanged from
                                          ! ocean to atmosphere
 C
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/param_sipc.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/param_sipc.h	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/param_sipc.h	(revision 98)
@@ -10,5 +10,5 @@
 	INTEGER jpbyteint,jpbyterea, jpbytecha  
         PARAMETER (jpbyteint = 4)
-	PARAMETER (jpbyterea = 4)
+	PARAMETER (jpbyterea = 8)
 	PARAMETER (jpbytecha = 1)
 C@
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/phyetat0.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/phyetat0.F	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/phyetat0.F	(revision 98)
@@ -1,4 +1,5 @@
       SUBROUTINE phyetat0 (fichnom,dtime,co2_ppm,solaire,
-     .            rlat,rlon,tsol,tsoil,deltat,qsol,snow,
+     .            rlat,rlon, pctsrf, tsol,tsoil,deltat,qsol,snow,
+     .           albe, evap, rain_fall, snow_fall, solsw, sollw,
      .           radsol,rugmer,agesno,clesphy0,
      .           zmea,zstd,zsig,zgam,zthe,zpic,zval,rugsrel,tabcntr0,
@@ -27,5 +28,11 @@
       REAL qsol(klon,nbsrf)
       REAL snow(klon,nbsrf)
+      REAL albe(klon,nbsrf)
+      REAL evap(klon,nbsrf)
       REAL radsol(klon)
+      REAL rain_fall(klon)
+      REAL snow_fall(klon)
+      REAL sollw(klon)
+      real solsw(klon)
       REAL rugmer(klon)
       REAL agesno(klon)
@@ -38,4 +45,6 @@
       REAL zval(klon)
       REAL rugsrel(klon)
+      REAL pctsrf(klon, nbsrf)
+      REAL fractint(klon)
 
       REAL t_ancien(klon,klev), q_ancien(klon,klev)
@@ -58,4 +67,5 @@
 c Ouvrir le fichier contenant l'etat initial:
 c
+      print*,'fichnom',fichnom
       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
       IF (ierr.NE.NF_NOERR) THEN
@@ -184,5 +194,123 @@
          CALL abort
       ENDIF
-c
+C
+C
+C Lecture du masque terre mer
+C
+      ierr = NF_INQ_VARID (nid, "masque", nvarid)
+      IF (ierr .EQ.  NF_NOERR) THEN
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zmasq)
+#else
+          ierr = NF_GET_VAR_REAL(nid, nvarid, zmasq)
+#endif
+          IF (ierr.NE.NF_NOERR) THEN
+              PRINT*, 'phyetat0: Lecture echouee pour <masque>'
+              CALL abort
+          ENDIF
+      else
+          PRINT*, 'phyetat0: Le champ <masque> est absent'
+          PRINT*, 'fichier startphy non compatible avec phyetat0'
+C      CALL abort
+      ENDIF
+C Lecture des fractions pour chaque sous-surface
+C
+C initialisation des sous-surfaces
+C
+      pctsrf = 0.
+C
+C fraction de terre
+C
+      ierr = NF_INQ_VARID (nid, "FTER", nvarid)
+      IF (ierr .EQ.  NF_NOERR) THEN
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, pctsrf(1 : klon,is_ter))
+#else
+          ierr = NF_GET_VAR_REAL(nid, nvarid, pctsrf(1 : klon,is_ter))
+#endif
+          IF (ierr.NE.NF_NOERR) THEN
+              PRINT*, 'phyetat0: Lecture echouee pour <FTER>'
+              CALL abort
+          ENDIF
+      else
+          PRINT*, 'phyetat0: Le champ <FTER> est absent'
+c$$$         CALL abort
+      ENDIF
+C
+C fraction de glace de terre
+C
+      ierr = NF_INQ_VARID (nid, "FLIC", nvarid)
+      IF (ierr .EQ.  NF_NOERR) THEN
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, pctsrf(1 : klon,is_lic))
+#else
+          ierr = NF_GET_VAR_REAL(nid, nvarid, pctsrf(1 : klon,is_lic))
+#endif
+          IF (ierr.NE.NF_NOERR) THEN
+              PRINT*, 'phyetat0: Lecture echouee pour <FLIC>'
+              CALL abort
+          ENDIF
+      else
+          PRINT*, 'phyetat0: Le champ <FLIC> est absent'
+c$$$         CALL abort
+      ENDIF
+C
+C fraction d'ocean
+C
+      ierr = NF_INQ_VARID (nid, "FOCE", nvarid)
+      IF (ierr .EQ.  NF_NOERR) THEN
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, pctsrf(1 : klon,is_oce))
+#else
+          ierr = NF_GET_VAR_REAL(nid, nvarid, pctsrf(1 : klon,is_oce))
+#endif
+          IF (ierr.NE.NF_NOERR) THEN
+              PRINT*, 'phyetat0: Lecture echouee pour <FOCE>'
+              CALL abort
+          ENDIF
+      else
+          PRINT*, 'phyetat0: Le champ <FOCE> est absent'
+c$$$         CALL abort
+      ENDIF
+C
+C fraction glace de mer
+C
+      ierr = NF_INQ_VARID (nid, "FSIC", nvarid)
+      IF (ierr .EQ.  NF_NOERR) THEN
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, pctsrf(1 : klon,is_sic))
+#else
+          ierr = NF_GET_VAR_REAL(nid, nvarid, pctsrf(1 : klon, is_sic))
+#endif
+          IF (ierr.NE.NF_NOERR) THEN
+              PRINT*, 'phyetat0: Lecture echouee pour <FSIC>'
+              CALL abort
+          ENDIF
+      else
+          PRINT*, 'phyetat0: Le champ <FSIC> est absent'
+c$$$         CALL abort
+      ENDIF
+C
+C  Verification de l'adequation entre le masque et les sous-surfaces
+C
+      fractint( 1 : klon) = pctsrf(1 : klon, is_ter) 
+     $    + pctsrf(1 : klon, is_lic)
+      DO i = 1 , klon
+        IF ( abs(fractint(i) - zmasq(i) ) .GT. EPSFRA ) THEN
+            WRITE(*,*) 'phyetat0: attention fraction terre pas ', 
+     $          'coherente ', i, zmasq(i), pctsrf(i, is_ter)
+     $          ,pctsrf(i, is_lic)
+        ENDIF 
+      END DO 
+      fractint (1 : klon) =  pctsrf(1 : klon, is_oce) 
+     $    + pctsrf(1 : klon, is_sic)
+      DO i = 1 , klon
+        IF ( abs( fractint(i) - (1. - zmasq(i))) .GT. EPSFRA ) THEN
+            WRITE(*,*) 'phyetat0 attention fraction ocean pas ', 
+     $          'coherente ', i, zmasq(i) , pctsrf(i, is_oce)
+     $          ,pctsrf(i, is_sic)
+        ENDIF 
+      END DO 
+C
 c Lecture des temperatures du sol:
 c
@@ -419,4 +547,226 @@
       ENDIF
 c
+c Lecture de albedo au sol:
+c
+      ierr = NF_INQ_VARID (nid, "ALBE", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <ALBE> est absent'
+         PRINT*, '          Mais je vais essayer de lire ALBE**'
+         DO nsrf = 1, nbsrf
+           IF (nsrf.GT.99) THEN
+             PRINT*, "Trop de sous-mailles"
+             CALL abort
+           ENDIF
+           WRITE(str2,'(i2.2)') nsrf
+           ierr = NF_INQ_VARID (nid, "ALBE"//str2, nvarid)
+           IF (ierr.NE.NF_NOERR) THEN
+              PRINT*, "phyetat0: Le champ <ALBE"//str2//"> est absent"
+              CALL abort
+           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))
+#endif
+           IF (ierr.NE.NF_NOERR) THEN
+             PRINT*, "phyetat0: Lecture echouee pour <ALBE"//str2//">"
+             CALL abort
+           ENDIF
+           xmin = 1.0E+20
+           xmax = -1.0E+20
+           DO i = 1, klon
+              xmin = MIN(snow(i,nsrf),xmin)
+              xmax = MAX(snow(i,nsrf),xmax)
+           ENDDO
+           PRINT*,'Neige du sol ALBE**:', nsrf, xmin, xmax
+         ENDDO
+      ELSE
+         PRINT*, 'phyetat0: Le champ <ALBE> est present'
+         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))
+#endif
+         IF (ierr.NE.NF_NOERR) THEN
+            PRINT*, "phyetat0: Lecture echouee pour <ALBE>"
+            CALL abort
+         ENDIF
+         xmin = 1.0E+20
+         xmax = -1.0E+20
+         DO i = 1, klon
+            xmin = MIN(albe(i,1),xmin)
+            xmax = MAX(albe(i,1),xmax)
+         ENDDO
+         PRINT*,'Neige du sol <ALBE>', xmin, xmax
+         DO nsrf = 2, nbsrf
+         DO i = 1, klon
+            albe(i,nsrf) = albe(i,1)
+         ENDDO
+         ENDDO
+      ENDIF
+
+c
+c Lecture de evaporation:  
+c
+      ierr = NF_INQ_VARID (nid, "EVAP", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <EVAP> est absent'
+         PRINT*, '          Mais je vais essayer de lire EVAP**'
+         DO nsrf = 1, nbsrf
+           IF (nsrf.GT.99) THEN
+             PRINT*, "Trop de sous-mailles"
+             CALL abort
+           ENDIF
+           WRITE(str2,'(i2.2)') nsrf
+           ierr = NF_INQ_VARID (nid, "EVAP"//str2, nvarid)
+           IF (ierr.NE.NF_NOERR) THEN
+              PRINT*, "phyetat0: Le champ <EVAP"//str2//"> est absent"
+              CALL abort
+           ENDIF
+#ifdef NC_DOUBLE
+           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, evap(1,nsrf))
+#else
+           ierr = NF_GET_VAR_REAL(nid, nvarid, evap(1,nsrf))
+#endif
+           IF (ierr.NE.NF_NOERR) THEN
+             PRINT*, "phyetat0: Lecture echouee pour <EVAP"//str2//">"
+             CALL abort
+           ENDIF
+           xmin = 1.0E+20
+           xmax = -1.0E+20
+           DO i = 1, klon
+              xmin = MIN(evap(i,nsrf),xmin)
+              xmax = MAX(evap(i,nsrf),xmax)
+           ENDDO
+           PRINT*,'Neige du sol EVAP**:', nsrf, xmin, xmax
+         ENDDO
+      ELSE
+         PRINT*, 'phyetat0: Le champ <EVAP> est present'
+         PRINT*, '          J ignore donc les autres EVAP**'
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, evap(1,1))
+#else
+         ierr = NF_GET_VAR_REAL(nid, nvarid, evap(1,1))
+#endif
+         IF (ierr.NE.NF_NOERR) THEN
+            PRINT*, "phyetat0: Lecture echouee pour <EVAP>"
+            CALL abort
+         ENDIF
+         xmin = 1.0E+20
+         xmax = -1.0E+20
+         DO i = 1, klon
+            xmin = MIN(evap(i,1),xmin)
+            xmax = MAX(evap(i,1),xmax)
+         ENDDO
+         PRINT*,'Neige du sol <EVAP>', xmin, xmax
+         DO nsrf = 2, nbsrf
+         DO i = 1, klon
+            evap(i,nsrf) = evap(i,1)
+         ENDDO
+         ENDDO
+      ENDIF
+c
+c Lecture precipitation liquide:
+c
+      ierr = NF_INQ_VARID (nid, "rain_f", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <rain_f> est absent'
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rain_fall)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, rain_fall)
+#endif
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Lecture echouee pour <rain_f>'
+         CALL abort
+      ENDIF
+      xmin = 1.0E+20
+      xmax = -1.0E+20
+      DO i = 1, klon
+         xmin = MIN(rain_fall(i),xmin)
+         xmax = MAX(rain_fall(i),xmax)
+      ENDDO
+      PRINT*,'Precipitation liquide rain_f:', xmin, xmax
+c
+c Lecture precipitation solide:
+c
+      ierr = NF_INQ_VARID (nid, "snow_f", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <snow_f> est absent'
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, snow_fall)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, snow_fall)
+#endif
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Lecture echouee pour <snow_f>'
+         CALL abort
+      ENDIF
+      xmin = 1.0E+20
+      xmax = -1.0E+20
+      DO i = 1, klon
+         xmin = MIN(snow_fall(i),xmin)
+         xmax = MAX(snow_fall(i),xmax)
+      ENDDO
+      PRINT*,'Precipitation solide snow_f:', xmin, xmax
+c
+c Lecture rayonnement solaire au sol:
+c
+      ierr = NF_INQ_VARID (nid, "solsw", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <solsw> est absent'
+         PRINT*, 'mis a zero'
+         solsw = 0.
+      ELSE
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, solsw)
+#else
+        ierr = NF_GET_VAR_REAL(nid, nvarid, solsw)
+#endif
+        IF (ierr.NE.NF_NOERR) THEN
+          PRINT*, 'phyetat0: Lecture echouee pour <solsw>'
+          CALL abort
+        ENDIF
+      ENDIF
+      xmin = 1.0E+20
+      xmax = -1.0E+20
+      DO i = 1, klon
+         xmin = MIN(solsw(i),xmin)
+         xmax = MAX(solsw(i),xmax)
+      ENDDO
+      PRINT*,'Rayonnement solaire au sol solsw:', xmin, xmax
+c
+c Lecture rayonnement IF au sol:
+c
+      ierr = NF_INQ_VARID (nid, "sollw", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <sollw> est absent'
+         PRINT*, 'mis a zero'
+         sollw = 0.
+      ELSE
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, sollw)
+#else
+        ierr = NF_GET_VAR_REAL(nid, nvarid, sollw)
+#endif
+        IF (ierr.NE.NF_NOERR) THEN
+          PRINT*, 'phyetat0: Lecture echouee pour <sollw>'
+          CALL abort
+        ENDIF
+      ENDIF
+      xmin = 1.0E+20
+      xmax = -1.0E+20
+      DO i = 1, klon
+         xmin = MIN(sollw(i),xmin)
+         xmax = MAX(sollw(i),xmax)
+      ENDDO
+      PRINT*,'Rayonnement IF au sol sollw:', xmin, xmax
+
+c
 c Lecture du rayonnement net au sol:
 c
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F	(revision 98)
@@ -1,4 +1,6 @@
       SUBROUTINE phyredem (fichnom,dtime,radpas,co2_ppm,solaire,
-     .           rlat,rlon,tsol,tsoil,deltat,qsol,snow,
+     .           rlat,rlon, pctsrf,tsol,tsoil,deltat,qsol,snow,
+     .           albedo, evap, rain_fall, snow_fall,
+     .           solsw, sollw,
      .           radsol,rugmer,agesno,
      .           zmea,zstd,zsig,zgam,zthe,zpic,zval,rugsrel,
@@ -29,4 +31,10 @@
       REAL qsol(klon,nbsrf)
       REAL snow(klon,nbsrf)
+      REAL albedo(klon,nbsrf)
+      REAL evap(klon,nbsrf)
+      REAL rain_fall(klon)
+      REAL snow_fall(klon)
+      real solsw(klon)
+      real sollw(klon)
       REAL radsol(klon)
       REAL rugmer(klon)
@@ -40,4 +48,5 @@
       REAL zval(klon)
       REAL rugsrel(klon)
+      REAL pctsrf(klon, nbsrf)
       REAL t_ancien(klon,klev), q_ancien(klon,klev)
 c
@@ -134,4 +143,91 @@
 #endif
 c
+C PB ajout du masque terre/mer
+C
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "masque", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "masque", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 16,
+     .                        "masque terre mer")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zmasq)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,zmasq)
+#endif      
+c BP ajout des fraction de chaque sous-surface
+C
+C 1. fraction de terre 
+C
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "FTER", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
+     .                        "fraction de continent")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pctsrf(1 : klon, is_ter))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,pctsrf(1 : klon, is_ter))
+#endif
+C 
+C 2. Fraction de glace de terre
+C 
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "FLIC", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 24,
+     .                        "fraction glace de terre")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pctsrf(1 : klon,is_lic))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,pctsrf(1 : klon, is_lic))
+#endif
+C
+C 3. fraction ocean
+C
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "FOCE", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 14,
+     .                        "fraction ocean")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pctsrf(1 : klon, is_oce))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,pctsrf(1 : klon, is_oce))
+#endif
+C
+C 4. Fraction glace de mer
+C
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "FSIC", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 18,
+     .                        "fraction glace mer")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pctsrf(1 : klon, is_sic))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,pctsrf(1 : klon, is_sic))
+#endif
+C
+C
 c
       DO nsrf = 1, nbsrf
@@ -227,4 +323,51 @@
         ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
+        ierr = NF_DEF_VAR (nid,"ALBE"//str2,NF_DOUBLE,1,idim2,nvarid)
+#else
+        ierr = NF_DEF_VAR (nid,"ALBE"//str2,NF_FLOAT,1,idim2,nvarid)
+#endif
+        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 23,
+     .                        "albedo de surface No."//str2)
+        ierr = NF_ENDDEF(nid)
+        ELSE
+        PRINT*, "Trop de sous-mailles"
+        CALL abort
+        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
+c
+      DO nsrf = 1, nbsrf
+        IF (nsrf.LE.99) THEN
+        WRITE(str2,'(i2.2)') nsrf
+        ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+        ierr = NF_DEF_VAR (nid,"EVAP"//str2,NF_DOUBLE,1,idim2,nvarid)
+#else
+        ierr = NF_DEF_VAR (nid,"EVAP"//str2,NF_FLOAT,1,idim2,nvarid)
+#endif
+        ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
+     .                        "Evaporation de surface No."//str2)
+        ierr = NF_ENDDEF(nid)
+        ELSE
+        PRINT*, "Trop de sous-mailles"
+        CALL abort
+        ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,evap(1,nsrf))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,evap(1,nsrf))
+#endif
+      ENDDO
+
+c
+      DO nsrf = 1, nbsrf
+        IF (nsrf.LE.99) THEN
+        WRITE(str2,'(i2.2)') nsrf
+        ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
         ierr = NF_DEF_VAR (nid,"SNOW"//str2,NF_DOUBLE,1,idim2,nvarid)
 #else
@@ -244,4 +387,5 @@
 #endif
       ENDDO
+
 c
       ierr = NF_REDEF (nid)
@@ -262,4 +406,64 @@
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "solsw", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "solsw", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 32,
+     .                        "Rayonnement solaire a la surface")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,solsw)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,solsw)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "sollw", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "sollw", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 27,
+     .                        "Rayonnement IF a la surface")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,sollw)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,sollw)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "rain_f", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "rain_f", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
+     .                        "precipitation liquide")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rain_fall)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rain_fall)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "snow_f", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "snow_f", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
+     .                        "precipitation solide")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,snow_fall)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,snow_fall)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
       ierr = NF_DEF_VAR (nid, "RUGMER", NF_DOUBLE, 1, idim2,nvarid)
 #else
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 97)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 98)
@@ -100,4 +100,6 @@
       REAL soilcap(klon,nbsrf), soilflux(klon,nbsrf)
       SAVE soilcap, soilflux
+      logical ok_veget
+      parameter (ok_veget = .false.)
 c======================================================================
 c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
@@ -132,9 +134,11 @@
       INTEGER iliq          ! indice de traceurs pour eau liquide
       PARAMETER (iliq=2)
-c
+
       INTEGER nvm           ! nombre de vegetations
       PARAMETER (nvm=8)
       REAL veget(klon,nvm)  ! couverture vegetale
       SAVE veget
+
+c
 c
 c Variables argument:
@@ -224,4 +228,7 @@
       SAVE ftsoil                 ! temperature dans le sol
 c
+      REAL fevap(klon,nbsrf)
+      SAVE fevap                 ! evaporation
+c
       REAL deltat(klon)
       SAVE deltat                 ! ecart avec la SST de reference
@@ -232,4 +239,7 @@
       REAL fsnow(klon,nbsrf)
       SAVE fsnow                  ! epaisseur neigeuse
+c
+      REAL falbe(klon,nbsrf)
+      SAVE falbe                  ! albedo par type de surface
 c
       REAL rugmer(klon)
@@ -379,9 +389,15 @@
       REAL cldemi(klon,klev)  ! emissivite infrarouge
 c
-      REAL fluxq(klon,klev)   ! flux turbulent d'humidite
-      REAL fluxt(klon,klev)   ! flux turbulent de chaleur
-      REAL fluxu(klon,klev)   ! flux turbulent de vitesse u
-      REAL fluxv(klon,klev)   ! flux turbulent de vitesse v
-c
+C§§§ PB 
+      REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
+      REAL fluxt(klon,klev, nbsrf)   ! flux turbulent de chaleur
+      REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
+      REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
+c
+      REAL zxfluxt(klon, klev)
+      REAL zxfluxq(klon, klev)
+      REAL zxfluxu(klon, klev)
+      REAL zxfluxv(klon, klev)
+C§§§
       REAL heat(klon,klev)    ! chauffage solaire
       REAL heat0(klon,klev)   ! chauffage solaire ciel clair
@@ -424,5 +440,4 @@
 c
       REAL zphi(klon,klev)
-      REAL zx_tmp_x(iim), zx_tmp_yjjmp1
       REAL zx_relief(iim,jjmp1)
       REAL zx_aire(iim,jjmp1)
@@ -561,21 +576,5 @@
 c
        IF (debut) THEN
-c
- 
-         IF (ok_oasis) THEN
-            PRINT*, "Attentions! les parametres suivants sont fixes:"
-            PRINT *,'***********************************************'
-            PRINT*, "npas, nexca, itimestep=", npas, nexca, itimestep
-            PRINT*, "Changer-les manuellement s il le faut"
-            PRINT *,'***********************************************'
-            CALL inicma( npas, nexca, itimestep)
-         ENDIF
-c
-         IF (ok_ocean) THEN
-            PRINT*, '************************'
-            PRINT*, 'SLAB OCEAN est active, prenez precautions !'
-            PRINT*, '************************'
-         ENDIF
-c
+
          DO k = 2, nvm          ! pas de vegetation
             DO i = 1, klon
@@ -588,4 +587,5 @@
          PRINT*, 'Pas de vegetation; desert partout'
 c
+c
 c Initialiser les compteurs:
 c
@@ -595,6 +595,7 @@
 c
          CALL phyetat0 ("startphy.nc",dtime,co2_ppm,solaire,
-     .         rlat,rlon,ftsol,ftsoil,deltat,fqsol,fsnow,
-     .        radsol,rugmer,agesno,clesphy0,
+     .       rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsol,fsnow,
+     .       falbe, fevap, rain_fall,snow_fall,sollw, solsw,
+     .       radsol,rugmer,agesno,clesphy0,
      .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
      .       t_ancien, q_ancien, ancien_ok )
@@ -646,11 +647,4 @@
          ENDIF
 c
-         IF (soil_model) THEN
-            DO nsrf = 1, nbsrf
-            CALL soil(dtime, nsrf, fsnow(1,nsrf),
-     .                ftsol(1,nsrf), ftsoil(1,1,nsrf),
-     .                soilcap(1,nsrf), soilflux(1,nsrf))
-            ENDDO
-         ENDIF
 c
          lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
@@ -716,4 +710,20 @@
      .                "ave(X)", zsto,zout)
 c
+         CALL histdef(nid_day, "tter", "Surface Temperature", "K",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_day, "tlic", "Surface Temperature", "K",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_day, "toce", "Surface Temperature", "K",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "ave(X)", zsto,zout)
+c
+         CALL histdef(nid_day, "tsic", "Surface Temperature", "K",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "ave(X)", zsto,zout)
+c
          CALL histdef(nid_day, "psol", "Surface Pressure", "Pa",
      .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
@@ -768,4 +778,35 @@
      .                "ave(X)", zsto,zout)
 c
+C §§§ PB flux pour chauqe sous surface
+C
+         DO nsrf = 1, nbsrf
+C
+           call histdef(nid_day, "pourc_"//clnsurf(nsrf), 
+     $         "Fraction"//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+
+           call histdef(nid_day, "sens_"//clnsurf(nsrf), 
+     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+c
+           call histdef(nid_day, "lat_"//clnsurf(nsrf), 
+     $         "Latent heat flux "//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+C
+           call histdef(nid_day, "taux_"//clnsurf(nsrf), 
+     $         "Zonal wind stress"//clnsurf(nsrf),"Pa",
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+
+           call histdef(nid_day, "tauy_"//clnsurf(nsrf), 
+     $         "Meridional xind stress "//clnsurf(nsrf), "Pa",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+C§§§
+         END DO 
+           
          CALL histdef(nid_day, "ruis", "Runoff", "mm/day",
      .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
@@ -956,4 +997,32 @@
      .                "ave(X)", zsto,zout)
 c
+         DO nsrf = 1, nbsrf
+C
+           call histdef(nid_mth, "pourc_"//clnsurf(nsrf), 
+     $         "Fraction "//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+C
+           call histdef(nid_mth, "sens_"//clnsurf(nsrf), 
+     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+c
+           call histdef(nid_mth, "lat_"//clnsurf(nsrf), 
+     $         "Latent heat flux "//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+C
+           call histdef(nid_mth, "taux_"//clnsurf(nsrf), 
+     $         "Zonal wind stress"//clnsurf(nsrf), "Pa",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+
+           call histdef(nid_mth, "tauy_"//clnsurf(nsrf), 
+     $         "Meridional xind stress "//clnsurf(nsrf), "Pa",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+         END DO
+C
          CALL histdef(nid_mth, "ruis", "Runoff", "mm/day",
      .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
@@ -1209,5 +1278,9 @@
 c Champs 2D:
 c
-         CALL histdef(nid_ins, "psol", "Surface Pressure", "Pa",
+        CALL histdef(nid_ins, "tsol", "Surface Temperature", "K",
+     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "inst(X)", zsto,zout)
+c
+        CALL histdef(nid_ins, "psol", "Surface Pressure", "Pa",
      .                iim,jjmp1,nhori, 1,1,1, -99, 32,
      .                "inst(X)", zsto,zout)
@@ -1256,4 +1329,38 @@
      .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
      .                "inst(X)", zsto,zout)
+
+         DO nsrf = 1, nbsrf
+C
+           call histdef(nid_ins, "pourc_"//clnsurf(nsrf), 
+     $         "Fraction"//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "inst(X)", zsto,zout)
+
+           call histdef(nid_ins, "sens_"//clnsurf(nsrf), 
+     $         "Sensible heat flux "//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "inst(X)", zsto,zout)
+c
+           call histdef(nid_ins, "tsol_"//clnsurf(nsrf), 
+     $         "Surface Temperature"//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "inst(X)", zsto,zout)
+c
+           call histdef(nid_ins, "lat_"//clnsurf(nsrf), 
+     $         "Latent heat flux "//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "inst(X)", zsto,zout)
+C
+           call histdef(nid_ins, "taux_"//clnsurf(nsrf), 
+     $         "Zonal wind stress"//clnsurf(nsrf),"Pa",
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "inst(X)", zsto,zout)
+
+           call histdef(nid_ins, "tauy_"//clnsurf(nsrf), 
+     $         "Meridional xind stress "//clnsurf(nsrf), "Pa",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "inst(X)", zsto,zout)
+C§§§
+         END DO 
 c
 c Champs 3D:
@@ -1306,19 +1413,4 @@
 cc         ENDDO
 c
-         IF (ok_oasis) THEN
-         DO i = 1, klon
-           oas_sols(i) = 0.0
-           oas_nsol(i) = 0.0
-           oas_rain(i) = 0.0
-           oas_snow(i) = 0.0
-           oas_evap(i) = 0.0
-           oas_ruis(i) = 0.0
-           oas_tsol(i) = 0.0
-           oas_fder(i) = 0.0
-           oas_albe(i) = 0.0
-           oas_taux(i) = 0.0
-           oas_tauy(i) = 0.0
-         ENDDO
-         ENDIF
 c
       ENDIF
@@ -1421,66 +1513,4 @@
          CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
       ENDIF
-cccccccccc
-      IF (ok_oasis .AND. MOD(itap-1,nexca).EQ.0) THEN
-C
-         CALL fromcpl(itap,jjmp1*iim,
-     .        cpl_sst,cpl_sic,cpl_alb_sst,cpl_alb_sic)
-         DO i = 1, iim-1 ! un seul point pour le pole nord
-            cpl_sst(i,1) = cpl_sst(iim,1)
-            cpl_sic(i,1) = cpl_sic(iim,1)
-            cpl_alb_sst(i,1) = cpl_alb_sst(iim,1)
-            cpl_alb_sic(i,1) = cpl_alb_sic(iim,1)
-         ENDDO
-         DO i = 2, iim ! un seul point pour le pole sud
-            cpl_sst(i,jjmp1) = cpl_sst(1,jjmp1)
-            cpl_sic(i,jjmp1) = cpl_sic(1,jjmp1)
-            cpl_alb_sst(i,jjmp1) = cpl_alb_sst(1,jjmp1)
-            cpl_alb_sic(i,jjmp1) = cpl_alb_sic(1,jjmp1)
-         ENDDO
-c
-         ig = 1
-         IF (pctsrf(ig,is_oce).GT.epsfra .OR.
-     .       pctsrf(ig,is_sic).GT.epsfra) THEN
-            pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
-     .                        - (cpl_sic(1,1)-pctsrf(ig,is_sic))
-            pctsrf(ig,is_sic) = cpl_sic(1,1)
-            lmt_sst(ig) = cpl_sst(1,1)
-         ENDIF
-         DO j = 2, jjm
-         DO i = 1, iim
-         ig = ig + 1
-         IF (pctsrf(ig,is_oce).GT.epsfra .OR.
-     .       pctsrf(ig,is_sic).GT.epsfra) THEN
-           pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
-     .                       - (cpl_sic(i,j)-pctsrf(ig,is_sic))
-           pctsrf(ig,is_sic) = cpl_sic(i,j)
-           lmt_sst(ig) = cpl_sst(i,j)
-         ENDIF
-         ENDDO
-         ENDDO
-         ig = ig + 1
-         IF (pctsrf(ig,is_oce).GT.epsfra .OR.
-     .       pctsrf(ig,is_sic).GT.epsfra) THEN
-            pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
-     .                        - (cpl_sic(1,jjmp1)-pctsrf(ig,is_sic))
-            pctsrf(ig,is_sic) = cpl_sic(1,jjmp1)
-            lmt_sst(ig) = cpl_sst(1,jjmp1)
-         ENDIF
-c
-      ENDIF ! ok_oasis
-cccccccccc
-c
- 
-      IF (ok_ocean) THEN
-         DO i = 1, klon
-            ftsol(i,is_oce) = lmt_sst(i) + deltat(i)
-         ENDDO
-
-      ELSE
-         DO i = 1, klon
-            ftsol(i,is_oce) = lmt_sst(i)
-         ENDDO
-
-      ENDIF
 c
 c Re-evaporer l'eau liquide nuageuse
@@ -1523,8 +1553,10 @@
 c
       CALL clmain(dtime,pctsrf,
-     e            t_seri,q_seri,u_seri,v_seri,soil_model,
-     e            ftsol,soilcap,soilflux,paprs,pplay,radsol,
-     e            fsnow,fqsol,
-     e            rlat, frugs,
+     e            t_seri,q_seri,u_seri,v_seri,ok_veget,
+     e            ftsol,paprs,pplay,radsol,
+     e            fsnow,fqsol,fevap,falbe,
+     e            rain_fall, snow_fall, solsw, sollw,
+     e            rlon, rlat, frugs,
+     e            debut, lafin,
      s            d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
      s            fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,rugmer,
@@ -1532,7 +1564,29 @@
      s            ycoefh,yu1,yv1) 
 c
-      DO i = 1, klon
-         sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
-         evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
+C§§§ PB
+C§§§ Incrementation des flux
+C§§
+      zxfluxt=0.
+      zxfluxq=0.
+      zxfluxu=0.
+      zxfluxv=0.
+      DO nsrf = 1, nbsrf
+        DO k = 1, klev
+          DO i = 1, klon
+            zxfluxt(i,k) = zxfluxt(i,k) + 
+     $          fluxt(i,k,nsrf) * pctsrf( i, nsrf)
+            zxfluxq(i,k) = zxfluxq(i,k) + 
+     $          fluxq(i,k,nsrf) * pctsrf( i, nsrf)
+            zxfluxu(i,k) = zxfluxu(i,k) + 
+     $          fluxu(i,k,nsrf) * pctsrf( i, nsrf)
+            zxfluxv(i,k) = zxfluxv(i,k) + 
+     $          fluxv(i,k,nsrf) * pctsrf( i, nsrf)
+          END DO 
+        END DO 
+      END DO 
+      DO i = 1, klon
+         sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
+c         evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
+         evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
          fder(i) = dsens(i) + devap(i)
       ENDDO
@@ -1551,4 +1605,10 @@
       DO i = 1, klon
          zxtsol(i) = 0.0
+         IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) + 
+     $       pctsrf(i, is_oce) + pctsrf(i, is_sic)  - 1.) .GT. EPSFRA) 
+     $       THEN 
+             WRITE(*,*) 'physiq : pb sous surface au point ', i, 
+     $           pctsrf(i, 1 : nbsrf)
+         ENDIF 
       ENDDO
       DO nsrf = 1, nbsrf
@@ -1568,14 +1628,4 @@
       ENDDO
 
-c
-c Appeler le modele du sol
-c
-      IF (soil_model) THEN
-         DO nsrf = 1, nbsrf
-         CALL soil(dtime, nsrf, fsnow(1,nsrf),
-     .             ftsol(1,nsrf), ftsoil(1,1,nsrf),
-     .             soilcap(1,nsrf), soilflux(1,nsrf))
-         ENDDO
-      ENDIF
 c
 c Calculer la derive du flux infrarouge
@@ -1623,5 +1673,5 @@
       ELSE IF (iflag_con.EQ.2) THEN
       CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
-     e            conv_t, conv_q, fluxq(1,1), omega,
+     e            conv_t, conv_q, zxfluxq(1,1), omega,
      s            d_t_con, d_q_con, rain_con, snow_con,
      s            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
@@ -1820,18 +1870,23 @@
       CALL albsno(veget,agesno,alb_neig)
       DO i = 1, klon
-         zx_alb_oce = alb_eau(i)
+         falbe(i,is_oce) = alb_eau(i)
          IF (pctsrf(i,is_oce).GT.epsfra .AND. ftsol(i,is_oce).LT.271.35)
-     .   zx_alb_oce = 0.6 ! pour slab_ocean
+     .   falbe(i,is_oce) = 0.6 ! pour slab_ocean
          zfra = MAX(0.0,MIN(1.0,fsnow(i,is_lic)/(fsnow(i,is_lic)+10.0)))
-         zx_alb_lic = alb_neig(i)*zfra + 0.6*(1.0-zfra)
+         falbe(i,is_lic) = alb_neig(i)*zfra + 0.6*(1.0-zfra)
          zfra = MAX(0.0,MIN(1.0,fsnow(i,is_ter)/(fsnow(i,is_ter)+10.0)))
-         zx_alb_ter = alb_neig(i)*zfra + lmt_alb(i)*(1.0-zfra)
+         falbe(i,is_ter) = alb_neig(i)*zfra + lmt_alb(i)*(1.0-zfra)
          zfra = MAX(0.0,MIN(1.0,fsnow(i,is_sic)/(fsnow(i,is_sic)+10.0)))
-         zx_alb_sic = alb_neig(i)*zfra + 0.6*(1.0-zfra)
-         albsol(i) = zx_alb_oce * pctsrf(i,is_oce)
-     .             + zx_alb_lic * pctsrf(i,is_lic)
-     .             + zx_alb_ter * pctsrf(i,is_ter)
-     .             + zx_alb_sic * pctsrf(i,is_sic)
-      ENDDO
+         falbe(i,is_sic) = alb_neig(i)*zfra + 0.6*(1.0-zfra)
+         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)
+      ENDDO
+c      DO nsrf = 1, nbsrf
+c        DO i = 1, klon
+c           albsol(i) = albsol(i) + falbe(i,nsrf)*pctsrf(i,nsrf)
+c        ENDDO
+c      ENDDO
       CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
      e            (dist, rmu0, fract, co2_ppm, solaire,
@@ -1856,6 +1911,6 @@
 c Calculer l'hydrologie de la surface
 c
-      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, evap,
-     .            agesno, ftsol,fqsol,fsnow, ruis)
+c      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
+c     .            agesno, ftsol,fqsol,fsnow, ruis)
 c
       DO i = 1, klon
@@ -2007,19 +2062,4 @@
 c Accumuler les variables a stocker dans les fichiers histoire:
 c
-      IF (ok_oasis) THEN ! couplage oasis
-      DO i = 1, klon
-        oas_sols(i) = oas_sols(i) + solsw(i)          / FLOAT(nexca)
-        oas_nsol(i) = oas_nsol(i) + (bils(i)-solsw(i))/ FLOAT(nexca)
-        oas_rain(i) = oas_rain(i) + rain_fall(i)      / FLOAT(nexca)
-        oas_snow(i) = oas_snow(i) + snow_fall(i)      / FLOAT(nexca)
-        oas_evap(i) = oas_evap(i) + evap(i)           / FLOAT(nexca)
-        oas_tsol(i) = oas_tsol(i) + zxtsol(i)         / FLOAT(nexca)
-        oas_fder(i) = oas_fder(i) + fder(i)           / FLOAT(nexca)
-        oas_albe(i) = oas_albe(i) + albsol(i)         / FLOAT(nexca)
-        oas_taux(i) = oas_taux(i) + fluxu(i,1)        / FLOAT(nexca)
-        oas_tauy(i) = oas_tauy(i) + fluxv(i,1)        / FLOAT(nexca)
-        oas_ruis(i) = oas_ruis(i) + ruis(i)       /FLOAT(nexca)/dtime
-      ENDDO
-      ENDIF
 c
 c
@@ -2043,4 +2083,21 @@
       CALL histwrite(nid_day,"tsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
+C
+      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_ter)
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d ,zx_tmp_2d)
+      CALL histwrite(nid_day,"tter",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+C
+      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_lic)
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
+      CALL histwrite(nid_day,"tlic",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+C
+      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_oce)
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
+      CALL histwrite(nid_day,"toce",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+C
+      zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, is_sic)
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
+      CALL histwrite(nid_day,"tsic",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+C
       DO i = 1, klon
          zx_tmp_fi2d(i) = paprs(i,1)
@@ -2085,21 +2142,50 @@
       CALL histwrite(nid_day,"ruis",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
-      DO i = 1, klon
-         zx_tmp_fi2d(i) = fluxu(i,1)
-      ENDDO
-      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
-      CALL histwrite(nid_day,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
-c
-      DO i = 1, klon
-         zx_tmp_fi2d(i) = fluxv(i,1)
-      ENDDO
-      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
-      CALL histwrite(nid_day,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
-c
-      DO i = 1, klon
-         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
-      ENDDO
-      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
-      CALL histwrite(nid_day,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+c      DO i = 1, klon
+c         zx_tmp_fi2d(i) = fluxu(i,1)
+c      ENDDO
+c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
+c      CALL histwrite(nid_day,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+c      DO i = 1, klon
+c         zx_tmp_fi2d(i) = fluxv(i,1)
+c      ENDDO
+c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
+c      CALL histwrite(nid_day,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      DO nsrf = 1, nbsrf
+C§§§
+        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_day,"pourc_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C 
+        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_day,"sens_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C
+        zx_tmp_fi2d(1 : klon) = - fluxq( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_day,"lat_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C
+        zx_tmp_fi2d(1 : klon) = - fluxu( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_day,"taux_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C      
+        zx_tmp_fi2d(1 : klon) = - fluxv( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_day,"tauy_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d)
+C 
+      END DO  
+C
+c$$$      DO i = 1, klon
+c$$$         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
+c$$$      ENDDO
+c$$$      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
+c$$$      CALL histwrite(nid_day,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldl,zx_tmp_2d)
@@ -2243,21 +2329,49 @@
       CALL histwrite(nid_mth,"ruis",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
-      DO i = 1, klon
-         zx_tmp_fi2d(i) = fluxu(i,1)
-      ENDDO
-      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
-      CALL histwrite(nid_mth,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
-c
-      DO i = 1, klon
-         zx_tmp_fi2d(i) = fluxv(i,1)
-      ENDDO
-      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
-      CALL histwrite(nid_mth,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
-c
-      DO i = 1, klon
-         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
-      ENDDO
-      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
-      CALL histwrite(nid_mth,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+c      DO i = 1, klon
+c         zx_tmp_fi2d(i) = fluxu(i,1)
+c      ENDDO
+c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
+c      CALL histwrite(nid_mth,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+c      DO i = 1, klon
+c         zx_tmp_fi2d(i) = fluxv(i,1)
+c      ENDDO
+c      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
+c      CALL histwrite(nid_mth,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      DO nsrf = 1, nbsrf
+C§§§
+        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_mth,"pourc_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C 
+        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_mth,"sens_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C
+        zx_tmp_fi2d(1 : klon) = - fluxq( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_mth,"lat_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C
+        zx_tmp_fi2d(1 : klon) = - fluxu( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_mth,"taux_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C      
+        zx_tmp_fi2d(1 : klon) = - fluxv( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_mth,"tauy_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d)
+C 
+      END DO  
+c$$$      DO i = 1, klon
+c$$$         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
+c$$$      ENDDO
+c$$$      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
+c$$$      CALL histwrite(nid_mth,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol,zx_tmp_2d)
@@ -2476,4 +2590,7 @@
       CALL histwrite(nid_ins,"psol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
+      CALL histwrite(nid_ins,"tsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
       CALL histwrite(nid_ins,"topl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
@@ -2508,4 +2625,38 @@
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, d_ts(1,is_sic),zx_tmp_2d)
       CALL histwrite(nid_ins,"dtsvdfi",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+
+      DO nsrf = 1, nbsrf
+C§§§
+        zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_ins,"pourc_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C 
+        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_ins,"sens_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C
+        zx_tmp_fi2d(1 : klon) = - fluxq( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_ins,"lat_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C
+        zx_tmp_fi2d(1 : klon) = ftsol( 1 : klon, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_ins,"tsol_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C
+        zx_tmp_fi2d(1 : klon) = - fluxu( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_ins,"taux_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C      
+        zx_tmp_fi2d(1 : klon) = - fluxv( 1 : klon, 1, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_ins,"tauy_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d)
+C 
+      END DO  
 
 c
@@ -2546,82 +2697,4 @@
       ENDIF
 c
-      IF (ok_oasis .AND. mod(itap,nexca).EQ.0) THEN
-c
-c Je ne traite pas le ruissellement, pour l'instant (qui m'aidera ?)
-         DO i = 1, klon
-            oas_ruisoce(i) = 0.0
-            oas_ruisriv(i) = 0.0
-         ENDDO
-c
-         ig = 1
-         DO i = 1, iim
-            z_sols(i,1) = oas_sols(ig)
-            z_nsol(i,1) = oas_nsol(ig)
-            z_rain(i,1) = oas_rain(ig)
-            z_snow(i,1) = oas_snow(ig)
-            z_evap(i,1) = oas_evap(ig)
-            z_ruisoce(i,1) = oas_ruisoce(ig)
-            z_ruisriv(i,1) = oas_ruisriv(ig)
-            z_tsol(i,1) = oas_tsol(ig)
-            z_fder(i,1) = oas_fder(ig)
-            z_albe(i,1) = oas_albe(ig)
-            z_taux(i,1) = oas_taux(ig)
-            z_tauy(i,1) = oas_tauy(ig)
-         ENDDO
-         DO j = 2, jjm
-         DO i = 1, iim
-            ig = ig + 1
-            z_sols(i,j) = oas_sols(ig)
-            z_nsol(i,j) = oas_nsol(ig)
-            z_rain(i,j) = oas_rain(ig)
-            z_snow(i,j) = oas_snow(ig)
-            z_evap(i,j) = oas_evap(ig)
-            z_ruisoce(i,j) = oas_ruisoce(ig)
-            z_ruisriv(i,j) = oas_ruisriv(ig)
-            z_tsol(i,j) = oas_tsol(ig)
-            z_fder(i,j) = oas_fder(ig)
-            z_albe(i,j) = oas_albe(ig)
-            z_taux(i,j) = oas_taux(ig)
-            z_tauy(i,j) = oas_tauy(ig)
-         ENDDO
-         ENDDO
-         ig = ig + 1
-         DO i = 1, iim
-            z_sols(i,jjmp1)    = oas_sols(ig)
-            z_nsol(i,jjmp1)    = oas_nsol(ig)
-            z_rain(i,jjmp1)    = oas_rain(ig)
-            z_snow(i,jjmp1)    = oas_snow(ig)
-            z_evap(i,jjmp1)    = oas_evap(ig)
-            z_ruisoce(i,jjmp1) = oas_ruisoce(ig)
-            z_ruisriv(i,jjmp1) = oas_ruisriv(ig)
-            z_tsol(i,jjmp1)    = oas_tsol(ig)
-            z_fder(i,jjmp1)    = oas_fder(ig)
-            z_albe(i,jjmp1)    = oas_albe(ig)
-            z_taux(i,jjmp1)    = oas_taux(ig)
-            z_tauy(i,jjmp1)    = oas_tauy(ig)
-         ENDDO
-c
-c Passer les champs au coupleur:
-c
-         CALL intocpl(itap,jjmp1*iim,
-     .                   z_sols, z_nsol,
-     .                   z_rain, z_snow, z_evap,
-     .                   z_ruisoce, z_ruisriv,
-     .                   z_tsol, z_fder, z_albe,
-     .                   z_taux, z_tauy)
-         DO i = 1, klon
-           oas_sols(i) = 0.0
-           oas_nsol(i) = 0.0
-           oas_rain(i) = 0.0
-           oas_snow(i) = 0.0
-           oas_evap(i) = 0.0
-           oas_ruis(i) = 0.0
-           oas_tsol(i) = 0.0
-           oas_fder(i) = 0.0
-           oas_albe(i) = 0.0
-           oas_taux(i) = 0.0
-           oas_tauy(i) = 0.0
-         ENDDO
-      ENDIF
 c
 c Ecrire la bande regionale (binaire grads)
@@ -2639,6 +2712,6 @@
          CALL ecriregs(84,bils)
          CALL ecriregs(84,pctsrf(1,is_sic))
-         CALL ecriregs(84,fluxu(1,1))
-         CALL ecriregs(84,fluxv(1,1))
+         CALL ecriregs(84,zxfluxu(1,1))
+         CALL ecriregs(84,zxfluxv(1,1))
          CALL ecriregs(84,ue)
          CALL ecriregs(84,ve)
@@ -2705,8 +2778,10 @@
 ccc         IF (ok_oasis) CALL quitcpl
          CALL phyredem ("restartphy.nc",dtime,radpas,co2_ppm,solaire,
-     .        rlat,rlon,ftsol,ftsoil,deltat,fqsol,fsnow,
-     .        radsol,rugmer,agesno,
-     .        zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
-     .        t_ancien, q_ancien)
+     .      rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsol, fsnow,
+     .      falbe, fevap, rain_fall, snow_fall,
+     .      solsw, sollw,
+     .      radsol,rugmer,agesno,
+     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
+     .      t_ancien, q_ancien)
       ENDIF
 
