Index: LMDZ.3.3/branches/rel-LF/libf/bibio/initdynav.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/bibio/initdynav.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/bibio/initdynav.F	(revision 177)
@@ -3,4 +3,5 @@
 
        USE IOIPSL
+       USE histcom
 
       implicit none
Index: LMDZ.3.3/branches/rel-LF/libf/bibio/initfluxsto.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/bibio/initfluxsto.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/bibio/initfluxsto.F	(revision 177)
@@ -4,4 +4,5 @@
 
        USE IOIPSL
+       USE histcom
 
       implicit none
Index: LMDZ.3.3/branches/rel-LF/libf/bibio/inithist.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/bibio/inithist.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/bibio/inithist.F	(revision 177)
@@ -3,4 +3,5 @@
 
        USE IOIPSL
+       USE histcom
 
       implicit none
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/calfis.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/calfis.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/calfis.F	(revision 177)
@@ -128,4 +128,5 @@
       REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqmx)
       REAL zdpsrf(ngridmx)
+      REAL zcufi(ngridmx),zcvfi(ngridmx)
 c
       REAL zsin(iim),zcos(iim),z1(iim)
@@ -177,12 +178,18 @@
          latfi(1)=rlatu(1)
          lonfi(1)=0.
+         zcufi(1) = cu(1,1)
+         zcvfi(1) = cv(1,1)
          DO j=2,jjm
             DO i=1,iim
                latfi((j-2)*iim+1+i)= rlatu(j)
                lonfi((j-2)*iim+1+i)= rlonv(i)
+               zcufi((j-2)*iim+1+i) = cu(i,j)
+               zcvfi((j-2)*iim+1+i) = cv(i,j)
      	    ENDDO
          ENDDO
          latfi(ngridmx)= rlatu(jjp1)
          lonfi(ngridmx)= 0.
+         zcufi(ngridmx) = cu(1,jjp1)
+         zcvfi(ngridmx) = cv(1,jjm)
          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
           PRINT*,'WARNING!!! vitesse verticale nulle dans la physique'
@@ -430,5 +437,5 @@
      ,             presnivs,clesphy0,  zufi, zvfi,ztfi, zqfi,  
 ccc     ,             pcvgu, pcvgv, pcvgt, pcvgq, 
-     ,             pvervel,
+     ,             pvervel, zcufi, zcvfi,
 C - sorties
      s             zdufi, zdvfi, zdtfi, zdqfi,zdpsrf              )
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/create_limit.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/create_limit.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/create_limit.F	(revision 177)
@@ -101,5 +101,4 @@
       REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
       REAL :: flic_tmp(iip1, jjp1)
-      REAL :: champint(iim, jjp1)
 c Diverses variables locales
       REAL time
@@ -153,6 +152,9 @@
       zmasq(:) = 0.
       CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmasq,0.0)
-      WHERE (zmasq(1 : klon) .LE. EPSFRA)
+      WHERE (zmasq(1 : klon) .LT. EPSFRA)
           zmasq(1 : klon) = 0.
+      END WHERE 
+      WHERE (1 - zmasq(1 : klon) .LT. EPSFRA)
+          zmasq(1 : klon) = 1.
       END WHERE 
 !      WRITE(*,*)zmasq
@@ -169,10 +171,17 @@
           END DO
       ENDIF 
-      DO i = 1, iim
+c$$$      DO i = 1, iim
+c$$$      DO j = 1, jjp1
+c$$$         mask(i,j) = masque(i,j)
+c$$$      ENDDO
+c$$$      ENDDO
+c$$$      CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, phy_nat0)
+      phy_nat0(1:klon) = zmasq(1:klon)
+      mask = 0.
       DO j = 1, jjp1
-         mask(i,j) = masque(i,j)
-      ENDDO
-      ENDDO
-      CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, phy_nat0)
+        DO i = 1, iim
+          IF ( masque(i,j) .GE. EPSFRA) mask (i,j) = 1
+        END DO 
+      END DO  
 C
 C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
@@ -275,8 +284,8 @@
      $    pctsrf(1:klon, is_lic))
 C adequation avec le maque terre/mer
-      WHERE (pctsrf(1 : klon, is_lic) .LE. EPSFRA ) 
+      WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA ) 
           pctsrf(1 : klon, is_lic) = 0. 
       END WHERE
-      WHERE (zmasq( 1 : klon) .LE. EPSFRA) 
+      WHERE (zmasq( 1 : klon) .LT. EPSFRA) 
           pctsrf(1 : klon, is_lic) = 0.
       END WHERE 
@@ -289,4 +298,8 @@
             ELSE 
                 pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
+                IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
+                    pctsrf(ji,is_ter) = 0.
+                    pctsrf(ji, is_lic) = zmasq(ji)
+                ENDIF 
             ENDIF 
         ENDIF 
@@ -557,5 +570,5 @@
 CPB  en attendant de mettre fraction de terre
 c
-          WHERE(phy_ice(1:klon) .GT. 1.) phy_ice(1 : klon) = 1.
+          WHERE(phy_ice(1:klon) .GE. 1.) phy_ice(1 : klon) = 1.
           WHERE(phy_ice(1:klon) .LT. EPSFRA) phy_ice(1 : klon) = 0.
 c 
@@ -564,15 +577,44 @@
             pctsrf_t(:,is_ter,k) = pctsrf(:,is_ter)
             pctsrf_t(:,is_lic,k) = pctsrf(:,is_lic)
+            pctsrf_t(1:klon,is_sic,k) =   phy_ice(1:klon) 
+     $            - pctsrf_t(1:klon,is_lic,k)
+c§§ Il y a des cas ou il y a de la glace dans landiceref et pas dans AMIP
+            WHERE (pctsrf_t(1:klon,is_sic,k) .LE. 0)
+                pctsrf_t(1:klon,is_sic,k) = 0.
+            END WHERE 
+            WHERE( 1. - zmasq(1:klon) .LT. EPSFRA)
+                pctsrf_t(1:klon,is_sic,k) = 0.
+                pctsrf_t(1:klon,is_oce,k) = 0.
+            END WHERE
             DO i = 1, klon
-              pctsrf_t(i,is_sic,k) = (1. - pctsrf_t(i,is_lic,k) -
-     .                               pctsrf_t(i,is_ter,k)) * phy_ice(i)
-              pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_lic,k) -
-     .                      pctsrf_t(i,is_ter,k) - pctsrf_t(i,is_sic,k)  
+c$$              pctsrf_t(i,is_sic,k) = (1. - pctsrf_t(i,is_lic,k) -
+c$$     .                               pctsrf_t(i,is_ter,k)) * phy_ice(i)
+c$$              pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_lic,k) -
+c$$     .                      pctsrf_t(i,is_ter,k) - pctsrf_t(i,is_sic,k)
+              IF ( 1. - zmasq(i) .GT. EPSFRA) THEN 
+                  IF ( pctsrf_t(i,is_sic,k) .GE. 1 - zmasq(i)) THEN
+                      pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
+                      pctsrf_t(i,is_oce,k) = 0.
+                  ELSE 
+                      pctsrf_t(i,is_oce,k) = 1 - zmasq(i) 
+     $                    - pctsrf_t(i,is_sic,k)
+                      IF (pctsrf_t(i,is_oce,k) .LT. EPSFRA) THEN
+                          pctsrf_t(i,is_oce,k) = 0.
+                          pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
+                      ENDIF 
+                  ENDIF
+              ENDIF  
               if (pctsrf_t(i,is_oce,k) .lt. 0.) then
-                      WRITE(*,*) 'pb sous maille au point : i,k '
-     $                    , i,k,pctsrf_t(:,is_oce,k)
+                  WRITE(*,*) 'pb sous maille au point : i,k '
+     $                , i,k,pctsrf_t(:,is_oce,k)
+              ENDIF
+              IF ( abs( pctsrf_t(i, is_ter,k) + pctsrf_t(i, is_lic,k) + 
+     $            pctsrf_t(i, is_oce,k) + pctsrf_t(i, is_sic,k)  - 1.) 
+     $            .GT. EPSFRA) THEN 
+                  WRITE(*,*) 'physiq : pb sous surface au point ', i, 
+     $                pctsrf_t(i, 1 : nbsrf,k), phy_ice(i)
               ENDIF 
             END DO
-          ELSE 
+        ELSE 
             DO i = 1, klon
               pctsrf_t(i,is_ter,k) = pctsrf(i,is_ter)
@@ -750,4 +792,5 @@
       ENDDO
 c
+      WHERE(phy_sst .LT. 271.35) phy_sst = 271.35
       ierr = NF_CLOSE(ncid)
 c
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/defrun_new.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/defrun_new.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/defrun_new.F	(revision 177)
@@ -75,5 +75,5 @@
 c-----------------------------------------------------------------------
 
-      OPEN( tapedef,file ='run.def',status='old',form='formatted')
+      OPEN( tapedef,file ='gcm.def',status='old',form='formatted')
 
 
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/etat0_netcdf.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/etat0_netcdf.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/etat0_netcdf.F	(revision 177)
@@ -186,6 +186,9 @@
       zmasq(:) = 0.
       CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmasq,0.0)
-      WHERE (zmasq(1 : klon) .LE. EPSFRA)
+      WHERE (zmasq(1 : klon) .LT. EPSFRA)
           zmasq(1 : klon) = 0.
+      END WHERE 
+      WHERE (1. - zmasq(1 : klon) .LT. EPSFRA)
+          zmasq(1 : klon) = 1.
       END WHERE 
       WRITE(*,*)zmasq
@@ -434,8 +437,9 @@
      $    pctsrf(1:klon, is_lic))
 C adequation avec le maque terre/mer
-      WHERE (pctsrf(1 : klon, is_lic) .LE. EPSFRA ) 
+c      zmasq(157) = 0.
+      WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA ) 
           pctsrf(1 : klon, is_lic) = 0. 
       END WHERE
-      WHERE (zmasq( 1 : klon) .LE. EPSFRA) 
+      WHERE (zmasq( 1 : klon) .LT. EPSFRA) 
           pctsrf(1 : klon, is_lic) = 0.
       END WHERE 
@@ -448,4 +452,8 @@
             ELSE 
                 pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
+                IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
+                    pctsrf(ji,is_ter) = 0.
+                    pctsrf(ji, is_lic) = zmasq(ji)
+                ENDIF 
             ENDIF 
         ENDIF 
@@ -558,6 +566,6 @@
       albe(:,is_sic) = 0.6
       evap(:,:) = 0.
-      qsolsrf(:,is_ter) = qsol
-      qsolsrf(:,is_lic) = qsol
+      qsolsrf(:,is_ter) = 150
+      qsolsrf(:,is_lic) = 150
       qsolsrf(:,is_oce) = 150.
       qsolsrf(:,is_sic) = 150.
@@ -577,5 +585,5 @@
       frugs(:,is_ter) = rugmer
       frugs(:,is_lic) = rugmer
-      frugs(:,is_sic) = rugmer
+      frugs(:,is_sic) = 0.001
 
       call physdem("startphy.nc",phystep,radpas, co2_ppm, solaire,
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F	(revision 177)
@@ -358,4 +358,6 @@
      $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
      $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
+      write(47,*)'dp apres caldyn'
+      write(47,*)dp
 
 c-----------------------------------------------------------------------
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/grid_noro.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/grid_noro.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/grid_noro.F	(revision 177)
@@ -332,5 +332,5 @@
 c angle theta:
            zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.
-c$$$           zphi(ii,jj)=zmea(ii,jj)
+           zphi(ii,jj)=zmea(ii,jj)
 c$$$           zmea(ii,jj)=zmea(ii,jj)
 c$$$           zpic(ii,jj)=zpic(ii,jj)
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/startvar.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/startvar.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/startvar.F	(revision 177)
@@ -1,2 +1,3 @@
+! $Header$
       MODULE startvar
     !
@@ -319,4 +320,6 @@
      $    phis, relief, zstd, zsig, zgam, zthe, zpic, zval, masque)
       phis = phis * 9.81
+!      write(*,*)'phis sortie grid_noro'
+!      write(*,*)phis
     !
     !PB supression ligne suivant pour masque avec % terre
Index: LMDZ.3.3/branches/rel-LF/libf/filtrez/parafilt.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/filtrez/parafilt.h	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/filtrez/parafilt.h	(revision 177)
@@ -4,5 +4,5 @@
 c       PARAMETER (nfilun=30,nfilus=30,nfilvn=30,nfilvs=30)
 c        PARAMETER (nfilun=6, nfilus=5, nfilvn=5, nfilvs=5)
-         PARAMETER (nfilun=15, nfilus=8, nfilvn=14, nfilvs=8)
+c        PARAMETER (nfilun=15, nfilus=8, nfilvn=14, nfilvs=8)
 c        PARAMETER (nfilun=24, nfilus=23, nfilvn=24, nfilvs=24)
 cmaf -debug  PARAMETER (nfilun=2, nfilus=1, nfilvn=2, nfilvs=2)
@@ -17,5 +17,5 @@
 c
 c 96 72 19 non-zoom:
-ccc      PARAMETER (nfilun=12, nfilus=11, nfilvn=12, nfilvs=12)
+      PARAMETER (nfilun=12, nfilus=11, nfilvn=12, nfilvs=12)
 c
 c        PARAMETER ( nfilun=20, nfilus=20, nfilvn=20, nfilvs=20 )
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 177)
@@ -3,6 +3,7 @@
      .                  jour, rmu0,
      .                  ok_veget, ocean, npas, nexca, ts,
-     .                  paprs,pplay,radsol,snow,qsol,evap,albe,
-     .                  rain_f, snow_f, solsw, sollw, fder,
+     .                  soil_model,ftsoil,
+     .                  paprs,pplay,radsol,snow,qsol,evap,albe,fluxlat,
+     .                  rain_f, snow_f, solsw, sollw, sollwdown, fder,
      .                  rlon, rlat, cufi, cvfi, rugos,
      .                  debut, lafin, agesno,rugoro,
@@ -45,4 +46,6 @@
 c rlat-----input-R- latitude en degree
 c rugos----input-R- longeur de rugosite (en m)
+c cufi-----input-R- resolution des mailles en x (m)
+c cvfi-----input-R- resolution des mailles en y (m)
 c
 c d_t------output-R- le changement pour "t"
@@ -56,5 +59,4 @@
 c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
 c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
-c rugmer---output-R- longeur de rugosite sur mer (m)
 c dflux_t derive du flux sensible
 c dflux_q derive du flux latent
@@ -70,5 +72,6 @@
 #include "dimphy.h"
 #include "indicesol.h"
-c
+c$$$ PB ajout pour soil
+#include "dimsoil.h"
 c
       REAL dtime
@@ -102,7 +105,10 @@
       REAL evap(klon,nbsrf)
       REAL albe(klon,nbsrf)
+c$$$ PB
+      REAL fluxlat(klon,nbsrf)
+C
       real rain_f(klon), snow_f(klon)
       REAL fder(klon)
-      REAL sollw(klon), solsw(klon)
+      REAL sollw(klon), solsw(klon), sollwdown(klon)
       REAL rugos(klon,nbsrf)
 C la nouvelle repartition des surfaces sortie de l'interface
@@ -113,4 +119,8 @@
       REAL zv1(klon)
 cAA
+c$$$ PB ajout pour soil
+      LOGICAL soil_model
+      REAL ftsoil(klon,nsoilmx,nbsrf)
+      REAL ytsoil(klon,nsoilmx)
 c======================================================================
       EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
@@ -121,7 +131,10 @@
       real ysnow(klon), yqsol(klon)
       real yrain_f(klon), ysnow_f(klon)
-      real ysollw(klon), ysolsw(klon)
+      real ysollw(klon), ysolsw(klon), ysollwdown(klon)
       real yfder(klon), ytaux(klon), ytauy(klon)
       REAL yrugm(klon), yrads(klon),yrugoro(klon)
+c$$$ PB
+      REAL yfluxlat(klon)
+C
       REAL y_d_ts(klon)
       REAL y_d_t(klon, klev), y_d_q(klon, klev)
@@ -197,4 +210,5 @@
       ysolsw = 0.0
       ysollw = 0.0
+      ysollwdown = 0.0
       yrugos = 0.0
       yu1 = 0.0
@@ -212,4 +226,5 @@
       y_flux_u = 0.0
       y_flux_v = 0.0
+      ytsoil = 0.0
 
       DO nsrf = 1, nbsrf
@@ -219,4 +234,5 @@
       END DO
 C§§§ PB
+      yfluxlat=0.
       flux_t = 0.
       flux_q = 0.
@@ -254,8 +270,9 @@
 
       pctsrf_pot = pctsrf
-      pctsrf_pot(:,is_sic) = pctsrf(:,is_oce)
+      pctsrf_pot(:,is_oce) = 1. - zmasq(:)
+      pctsrf_pot(:,is_sic) = 1. - zmasq(:)
 
       DO 99999 nsrf = 1, nbsrf
-      totalflu = radsol
+c$$$   PB   totalflu = radsol
 
 c chercher les indices:
@@ -274,4 +291,5 @@
       ENDDO
 c
+      write(*,*)'CLMAIN, nsrf, knon =',nsrf, knon
       IF (knon.EQ.0) GOTO 99999
       DO j = 1, knon
@@ -290,11 +308,25 @@
         ysolsw(j) = solsw(i)
         ysollw(j) = sollw(i)
+        ysollwdown(j) = sollwdown(i)
         yrugos(j) = rugos(i,nsrf)
         yrugoro(j) = rugoro(i)
         yu1(j) = u1lay(i)
         yv1(j) = v1lay(i)
-        yrads(j) = totalflu(i)
+c$$$    PB    yrads(j) = totalflu(i)
+        yrads(j) = (1 - albe(i,nsrf)) 
+     $      /(1 - pctsrf(i,is_ter) * albe(i,is_ter) 
+     $      - pctsrf(i, is_lic) *albe(i,is_lic) 
+     $      - pctsrf(i, is_oce) *albe(i,is_oce)
+     $      - pctsrf(i, is_sic) *albe(i,is_sic) 
+     $      ) * solsw(i) + sollw(i)
         ypaprs(j,klev+1) = paprs(i,klev+1)
-      ENDDO
+      END DO
+c$$$ PB ajour pour soil
+      DO k = 1, nsoilmx
+        DO j = 1, knon
+          i = ni(j)
+          ytsoil(j,k) = ftsoil(i,k,nsrf)
+        END DO  
+      END DO 
       DO k = 1, klev
       DO j = 1, knon
@@ -339,4 +371,5 @@
      e          rlon, rlat, cufi, cvfi,
      e          knon, nsrf, ni, pctsrf,
+     e          soil_model, ytsoil,
      e          ok_veget, ocean, npas, nexca,
      e          rmu0, yrugos, yrugoro,
@@ -345,5 +378,6 @@
      e          ydelp,yrads, yevap,yalb, ysnow, yqsol, 
      e          yrain_f, ysnow_f, yfder, ytaux, ytauy, 
-     e          ysollw, ysolsw,
+c$$$     e          ysollw, ysolsw,
+     e          ysollw, ysollwdown, ysolsw,yfluxlat,
      s          pctsrf_new, agesno,
      s          y_d_t, y_d_q, y_d_ts, yz0_new,
@@ -384,4 +418,5 @@
         ENDDO
       ENDDO
+      
 
       evap(:,nsrf) = - flux_q(:,1,nsrf)
@@ -394,5 +429,7 @@
          qsol(i,nsrf) = yqsol(j)
          rugos(i,nsrf) = yz0_new(j)
-         rugmer(i) = yrugm(j)
+         fluxlat(i,nsrf) = yfluxlat(j)
+c$$$ pb         rugmer(i) = yrugm(j)
+         IF (nsrf .EQ. is_oce) rugmer(i) = yrugm(j)
          cdragh(i) = cdragh(i) + ycoefh(j,1)
          cdragm(i) = cdragm(i) + ycoefm(j,1)
@@ -401,5 +438,12 @@
          zu1(i) = zu1(i) + yu1(j)
          zv1(i) = zv1(i) + yv1(j)
-      ENDDO
+      END DO
+c$$$ PB ajout pour soil
+      DO k = 1, nsoilmx
+        DO j = 1, knon
+          i = ni(j)
+          ftsoil(i, k, nsrf) = ytsoil(j,k)
+        END DO 
+      END DO 
 c
 #ifdef CRAY
@@ -438,4 +482,5 @@
      e                rlon, rlat, cufi, cvfi, 
      e                knon, nisurf, knindex, pctsrf,
+     $                soil_model,tsoil,
      e                ok_veget, ocean, npas, nexca,
      e                rmu0, rugos, rugoro,
@@ -444,5 +489,6 @@
      e                delp,radsol,evap,albedo,snow,qsol, 
      e                precip_rain, precip_snow, fder, taux, tauy,
-     e                lwdown, swdown, 
+c$$$     e                lwdown, swdown, 
+     $                sollw, sollwdown, swdown,fluxlat, 
      s                pctsrf_new, agesno,
      s                d_t, d_q, d_ts, z0_new, 
@@ -462,4 +508,5 @@
 #include "FCTTRE.h"
 #include "indicesol.h"
+#include "dimsoil.h"
 c Arguments:
       INTEGER knon
@@ -542,9 +589,12 @@
       real fder(klon), taux(klon), tauy(klon)
       real temp_air(klon), spechum(klon)
-      real hum_air(klon), ccanopy(klon)
+      real epot_air(klon), ccanopy(klon)
       real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
       real petBcoef(klon), peqBcoef(klon)
-      real lwdown(klon), swnet(klon), swdown(klon)
+      real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)
       real p1lay(klon)
+c$$$C PB ajout pour soil
+      LOGICAL soil_model
+      REAL tsoil(klon, nsoilmx)
 
 ! Parametres de sortie
@@ -553,4 +603,5 @@
       real emis_new(klon), z0_new(klon)
       real pctsrf_new(klon,nbsrf)
+      
 c
 
@@ -669,4 +720,5 @@
         tq_cdrag=coef(:,1)
         temp_air=t(:,1)
+        epot_air=local_h(:,1)
         spechum=q(:,1)
         p1lay = pplay(:,1)
@@ -675,13 +727,13 @@
 c      enddo
 c En attendant mieux
-      hum_air = 0.
-      ccanopy = 0.
+      ccanopy = 365.
 
       CALL interfsurf(itime, dtime, jour, rmu0,
-     e klon, iim, jjm, nisurf, knon, knindex, pctsrf, rlon, rlat,
-     e debut, lafin, ok_veget,
-     e zlev1,  u1lay, v1lay, temp_air, spechum, hum_air, ccanopy, 
+     e klon, iim, jjm, nisurf, knon, knindex, pctsrf, 
+     e rlon, rlat, cufi, cvfi, 
+     e debut, lafin, ok_veget, soil_model, nsoilmx,tsoil,
+     e zlev1,  u1lay, v1lay, temp_air, spechum, epot_air, ccanopy, 
      e tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,
-     e precip_rain, precip_snow, lwdown, swnet, swdown,
+     e precip_rain, precip_snow, sollw, sollwdown, swnet, swdown,
      e fder, taux, tauy, rugos, rugoro,
      e albedo, snow, qsol,
@@ -691,4 +743,5 @@
      s tsol_rad, tsurf_new, alb_new, emis_new, z0_new, 
      s pctsrf_new, agesno)
+
 
       do i = 1, knon
@@ -1044,4 +1097,5 @@
          ENDIF
       ENDDO
+
 c
 c Calculer les coefficients turbulents dans l'atmosphere
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/hgardfou.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/hgardfou.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/hgardfou.F	(revision 177)
@@ -51,5 +51,6 @@
          jbad = 0
          DO i = 1, klon
-         IF (zt(i).LT.100.0) THEN
+!         IF (zt(i).LT.100.0) THEN
+         IF (zt(i).LT.50.0) THEN
             jbad = jbad + 1
             jadrs(jbad) = i
@@ -91,5 +92,6 @@
          jbad = 0
          DO i = 1, klon
-         IF (zt(i).LT.100.0) THEN
+!         IF (zt(i).LT.100.0) THEN
+         IF (zt(i).LT.50.0) THEN
             jbad = jbad + 1
             jadrs(jbad) = i
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 177)
@@ -41,5 +41,5 @@
 ! run_off      ruissellement total
   real, allocatable, dimension(:),save    :: run_off
-
+  real, allocatable, dimension(:),save    :: coastalflow, riverflow
 
 
@@ -49,9 +49,10 @@
 !
   SUBROUTINE interfsurf_hq(itime, dtime, jour, rmu0, &
-      & klon, iim, jjm, nisurf, knon, knindex, pctsrf, rlon, rlat, &
-      & debut, lafin, ok_veget, &
-      & zlev,  u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, & 
+      & klon, iim, jjm, nisurf, knon, knindex, pctsrf, &
+      & rlon, rlat, cufi, cvfi,&
+      & debut, lafin, ok_veget, soil_model, nsoilmx, tsoil,&
+      & zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, & 
       & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
-      & precip_rain, precip_snow, lwdown, swnet, swdown, &
+      & precip_rain, precip_snow, sollw, sollwdown, swnet, swdown, &
       & fder, taux, tauy, rugos, rugoro, &
       & albedo, snow, qsol, &
@@ -84,4 +85,5 @@
 !   rlon         longitudes
 !   rlat         latitudes
+!   cufi,cvfi    resolution des mailles en x et y (m)
 !   debut        logical: 1er appel a la physique
 !   lafin        logical: dernier appel a la physique
@@ -93,5 +95,5 @@
 !   temp_air     temperature de l'air 1ere couche
 !   spechum      humidite specifique 1ere couche
-!   hum_air      humidite de l'air
+!   epot_air     temp potentielle de l'air
 !   ccanopy      concentration CO2 canopee
 !   tq_cdrag     cdrag
@@ -102,5 +104,6 @@
 !   precip_rain  precipitation liquide
 !   precip_snow  precipitation solide
-!   lwdown       flux IR entrant a la surface
+!   sollw        flux IR net a la surface
+!   sollwdown    flux IR descendant a la surface
 !   swnet        flux solaire net
 !   swdown       flux solaire entrant a la surface
@@ -142,14 +145,17 @@
   logical, intent(IN) :: debut, lafin, ok_veget
   real, dimension(klon), intent(IN) :: rlon, rlat
+  real, dimension(klon), intent(IN) :: cufi, cvfi
+  real, dimension(klon), intent(INOUT) :: tq_cdrag
   real, dimension(klon), intent(IN) :: zlev
   real, dimension(klon), intent(IN) :: u1_lay, v1_lay
   real, dimension(klon), intent(IN) :: temp_air, spechum
-  real, dimension(klon), intent(IN) :: hum_air, ccanopy
-  real, dimension(klon), intent(IN) :: tq_cdrag, petAcoef, peqAcoef
+  real, dimension(klon), intent(IN) :: epot_air, ccanopy
+  real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
   real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
   real, dimension(klon), intent(IN) :: precip_rain, precip_snow
-  real, dimension(klon), intent(IN) :: lwdown, swnet, swdown, ps, albedo
+  real, dimension(klon), intent(IN) :: sollw, sollwdown, swnet, swdown
+  real, dimension(klon), intent(IN) :: ps, albedo
   real, dimension(klon), intent(IN) :: tsurf, p1lay
-  real, dimension(klon), intent(IN) :: radsol
+  real, dimension(klon), intent(INOUT) :: radsol
   real, dimension(klon), intent(IN) :: zmasq
   real, dimension(klon), intent(IN) :: fder, taux, tauy, rugos, rugoro
@@ -157,5 +163,10 @@
   integer              :: npas, nexca ! nombre et pas de temps couplage
   real, dimension(klon), intent(INOUT) :: evap, snow, qsol
-
+!! PB ajout pour soil
+  logical          :: soil_model
+  integer          :: nsoilmx
+  REAL, DIMENSION(klon, nsoilmx) :: tsoil
+  REAL, dimension(klon)          :: soilcap
+  REAL, dimension(klon)          :: soilflux
 ! Parametres de sortie
   real, dimension(klon), intent(OUT):: fluxsens, fluxlat
@@ -173,5 +184,6 @@
   logical              :: check = .true.
   real, dimension(klon):: cal, beta, dif_grnd, capsol
-  real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=1./86400.*5.
+!!$PB  real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5.
+  real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5.
   real, parameter      :: calsno=1./(2.3867e+06*.15)
   real, dimension(klon):: alb_ice
@@ -180,4 +192,5 @@
   real, dimension(klon):: alb_neig, alb_eau
   real, DIMENSION(klon):: zfra
+  logical              :: cumul = .false.
 
   if (check) write(*,*) 'Entree ', modname
@@ -228,13 +241,37 @@
 !
 ! allocation du run-off
-    if (.not. allocated(run_off)) then
+    if (.not. allocated(coastalflow)) then
+      allocate(coastalflow(knon), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation coastalflow'
+        call abort_gcm(modname,abort_message,1)
+      endif
+      allocate(riverflow(knon), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation riverflow'
+        call abort_gcm(modname,abort_message,1)
+      endif
       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
+        abort_message='Pb allocation runoff'
+        call abort_gcm(modname,abort_message,1)
+      endif
+    else if (size(coastalflow) /= knon) then
       write(*,*)'Bizarre, le nombre de points continentaux'
-      write(*,*)'a change entre deux appels. Je continue ...'
+      write(*,*)'a change entre deux appels. J''arrete ...'
+      abort_message='voir ci-dessus'
+      call abort_gcm(modname,abort_message,1)
+      deallocate(coastalflow, stat = error)
+      allocate(coastalflow(knon), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation coastalflow'
+        call abort_gcm(modname,abort_message,1)
+      endif
+      deallocate(riverflow, stat = error)
+      allocate(riverflow(knon), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation riverflow'
+        call abort_gcm(modname,abort_message,1)
+      endif
       deallocate(run_off, stat = error)
       allocate(run_off(knon), stat = error)
@@ -248,25 +285,9 @@
 ! Calcul age de la neige
 !
-
-  CALL albsno(klon,agesno,alb_neig_grid)  
-  
-  
+!!$ PB ATTENTION changement ordre des appels
+!!$   CALL albsno(klon,agesno,alb_neig_grid)  
+ 
  
     if (.not. ok_veget) then
-!
-! calcul snow et qsol, hydrol adapté
-!
-      call calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
-!      if (check) write(*,*)'Sortie calbeta'
-!      if (check) write(*,*)'RCPD = ',RCPD,' capsol = '
-!      if (check) write(*,*)capsol
-      cal = RCPD * capsol
-      call calcul_fluxs( klon, knon, nisurf, dtime, &
-     &   tsurf, p1lay, cal, beta, tq_cdrag, 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 
@@ -283,15 +304,40 @@
        alb_new = alb_neig*zfra + alb_new*(1.0-zfra)
        z0_new = SQRT(z0_new**2+rugoro**2)
+!
+       CALL albsno(klon,agesno,alb_neig_grid)  
+  
+! calcul snow et qsol, hydrol adapté
+!
+       IF (soil_model) THEN 
+           CALL soil(dtime, nisurf, snow, tsurf, tsoil,soilcap, soilflux)
+           cal = RCPD / soilcap
+           radsol = radsol + soilflux
+       ELSE 
+           CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
+!      if (check) write(*,*)'Sortie calbeta'
+!      if (check) write(*,*)'RCPD = ',RCPD,' capsol = '
+!      if (check) write(*,*)capsol
+           cal = RCPD * capsol
+!!$      cal = capsol
+       ENDIF
+       CALL calcul_fluxs( klon, knon, nisurf, dtime, &
+     &   tsurf, p1lay, cal, beta, tq_cdrag, 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
+        CALL albsno(klon,agesno,alb_neig_grid)  
 !
 !  appel a sechiba
 !
       call interfsol(itime, klon, dtime, nisurf, knon, &
-     &  knindex, rlon, rlat, &
+     &  knindex, rlon, rlat, cufi, cvfi, iim, jjm, pctsrf, &
      &  debut, lafin, ok_veget, &
-     &  zlev,  u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, & 
+     &  zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, & 
      &  tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
-     &  precip_rain, precip_snow, lwdown, swnet, swdown, &
-     &  tsurf, p1lay, ps, radsol, &
+     &  precip_rain, precip_snow, sollwdown, swnet, swdown, &
+     &  tsurf, p1lay/100., ps, radsol, &
      &  evap, fluxsens, fluxlat, &              
      &  tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
@@ -317,9 +363,11 @@
       endif
 
-      call interfoce(itime, dtime, &
+      cumul = .false.
+
+      call interfoce(itime, dtime, cumul, &
       & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
       & ocean, npas, nexca, debut, lafin, &
-      & swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
-      & fder, albedo, taux, tauy, zmasq, &
+      & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
+      & fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
       & tsurf_new, alb_new, alb_ice, pctsrf_new)
 
@@ -345,4 +393,25 @@
      &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
      &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+
+!
+! 2eme appel a interfoce pour le cumul des champs (en particulier
+! fluxsens et fluxlat calcules dans calcul_fluxs)
+!
+    if (ocean == 'couple') then
+
+      cumul = .true.
+
+      call interfoce(itime, dtime, cumul, &
+      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
+      & ocean, npas, nexca, debut, lafin, &
+      & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
+      & fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
+      & tsurf_new, alb_new, alb_ice, pctsrf_new)
+
+!    else if (ocean == 'slab  ') then
+!      call interfoce(nisurf)
+
+    endif
+
 !
 ! calcul albedo
@@ -370,12 +439,15 @@
     if (ocean == 'couple') then
 
-      call interfoce(itime, dtime, &
+      cumul =.false.
+
+      call interfoce(itime, dtime, cumul, &
       & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
       & ocean, npas, nexca, debut, lafin, &
-      & swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
-      & fder, albedo, taux, tauy, zmasq, &
+      & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
+      & fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
       & tsurf_new, alb_new, alb_ice, pctsrf_new)
 
       tsurf_temp = tsurf_new
+      cal = 0.
       dif_grnd = 0.
       beta = 1.0
@@ -390,10 +462,9 @@
 
       tsurf_temp = tsurf
-      dif_grnd = 1.0 / tau_gl
-      beta = 1.0
-    endif
-
-    cal = calice
-    where (snow > 0.0) cal = calsno
+           dif_grnd = 1.0 / tau_gl
+           beta = 1.0
+           cal = RCPD * calice
+           WHERE (snow > 0.0) cal = RCPD * calsno
+    endif
 
     call calcul_fluxs( klon, knon, nisurf, dtime, &
@@ -405,4 +476,23 @@
 
 !
+! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
+!
+    if (ocean == 'couple') then
+
+      cumul =.true.
+
+      call interfoce(itime, dtime, cumul, &
+      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
+      & ocean, npas, nexca, debut, lafin, &
+      & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
+      & fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
+      & tsurf_new, alb_new, alb_ice, pctsrf_new)
+
+!    else if (ocean == 'slab  ') then
+!      call interfoce(nisurf)
+
+    endif
+
+!
 ! calcul albedo
 !
@@ -413,5 +503,5 @@
     alb_new = alb_neig*zfra + 0.6 * (1.0-zfra)
 
-    z0_new = rugos
+    z0_new = 0.001
 
   else if (nisurf == is_lic) then
@@ -423,7 +513,12 @@
 !
 !    call interfsol(nisurf)
-
-    cal = calice
-    where (snow > 0.0) cal = calsno
+    IF (soil_model) THEN 
+        CALL soil(dtime, nisurf, snow, tsurf, tsoil,soilcap, soilflux)
+        cal = RCPD / soilcap
+        radsol = radsol + soilflux
+    ELSE 
+        cal = RCPD * calice
+        WHERE (snow > 0.0) cal = RCPD * calsno
+    ENDIF 
     beta = 1.0
     dif_grnd = 0.0
@@ -490,7 +585,7 @@
 !
   SUBROUTINE interfsol(itime, klon, dtime, nisurf, knon, &
-     & knindex, rlon, rlat, &
+     & knindex, rlon, rlat, cufi, cvfi, iim, jjm, pctsrf, &
      & debut, lafin, ok_veget, &
-     & zlev,  u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, & 
+     & zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, & 
      & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
      & precip_rain, precip_snow, lwdown, swnet, swdown, &
@@ -498,4 +593,6 @@
      & evap, fluxsens, fluxlat, &              
      & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
+
+  USE intersurf
 
 ! Cette routine sert d'interface entre le modele atmospherique et le 
@@ -513,4 +610,5 @@
 !   rlon         longitudes de la grille entiere
 !   rlat         latitudes de la grille entiere
+!   pctsrf       tableau des fractions de surface de chaque maille
 !   debut        logical: 1er appel a la physique (lire les restart)
 !   lafin        logical: dernier appel a la physique (ecrire les restart)
@@ -522,5 +620,5 @@
 !   temp_air     temperature de l'air 1ere couche
 !   spechum      humidite specifique 1ere couche
-!   hum_air      humidite de l'air
+!   epot_air     temp pot de l'air
 !   ccanopy      concentration CO2 canopee
 !   tq_cdrag     cdrag
@@ -531,5 +629,5 @@
 !   precip_rain  precipitation liquide
 !   precip_snow  precipitation solide
-!   lwdown       flux IR entrant a la surface
+!   lwdown       flux IR descendant a la surface
 !   swnet        flux solaire net
 !   swdown       flux solaire entrant a la surface
@@ -560,12 +658,16 @@
   integer, intent(IN) :: nisurf
   integer, intent(IN) :: knon
+  integer, intent(IN) :: iim, jjm
   integer, dimension(klon), intent(IN) :: knindex
   logical, intent(IN) :: debut, lafin, ok_veget
+  real, dimension(klon,nbsrf), intent(IN) :: pctsrf
   real, dimension(klon), intent(IN) :: rlon, rlat
+  real, dimension(klon), intent(IN) :: cufi, cvfi
   real, dimension(klon), intent(IN) :: zlev
   real, dimension(klon), intent(IN) :: u1_lay, v1_lay
   real, dimension(klon), intent(IN) :: temp_air, spechum
-  real, dimension(klon), intent(IN) :: hum_air, ccanopy
-  real, dimension(klon), intent(IN) :: tq_cdrag, petAcoef, peqAcoef
+  real, dimension(klon), intent(IN) :: epot_air, ccanopy
+  real, dimension(klon), intent(INOUT) :: tq_cdrag
+  real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
   real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
   real, dimension(klon), intent(IN) :: precip_rain, precip_snow
@@ -581,5 +683,5 @@
 ! Local
 !
-  integer              :: ii
+  integer              :: ii, ij, jj, igrid, ireal, i, index
   integer              :: error
   character (len = 20) :: modname = 'interfsol'
@@ -595,11 +697,28 @@
 ! pts voisins
   integer,allocatable, dimension(:,:), save :: neighbours
+! fractions continents
+  real,allocatable, dimension(:), save :: contfrac
 ! resolution de la grille
   real, allocatable, dimension (:,:), save :: resolution
+! correspondance point n -> indices (i,j)
+  integer, allocatable, dimension(:,:), save :: correspond
+! offset pour calculer les point voisins
+  integer, dimension(8,3), save :: off_ini
+  integer, dimension(8), save :: offset
 ! Identifieurs des fichiers restart et histoire
   integer, save          :: rest_id, hist_id 
   integer, save          :: rest_id_stom, hist_id_stom
-
+! 
+  real, allocatable, dimension (:,:), save :: lon_scat, lat_scat  
+
+  logical                :: lrestart_read = .true. , lrestart_write = .true.
+
+  real, dimension(klon):: qsurf
   real, dimension(klon):: snow, qsol
+  real                 :: date0 = 0.
+  real, dimension(knon,2) :: albedo_out
+! Pb de nomenclature
+  real, dimension(klon) :: petA_orc, peqA_orc
+  real, dimension(klon) :: petB_orc, peqB_orc
 
   if (check) write(*,*)'Entree ', modname
@@ -607,47 +726,160 @@
 
 ! initialisation
-!    if (debut) then
-!      !
-!      ! Configuration de parametres specifiques a la SSL
-!      !
-!      call intsurf_config(control_in)
-!      !
-!      ! Allouer et initialiser le tableau de coordonnees du sol
-!      !
-!      if (( .not. allocated(lalo))) then
-!        allocate(lalo(knon,2), stat = error)
-!        if (error /= 0) then
-!          abort_message='Pb allocation lalo'
-!          call abort_gcm(modname,abort_message,1)
-!        endif      
-!      endif
-!      do ii = 1, knon
-!        lalo(ii,1) = rlat(knindex(ii))
-!        lalo(ii,2) = rlon(knindex(ii))
-!      enddo
-      !-
-      !- Compute variable to help describe the grid
-      !- once the points are gathered.
-      !-
-!      IF ( (.NOT.ALLOCATED(neighbours))) THEN
-!        ALLOCATE(neighbours(knon,4), stat = error)
-!        if (error /= 0) then
-!          abort_message='Pb allocation neighbours'
-!          call abort_gcm(modname,abort_message,1)
-!        endif
-!      ENDIF
-!      IF ( (.NOT.ALLOCATED(resolution))) THEN
-!        ALLOCATE(resolution(knon,2), stat = error)
-!        if (error /= 0) then
-!          abort_message='Pb allocation resolution'
-!          call abort_gcm(modname,abort_message,1)
-!        endif
-!      ENDIF
-
-! call grid_stuff
-! call sechiba_restart_init
-! call sechiba_history_init
-
-!    endif                          ! (fin debut) 
+  if (debut) then
+
+!
+!  Initialisation des offset    
+!
+! offset bord ouest
+   off_ini(1,1) = - iim  ; off_ini(2,1) = - iim + 1; off_ini(3,1) = 1
+   off_ini(4,1) = iim + 1; off_ini(5,1) = iim      ; off_ini(6,1) = 2 * iim - 1
+   off_ini(7,1) = iim -1 ; off_ini(8,1) = - 1
+! offset point normal
+   off_ini(1,2) = - iim  ; off_ini(2,2) = - iim + 1; off_ini(3,2) = 1
+   off_ini(4,2) = iim + 1; off_ini(5,2) = iim      ; off_ini(6,2) = iim - 1
+   off_ini(7,2) = -1     ; off_ini(8,2) = - iim - 1
+! offset bord   est
+   off_ini(1,3) = - iim; off_ini(2,3) = - 2 * iim + 1; off_ini(3,3) = - iim + 1
+   off_ini(4,3) =  1   ; off_ini(5,3) = iim          ; off_ini(6,3) = iim - 1
+   off_ini(7,3) = -1   ; off_ini(8,3) = - iim - 1
+!
+! Initialisation des correspondances point -> indices i,j
+!
+    if (( .not. allocated(correspond))) then
+      allocate(correspond(iim,jjm+1), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation correspond'
+        call abort_gcm(modname,abort_message,1)
+      endif      
+    endif
+!
+! Attention aux poles
+!
+    do igrid = 1, knon
+      index = knindex(igrid)
+      ij = index - int((index-1)/iim)*iim - 1
+      jj = 2 + int((index-1)/iim)
+      if (mod(index,iim) == 1 ) then
+        jj = 1 + int((index-1)/iim)
+        ij = iim
+      endif
+      correspond(ij,jj) = igrid
+    enddo
+!
+! Allouer et initialiser le tableau de coordonnees du sol
+!
+    if ((.not. allocated(lalo))) then
+      allocate(lalo(knon,2), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation lalo'
+        call abort_gcm(modname,abort_message,1)
+      endif      
+    endif
+    if ((.not. allocated(lon_scat))) then
+      allocate(lon_scat(iim,jjm), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation lon_scat'
+        call abort_gcm(modname,abort_message,1)
+      endif      
+    endif
+    if ((.not. allocated(lat_scat))) then
+      allocate(lat_scat(iim,jjm), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation lat_scat'
+        call abort_gcm(modname,abort_message,1)
+      endif      
+    endif
+    lon_scat = 0.
+    lat_scat = 0.
+    do igrid = 1, knon
+      index = knindex(igrid)
+      lalo(igrid,2) = rlon(index)
+      lalo(igrid,1) = rlat(index)
+      ij = index - int((index-1)/iim)*iim - 1
+      jj = 2 + int((index-1)/iim)
+      if (mod(index,iim) == 1 ) then
+        jj = 1 + int((index-1)/iim)
+        ij = iim
+      endif
+      lon_scat(ij,jj) = rlon(index)
+      lat_scat(ij,jj) = rlat(index)
+    enddo
+    index = 1
+    do jj = 2, jjm
+      do ij = 1, iim
+        index = index + 1
+        lon_scat(ij,jj) = rlon(index)
+        lat_scat(ij,jj) = rlat(index)
+      enddo
+    enddo
+    lon_scat(:,1) = lon_scat(:,2)
+    lat_scat(:,1) = rlat(1)
+
+!
+! Allouer et initialiser le tableau des voisins et des fraction de continents
+!
+    if ( (.not.allocated(neighbours))) THEN
+      allocate(neighbours(knon,8), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation neighbours'
+        call abort_gcm(modname,abort_message,1)
+      endif
+    endif
+    neighbours = 0.
+    if (( .not. allocated(contfrac))) then
+      allocate(contfrac(knon), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation contfrac'
+        call abort_gcm(modname,abort_message,1)
+      endif      
+    endif
+
+    do igrid = 1, knon
+      ireal = knindex(igrid)
+      contfrac(igrid) = pctsrf(ireal,is_ter)
+      if (mod(ireal - 2, iim) == 0) then
+        offset = off_ini(:,1)
+      else if(mod(ireal - 1, iim) == 0) then
+        offset = off_ini(:,3)
+      else
+        offset = off_ini(:,2)
+      endif 
+      if (ireal == 98) write (*,*) offset
+      do i = 1, 8
+        index = ireal + offset(i)
+        if (index <= 1) index = 1
+        if (index >= klon) index = klon
+        if (pctsrf(index, is_ter) > EPSFRA) then
+          ij = index - int((index-1)/iim)*iim - 1
+          jj = 2 + int((index-1)/iim)
+          if (mod(index,iim) == 1 ) then
+            jj = 1 + int((index-1)/iim)
+            ij = iim
+          endif
+!          write(*,*)'correspond',igrid, ireal,index,ij,jj
+          if ( ij >= 1 .and. ij <= iim .and. jj >= 1 .and. jj <= jjm) then
+!          write(*,*)'correspond',igrid, ireal,index,ij,jj
+            neighbours(igrid, i) = correspond(ij, jj)
+          endif
+        endif
+      enddo
+    enddo
+
+!
+!  Allocation et calcul resolutions
+    IF ( (.NOT.ALLOCATED(resolution))) THEN
+      ALLOCATE(resolution(knon,2), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation resolution'
+        call abort_gcm(modname,abort_message,1)
+      endif
+    ENDIF
+    do igrid = 1, knon
+      ij = knindex(igrid)
+      resolution(igrid,1) = cufi(ij)
+      resolution(igrid,2) = cvfi(ij)
+    enddo  
+
+  endif                          ! (fin debut) 
 
 ! 
@@ -655,17 +887,20 @@
 !
 
-!    call sechiba_main(itime, klon, knon, knindex, dtime, &
-!     & debut, lafin, coupling, control_in, &
-!     & lalo, neighbours, resolution,&
-!     & 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, &
-!     & evap, fluxsens, fluxlat, &
-!     & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, &
-!     & rest_id, hist_id, rest_id_stom, hist_id_stom)
-
-!
-! Sauvegarde dans fichiers histoire
-!
+  petA_orc = petBcoef * dtime
+  petB_orc = petAcoef
+  peqA_orc = peqBcoef * dtime
+  peqB_orc = peqAcoef
+
+  call intersurf_main (itime, iim, jjm, knon, knindex, dtime, &
+     & lrestart_read, lrestart_write, lalo, &
+     & contfrac, neighbours, resolution, date0, &
+     & zlev,  u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
+     & tq_cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
+     & precip_rain, precip_snow, lwdown, swnet, swdown, p1lay, &
+     & evap, fluxsens, fluxlat, coastalflow, riverflow, &
+     & tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
+     & lon_scat, lat_scat)
+
+  alb_new(:) = albedo_out(:,1)
 
   END SUBROUTINE interfsol
@@ -673,9 +908,9 @@
 !#########################################################################
 !
-  SUBROUTINE interfoce_cpl(itime, dtime, &
+  SUBROUTINE interfoce_cpl(itime, dtime, cumul, &
       & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
       & ocean, npas, nexca, debut, lafin, &
       & swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
-      & fder, albsol, taux, tauy, zmasq, &
+      & fluxlat, fluxsens, fder, albsol, taux, tauy, zmasq, &
       & tsurf_new, alb_new, alb_ice, pctsrf_new)
 
@@ -712,5 +947,5 @@
 !   nexca        frequence de couplage
 !   swdown       flux solaire entrant a la surface
-!   lwdown       flux IR entrant a la surface
+!   lwdown       flux IR net a la surface
 !   precip_rain  precipitation liquide
 !   precip_snow  precipitation solide
@@ -750,5 +985,6 @@
   INTEGER              :: nexca, npas, kstep
   real, dimension(klon), intent(IN) :: zmasq
-
+  real, dimension(klon), intent(IN) :: fluxlat, fluxsens
+  logical, intent(IN)               :: cumul
   real, dimension(klon), intent(INOUT) :: evap
 
@@ -775,7 +1011,8 @@
   real, dimension(iim, jjm+1) :: wri_sol_ice, wri_sol_sea, wri_nsol_ice 
   real, dimension(iim, jjm+1) :: wri_nsol_sea, wri_fder_ice, wri_evap_ice
-  real, dimension(iim, jjm+1) :: wri_evap_sea
-  real, dimension(iim, jjm+1) :: wri_rain, wri_snow, wri_taux
-  real, dimension(iim, jjm+1) :: wri_tauy, wri_rriv, wri_rcoa
+  REAL, DIMENSION(iim, jjm+1) :: wri_evap_sea, wri_rcoa, wri_rriv
+  REAL, DIMENSION(iim, jjm+1) :: wri_rain, wri_snow, wri_taux, wri_tauy
+  REAL, DIMENSION(iim, jjm+1) :: wri_tauxx, wri_tauyy, wri_tauzz
+  REAL, DIMENSION(iim, jjm+1) :: tmp_lon, tmp_lat
 ! variables relues par le coupleur
 ! read_sic = fraction de glace
@@ -871,30 +1108,34 @@
   cpl_index = 1
   if (nisurf == is_sic) cpl_index = 2
-  do ig = 1, knon
-    cpl_sols(ig,cpl_index) = cpl_sols(ig,cpl_index) &
-     &                               + swdown(ig)      / FLOAT(nexca)
-    cpl_nsol(ig,cpl_index) = cpl_nsol(ig,cpl_index) &
-     &                               + lwdown(ig)      / FLOAT(nexca)
-    cpl_rain(ig,cpl_index) = cpl_rain(ig,cpl_index) &
-     &                               + precip_rain(ig) / FLOAT(nexca)
-    cpl_snow(ig,cpl_index) = cpl_snow(ig,cpl_index) &
-     &                               + precip_snow(ig) / FLOAT(nexca)
-    cpl_evap(ig,cpl_index) = cpl_evap(ig,cpl_index) &
-     &                               + evap(ig)        / FLOAT(nexca)
-    cpl_tsol(ig,cpl_index) = cpl_tsol(ig,cpl_index) &
-     &                               + tsurf(ig)       / FLOAT(nexca)
-    cpl_fder(ig,cpl_index) = cpl_fder(ig,cpl_index) &
-     &                               + fder(ig)        / FLOAT(nexca)
-    cpl_albe(ig,cpl_index) = cpl_albe(ig,cpl_index) &
-     &                               + albsol(ig)      / FLOAT(nexca)
-    cpl_taux(ig,cpl_index) = cpl_taux(ig,cpl_index) &
-     &                               + taux(ig)        / FLOAT(nexca)
-    cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) &
-     &                               + tauy(ig)        / FLOAT(nexca)
-    cpl_rriv(ig,cpl_index) = cpl_rriv(ig,cpl_index) &
-     &                               + 0.     / FLOAT(nexca)/dtime
-    cpl_rcoa(ig,cpl_index) = cpl_rcoa(ig,cpl_index) &
-     &                               + 0.     / FLOAT(nexca)/dtime
-  enddo
+  if (cumul) then
+    do ig = 1, knon
+      if (check) write(*,*) modname, 'cumul des champs'
+      cpl_sols(ig,cpl_index) = cpl_sols(ig,cpl_index) &
+       &                          + swdown(ig)      / FLOAT(nexca)
+      cpl_nsol(ig,cpl_index) = cpl_nsol(ig,cpl_index) &
+       &                          + (lwdown(ig) + fluxlat(ig) +fluxsens(ig))&
+       &                                / FLOAT(nexca)
+      cpl_rain(ig,cpl_index) = cpl_rain(ig,cpl_index) &
+       &                          + precip_rain(ig) / FLOAT(nexca)
+      cpl_snow(ig,cpl_index) = cpl_snow(ig,cpl_index) &
+       &                          + precip_snow(ig) / FLOAT(nexca)
+      cpl_evap(ig,cpl_index) = cpl_evap(ig,cpl_index) &
+       &                          + evap(ig)        / FLOAT(nexca)
+      cpl_tsol(ig,cpl_index) = cpl_tsol(ig,cpl_index) &
+       &                          + tsurf(ig)       / FLOAT(nexca)
+      cpl_fder(ig,cpl_index) = cpl_fder(ig,cpl_index) &
+       &                          + fder(ig)        / FLOAT(nexca)
+      cpl_albe(ig,cpl_index) = cpl_albe(ig,cpl_index) &
+       &                          + albsol(ig)      / FLOAT(nexca)
+      cpl_taux(ig,cpl_index) = cpl_taux(ig,cpl_index) &
+       &                          + taux(ig)        / FLOAT(nexca)
+      cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) &
+       &                          + tauy(ig)        / FLOAT(nexca)
+      cpl_rriv(ig,cpl_index) = cpl_rriv(ig,cpl_index) &
+       &                          + 0.     / FLOAT(nexca)/dtime
+      cpl_rcoa(ig,cpl_index) = cpl_rcoa(ig,cpl_index) &
+       &                          + 0.     / FLOAT(nexca)/dtime
+    enddo
+  endif
 
   if (mod(itime, nexca) == 1) then
@@ -904,5 +1145,5 @@
 ! Si le domaine considere est l'ocean, on lit les champs venant du coupleur
 !
-    if (nisurf == is_oce) then
+    if (nisurf == is_oce .and. .not. cumul) then
       if (check) write(*,*)'rentree fromcpl, itime-1 = ',itime-1
       call fromcpl(itime-1,(jjm+1)*iim,                                  &
@@ -1010,5 +1251,5 @@
 ! Si le domaine considere est la banquise, on envoie les champs au coupleur
 !
-    if (nisurf == is_sic) then
+    if (nisurf == is_sic .and. cumul) then
       wri_rain = 0.; wri_snow = 0.; wri_rcoa = 0.; wri_rriv = 0. 
       wri_taux = 0.; wri_tauy = 0.
@@ -1038,9 +1279,23 @@
       &            tmp_tauy(:,:,2) * tamp_srf(:,:,2) / deno
       endwhere
+!
+! on passe les coordonnées de la grille
+!
+      CALL gath2cpl(rlon(1), tmp_lon(1,1), klon, knon,iim,jjm, knindex) 
+      CALL gath2cpl(rlat(1), tmp_lat(1,1), klon, knon,iim,jjm, knindex) 
+      DO i = 1, iim
+        tmp_lon(i,1) = rlon(i+1)
+        tmp_lon(i,jjm + 1) = rlon(i+1)
+      ENDDO
+!
+! calcul 3 coordonnées du vent
+!
+      CALL atm2geo (iim , jjm + 1, wri_taux, wri_tauy, tmp_lon, tmp_lat, &
+         & wri_tauxx, wri_tauyy, wri_tauzz )
 
       call intocpl(itime, (jjm+1)*iim, wri_sol_ice, wri_sol_sea, wri_nsol_ice,&
       & wri_nsol_sea, wri_fder_ice, wri_evap_ice, wri_evap_sea, wri_rain, &
-      & wri_snow, wri_rcoa, wri_rriv, wri_taux, wri_tauy, wri_taux, wri_tauy, &
-      & lafin ) 
+      & wri_snow, wri_rcoa, wri_rriv, wri_tauxx, wri_tauyy, wri_tauzz, &
+      & wri_tauxx, wri_tauyy, wri_tauzz,lafin ) 
       cpl_sols = 0.; cpl_nsol = 0.; cpl_rain = 0.; cpl_snow = 0.
       cpl_evap = 0.; cpl_tsol = 0.; cpl_fder = 0.; cpl_albe = 0.
@@ -1572,5 +1827,9 @@
   real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
   real                  :: bilan_f, fq_fonte
+  REAL                  :: subli, fsno
   real, parameter :: t_grnd = 271.35, t_coup = 273.15
+!! PB temporaire en attendant mieux pour le modele de neige
+  REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
+!
   logical         :: check = .true.
   character (len = 20)  :: modname = 'calcul_fluxs'
@@ -1578,4 +1837,5 @@
   real            :: max_eau_sol = 150.0
   character (len = 80) :: abort_message 
+  logical,save         :: first = .t.,second=.f.
 
   if (check) write(*,*)'Entree ', modname,' surface = ',nisurf
@@ -1590,4 +1850,14 @@
 ! Traitement neige et humidite du sol
 !
+!    if (first .and. nisurf == is_ter) then
+!      do i = 1, knon
+!        write(67,*)i, snow(i), precip_snow(i), evap(i)
+!      enddo
+!    endif
+!    if (second .and. nisurf == is_ter) then
+!      do i = 1, knon
+!        write(77,*)i, snow(i), precip_snow(i), evap(i)
+!      enddo
+!    endif
     if (nisurf == is_oce) then
       snow = 0.
@@ -1599,4 +1869,5 @@
       qsol = qsol + (precip_rain - evap) * dtime
     endif 
+    IF (nisurf /= is_ter) qsol = max_eau_sol
 
 
@@ -1635,4 +1906,13 @@
      & * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
      & * p1lay(i)/(RD*t1lay(i))
+
+!    if (first .and. nisurf == is_ter) then
+!      write(43,*) &
+!     &i,zx_coef(i),coef1lay(i),u1lay(i),v1lay(i),p1lay(i),t1lay(i)
+!    endif
+!    if (second .and. nisurf == is_ter) then
+!      write(53,*) &
+!     &i,zx_coef(i),coef1lay(i),u1lay(i),v1lay(i),p1lay(i),t1lay(i)
+!    endif
     
   ENDDO
@@ -1664,4 +1944,16 @@
     zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
     zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
+
+!    if (first .and. nisurf == is_ter) then
+!      write(41,*) &
+!     &            i,zx_k1(i),petAcoef(i),petBcoef(i),&
+!     & zx_oh(i),zx_mh(i),zx_nh(i)
+!    endif
+
+!    if (second .and. nisurf == is_ter) then
+!      write(51,*) &
+!     &            i,zx_k1(i),petAcoef(i),petBcoef(i),&
+!     & zx_oh(i),zx_mh(i),zx_nh(i)
+!    endif
 
 ! Tsurface
@@ -1691,4 +1983,16 @@
     dflux_s(i) = zx_nh(i)
     dflux_l(i) = (zx_sl(i) * zx_nq(i))
+
+!    if (first .and. nisurf == is_ter) then
+!      write(42,*) &
+!     &            i,tsurf_new(i),d_ts(i),zx_h_ts(i),zx_q_0(i),&
+!     & evap(i),fluxsens(i)
+!    endif
+
+!    if (second .and. nisurf == is_ter) then
+!      write(52,*) &
+!     &            i,tsurf_new(i),d_ts(i),zx_h_ts(i),zx_q_0(i),&            
+!     & evap(i),fluxsens(i) 
+!    endif
 !
 ! en cas de fonte de neige
@@ -1707,4 +2011,13 @@
     qsol(i) = min(qsol(i), max_eau_sol) 
   ENDDO
+
+!  if (nisurf == is_oce .and. second) then
+!    second =.f.
+!  endif
+
+!  if (nisurf == is_oce .and. first) then
+!    first = .f.
+!    second = .t.
+!  endif
 
   END SUBROUTINE calcul_fluxs
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.F	(revision 177)
@@ -80,8 +80,14 @@
       cl_writ(10)='CORUNCOA'
       cl_writ(11)='CORIVFLU'
-      cl_writ(12)='COZOTAUX'
-      cl_writ(13)='COZOTAUV'
-      cl_writ(14)='COMETAUY'
-      cl_writ(15)='COMETAUU'
+c$$$      cl_writ(12)='COZOTAUX'
+c$$$      cl_writ(13)='COZOTAUV'
+c$$$      cl_writ(14)='COMETAUY'
+c$$$      cl_writ(15)='COMETAUU'
+      cl_writ(12)='COTAUXXU'
+      cl_writ(13)='COTAUYYU'
+      cl_writ(14)='COTAUZZU'
+      cl_writ(15)='COTAUXXV'
+      cl_writ(16)='COTAUYYV'
+      cl_writ(17)='COTAUZZV'
 c
 c     Define files name for fields exchanged from atmos to coupler,
@@ -103,5 +109,7 @@
       cl_f_writ(14)='flxatmos'
       cl_f_writ(15)='flxatmos'
-c      cl_f_writ(16)='flxatmos'
+      cl_f_writ(16)='flxatmos'
+      cl_f_writ(17)='flxatmos'
+
 c
 c
@@ -294,5 +302,5 @@
       SUBROUTINE intocpl(kt, imjm, fsolice, fsolwat, fnsolice, fnsolwat,
      $    fnsicedt, evice, evwat, lpre, spre, dirunoff, rivrunoff,
-     $    tauxu, tauxv, tauyv, tauyu, last) 
+     $    tauxx_u, tauyy_u, tauzz_u, tauxx_v, tauyy_v, tauzz_v,last) 
 c ======================================================================
 c S. Valcke (02/99) adapted From L.Z.X Li: this subroutine provides the 
@@ -306,8 +314,7 @@
       REAL fsolice(imjm)
       REAL fsolwat(imjm)
-      REAL fnsolice(imjm)
       REAL fnsolwat(imjm)
-      REAL fnsicedt(imjm)
-      REAL ictemp(imjm)
+      REAL fnsolice(imjm) 
+      REAL fnsicedt(imjm) 
       REAL evice(imjm)
       REAL evwat(imjm)
@@ -316,8 +323,14 @@
       REAL dirunoff(imjm)
       REAL rivrunoff(imjm)
-      REAL tauxu(imjm)
-      REAL tauxv(imjm)
-      REAL tauyu(imjm)
-      REAL tauyv(imjm)
+c$$$      REAL tauxu(imjm)
+c$$$      REAL tauxv(imjm)
+c$$$      REAL tauyu(imjm)
+c$$$      REAL tauyv(imjm)
+      REAL tauxx_u(imjm)
+      REAL tauxx_v(imjm)
+      REAL tauyy_u(imjm)
+      REAL tauyy_v(imjm)
+      REAL tauzz_u(imjm)
+      REAL tauzz_v(imjm)
       LOGICAL last
 c
@@ -431,16 +444,34 @@
      $          CALL locwrite(cl_writ(jf),rivrunoff, imjm,
      $          file_unit_field(jf), ierror, nuout) 
+c$$$            IF (jf.eq.12)
+c$$$     $          CALL locwrite(cl_writ(jf),tauxu, imjm,
+c$$$     $          file_unit_field(jf),ierror, nuout) 
+c$$$            IF (jf.eq.13)
+c$$$     $          CALL locwrite(cl_writ(jf),tauxv, imjm,
+c$$$     $          file_unit_field(jf),ierror, nuout) 
+c$$$            IF (jf.eq.14)
+c$$$     $          CALL locwrite(cl_writ(jf),tauyv, imjm,
+c$$$     $          file_unit_field(jf),ierror, nuout) 
+c$$$            IF (jf.eq.15)
+c$$$     $          CALL locwrite(cl_writ(jf),tauyu, imjm,
+c$$$     $          file_unit_field(jf), ierror, nuout) 
             IF (jf.eq.12)
-     $          CALL locwrite(cl_writ(jf),tauxu, imjm,
-     $          file_unit_field(jf),ierror, nuout) 
+     $          CALL locwrite(cl_writ(jf),tauxx_u, imjm,
+     $          file_unit_field(jf),ierror)
             IF (jf.eq.13)
-     $          CALL locwrite(cl_writ(jf),tauxv, imjm,
-     $          file_unit_field(jf),ierror, nuout) 
+     $          CALL locwrite(cl_writ(jf),tauyy_u, imjm,
+     $          file_unit_field(jf),ierror)
             IF (jf.eq.14)
-     $          CALL locwrite(cl_writ(jf),tauyv, imjm,
-     $          file_unit_field(jf),ierror, nuout) 
+     $          CALL locwrite(cl_writ(jf),tauzz_u, imjm,
+     $          file_unit_field(jf),ierror)
             IF (jf.eq.15)
-     $          CALL locwrite(cl_writ(jf),tauyu, imjm,
-     $          file_unit_field(jf), ierror, nuout) 
+     $          CALL locwrite(cl_writ(jf),tauxx_v, imjm,
+     $          file_unit_field(jf),ierror)
+            IF (jf.eq.16)
+     $          CALL locwrite(cl_writ(jf),tauyy_v, imjm,
+     $          file_unit_field(jf),ierror)
+            IF (jf.eq.17)
+     $          CALL locwrite(cl_writ(jf),tauzz_v, imjm,
+     $          file_unit_field(jf),ierror)
           END DO 
 C
@@ -484,8 +515,14 @@
           IF (jn.eq.10) CALL CLIM_Export(cl_writ(jn),kt,dirunoff, info)
           IF (jn.eq.11) CALL CLIM_Export(cl_writ(jn),kt,rivrunoff,info)
-          IF (jn.eq.12) CALL CLIM_Export(cl_writ(jn), kt, tauxu, info)
-          IF (jn.eq.13) CALL CLIM_Export(cl_writ(jn), kt, tauxv, info)
-          IF (jn.eq.14) CALL CLIM_Export(cl_writ(jn), kt, tauyv, info)
-          IF (jn.eq.15) CALL CLIM_Export(cl_writ(jn), kt, tauyu, info)
+c$$$          IF (jn.eq.12) CALL CLIM_Export(cl_writ(jn), kt, tauxu, info)
+c$$$          IF (jn.eq.13) CALL CLIM_Export(cl_writ(jn), kt, tauxv, info)
+c$$$          IF (jn.eq.14) CALL CLIM_Export(cl_writ(jn), kt, tauyv, info)
+c$$$          IF (jn.eq.15) CALL CLIM_Export(cl_writ(jn), kt, tauyu, info)
+          IF (jn.eq.12) CALL CLIM_Export(cl_writ(jn), kt, tauxx_u, info)
+          IF (jn.eq.13) CALL CLIM_Export(cl_writ(jn), kt, tauyy_u, info)
+          IF (jn.eq.14) CALL CLIM_Export(cl_writ(jn), kt, tauzz_u, info)
+          IF (jn.eq.15) CALL CLIM_Export(cl_writ(jn), kt, tauxx_v, info)
+          IF (jn.eq.16) CALL CLIM_Export(cl_writ(jn), kt, tauyy_v, info)
+          IF (jn.eq.17) CALL CLIM_Export(cl_writ(jn), kt, tauzz_v, info)
           
             IF (info .NE. CLIM_Ok) THEN
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.true
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.true	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/oasis.true	(revision 177)
@@ -80,8 +80,14 @@
       cl_writ(10)='CORUNCOA'
       cl_writ(11)='CORIVFLU'
-      cl_writ(12)='COZOTAUX'
-      cl_writ(13)='COZOTAUV'
-      cl_writ(14)='COMETAUY'
-      cl_writ(15)='COMETAUU'
+c$$$      cl_writ(12)='COZOTAUX'
+c$$$      cl_writ(13)='COZOTAUV'
+c$$$      cl_writ(14)='COMETAUY'
+c$$$      cl_writ(15)='COMETAUU'
+      cl_writ(12)='COTAUXXU'
+      cl_writ(13)='COTAUYYU'
+      cl_writ(14)='COTAUZZU'
+      cl_writ(15)='COTAUXXV'
+      cl_writ(16)='COTAUYYV'
+      cl_writ(17)='COTAUZZV'
 c
 c     Define files name for fields exchanged from atmos to coupler,
@@ -103,5 +109,7 @@
       cl_f_writ(14)='flxatmos'
       cl_f_writ(15)='flxatmos'
-c      cl_f_writ(16)='flxatmos'
+      cl_f_writ(16)='flxatmos'
+      cl_f_writ(17)='flxatmos'
+
 c
 c
@@ -294,5 +302,5 @@
       SUBROUTINE intocpl(kt, imjm, fsolice, fsolwat, fnsolice, fnsolwat,
      $    fnsicedt, evice, evwat, lpre, spre, dirunoff, rivrunoff,
-     $    tauxu, tauxv, tauyv, tauyu, last) 
+     $    tauxx_u, tauyy_u, tauzz_u, tauxx_v, tauyy_v, tauzz_v,last) 
 c ======================================================================
 c S. Valcke (02/99) adapted From L.Z.X Li: this subroutine provides the 
@@ -306,8 +314,7 @@
       REAL fsolice(imjm)
       REAL fsolwat(imjm)
-      REAL fnsolice(imjm)
       REAL fnsolwat(imjm)
-      REAL fnsicedt(imjm)
-      REAL ictemp(imjm)
+      REAL fnsolice(imjm) 
+      REAL fnsicedt(imjm) 
       REAL evice(imjm)
       REAL evwat(imjm)
@@ -316,8 +323,14 @@
       REAL dirunoff(imjm)
       REAL rivrunoff(imjm)
-      REAL tauxu(imjm)
-      REAL tauxv(imjm)
-      REAL tauyu(imjm)
-      REAL tauyv(imjm)
+c$$$      REAL tauxu(imjm)
+c$$$      REAL tauxv(imjm)
+c$$$      REAL tauyu(imjm)
+c$$$      REAL tauyv(imjm)
+      REAL tauxx_u(imjm)
+      REAL tauxx_v(imjm)
+      REAL tauyy_u(imjm)
+      REAL tauyy_v(imjm)
+      REAL tauzz_u(imjm)
+      REAL tauzz_v(imjm)
       LOGICAL last
 c
@@ -431,16 +444,34 @@
      $          CALL locwrite(cl_writ(jf),rivrunoff, imjm,
      $          file_unit_field(jf), ierror, nuout) 
+c$$$            IF (jf.eq.12)
+c$$$     $          CALL locwrite(cl_writ(jf),tauxu, imjm,
+c$$$     $          file_unit_field(jf),ierror, nuout) 
+c$$$            IF (jf.eq.13)
+c$$$     $          CALL locwrite(cl_writ(jf),tauxv, imjm,
+c$$$     $          file_unit_field(jf),ierror, nuout) 
+c$$$            IF (jf.eq.14)
+c$$$     $          CALL locwrite(cl_writ(jf),tauyv, imjm,
+c$$$     $          file_unit_field(jf),ierror, nuout) 
+c$$$            IF (jf.eq.15)
+c$$$     $          CALL locwrite(cl_writ(jf),tauyu, imjm,
+c$$$     $          file_unit_field(jf), ierror, nuout) 
             IF (jf.eq.12)
-     $          CALL locwrite(cl_writ(jf),tauxu, imjm,
-     $          file_unit_field(jf),ierror, nuout) 
+     $          CALL locwrite(cl_writ(jf),tauxx_u, imjm,
+     $          file_unit_field(jf),ierror)
             IF (jf.eq.13)
-     $          CALL locwrite(cl_writ(jf),tauxv, imjm,
-     $          file_unit_field(jf),ierror, nuout) 
+     $          CALL locwrite(cl_writ(jf),tauyy_u, imjm,
+     $          file_unit_field(jf),ierror)
             IF (jf.eq.14)
-     $          CALL locwrite(cl_writ(jf),tauyv, imjm,
-     $          file_unit_field(jf),ierror, nuout) 
+     $          CALL locwrite(cl_writ(jf),tauzz_u, imjm,
+     $          file_unit_field(jf),ierror)
             IF (jf.eq.15)
-     $          CALL locwrite(cl_writ(jf),tauyu, imjm,
-     $          file_unit_field(jf), ierror, nuout) 
+     $          CALL locwrite(cl_writ(jf),tauxx_v, imjm,
+     $          file_unit_field(jf),ierror)
+            IF (jf.eq.16)
+     $          CALL locwrite(cl_writ(jf),tauyy_v, imjm,
+     $          file_unit_field(jf),ierror)
+            IF (jf.eq.17)
+     $          CALL locwrite(cl_writ(jf),tauzz_v, imjm,
+     $          file_unit_field(jf),ierror)
           END DO 
 C
@@ -484,8 +515,14 @@
           IF (jn.eq.10) CALL CLIM_Export(cl_writ(jn),kt,dirunoff, info)
           IF (jn.eq.11) CALL CLIM_Export(cl_writ(jn),kt,rivrunoff,info)
-          IF (jn.eq.12) CALL CLIM_Export(cl_writ(jn), kt, tauxu, info)
-          IF (jn.eq.13) CALL CLIM_Export(cl_writ(jn), kt, tauxv, info)
-          IF (jn.eq.14) CALL CLIM_Export(cl_writ(jn), kt, tauyv, info)
-          IF (jn.eq.15) CALL CLIM_Export(cl_writ(jn), kt, tauyu, info)
+c$$$          IF (jn.eq.12) CALL CLIM_Export(cl_writ(jn), kt, tauxu, info)
+c$$$          IF (jn.eq.13) CALL CLIM_Export(cl_writ(jn), kt, tauxv, info)
+c$$$          IF (jn.eq.14) CALL CLIM_Export(cl_writ(jn), kt, tauyv, info)
+c$$$          IF (jn.eq.15) CALL CLIM_Export(cl_writ(jn), kt, tauyu, info)
+          IF (jn.eq.12) CALL CLIM_Export(cl_writ(jn), kt, tauxx_u, info)
+          IF (jn.eq.13) CALL CLIM_Export(cl_writ(jn), kt, tauyy_u, info)
+          IF (jn.eq.14) CALL CLIM_Export(cl_writ(jn), kt, tauzz_u, info)
+          IF (jn.eq.15) CALL CLIM_Export(cl_writ(jn), kt, tauxx_v, info)
+          IF (jn.eq.16) CALL CLIM_Export(cl_writ(jn), kt, tauyy_v, info)
+          IF (jn.eq.17) CALL CLIM_Export(cl_writ(jn), kt, tauzz_v, info)
           
             IF (info .NE. CLIM_Ok) THEN
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/orografi.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/orografi.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/orografi.F	(revision 177)
@@ -220,4 +220,5 @@
 c
       real   ztau(klon,klev+1),
+     $       ztauf(klon,klev+1),
      *       zstab(klon,klev+1),
      *       zvph(klon,klev+1),
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/phyredem.F	(revision 177)
@@ -648,4 +648,20 @@
 #endif
 c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "RUGMER", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "RUGMER", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 28,
+     .                        "Longueur de rugosite sur mer")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,frugs(1,is_oce))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,frugs(1,is_oce))
+#endif
+c
+c
       ierr = NF_CLOSE(nid)
 c
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 177)
@@ -3,7 +3,9 @@
      .            paprs,pplay,pphi,pphis,paire,presnivs,clesphy0,
      .            u,v,t,qx,
-     .            omega,
+     .            omega, cufi, cvfi,
      .            d_u, d_v, d_t, d_qx, d_ps)
       USE ioipsl
+      USE histcom
+
       IMPLICIT none
 c======================================================================
@@ -89,5 +91,5 @@
 ccc      PARAMETER (soil_model=.FALSE.)
       logical ok_veget
-      parameter (ok_veget = .false.)
+      parameter (ok_veget = .true.)
 c======================================================================
 c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
@@ -217,4 +219,6 @@
       REAL fevap(klon,nbsrf)
       SAVE fevap                 ! evaporation
+      REAL fluxlat(klon,nbsrf)
+      SAVE fluxlat
 c
       REAL deltat(klon)
@@ -339,5 +343,4 @@
 cAA
       EXTERNAL hgardfou  ! verifier les temperatures
-      EXTERNAL hydrol    ! hydrologie du sol
       EXTERNAL nuage     ! calculer les proprietes radiatives
       EXTERNAL o3cm      ! initialiser l'ozone
@@ -379,9 +382,10 @@
       REAL cool0(klon,klev)   ! refroidissement infrarouge ciel clair
       REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
+      real sollwdown(klon)    ! downward LW flux at surface
       REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
       REAL albpla(klon)
 c Le rayonnement n'est pas calcule tous les pas, il faut donc
 c                      sauvegarder les sorties du rayonnement
-      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw
+      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
       SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
       INTEGER itaprad
@@ -492,4 +496,5 @@
       logical ok_sync
 
+
 c
 c Declaration des constantes et des fonctions thermodynamiques
@@ -538,5 +543,5 @@
          CALL phyetat0 ("startphy.nc",dtime,co2_ppm,solaire,
      .       rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsol,fsnow,
-     .       falbe, fevap, rain_fall,snow_fall,solsw, sollw,
+     .       falbe, fevap, rain_fall,snow_fall,solsw, sollwdown,
      .       fder,radsol,frugs,agesno,clesphy0,
      .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
@@ -718,4 +723,8 @@
      .                "ave(X)", zsto,zout)
 c
+         CALL histdef(nid_day, "solldown", "Down. IR rad. at surface", 
+     .                "W/m2", iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "ave(X)", zsto,zout)
+c
          CALL histdef(nid_day, "bils", "Surf. total heat flux", "W/m2",
      .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
@@ -746,5 +755,10 @@
      $         iim,jjmp1,nhori, 1,1,1, -99, 32,
      $         "ave(X)", zsto,zout)
-
+C
+           call histdef(nid_day, "tsol_"//clnsurf(nsrf), 
+     $         "Fraction"//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+C
            call histdef(nid_day, "sens_"//clnsurf(nsrf), 
      $         "Sensible heat flux "//clnsurf(nsrf), "W/m2",  
@@ -766,4 +780,15 @@
      $         iim,jjmp1,nhori, 1,1,1, -99, 32,
      $         "ave(X)", zsto,zout)
+C
+           call histdef(nid_day, "albe_"//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, "rugs_"//clnsurf(nsrf), 
+     $         "Latent heat flux "//clnsurf(nsrf), "W/m2",  
+     $         iim,jjmp1,nhori, 1,1,1, -99, 32,
+     $         "ave(X)", zsto,zout)
+
 C§§§
          END DO 
@@ -921,4 +946,8 @@
      .                "ave(X)", zsto,zout)
 c
+         CALL histdef(nid_mth, "solldown", "Down. IR rad. at surface", 
+     .                "W/m2", iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "ave(X)", zsto,zout)
+c
          CALL histdef(nid_mth, "tops0", "Solar rad. at TOA", "W/m2",
      .                iim,jjmp1,nhori, 1,1,1, -99, 32,
@@ -964,4 +993,9 @@
      $         "ave(X)", zsto,zout)
 C
+           call histdef(nid_mth, "tsol_"//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",  
@@ -981,4 +1015,14 @@
            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)
+c
+           call histdef(nid_mth, "albe_"//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, "rugs_"//clnsurf(nsrf), 
+     $         "Latent heat flux "//clnsurf(nsrf), "W/m2",  
      $         iim,jjmp1,nhori, 1,1,1, -99, 32,
      $         "ave(X)", zsto,zout)
@@ -1258,4 +1302,8 @@
      .                "inst(X)", zsto,zout)
 c
+         CALL histdef(nid_ins, "solldown", "Down. IR rad. at surface", 
+     .                "W/m2", iim,jjmp1,nhori, 1,1,1, -99, 32, 
+     .                "ave(X)", zsto,zout)
+c
          CALL histdef(nid_ins, "bils", "Surf. total heat flux", "W/m2",
      .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
@@ -1526,5 +1574,5 @@
       DO nsrf = 1, nbsrf
       DO i = 1, klon
-         zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
+            zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
       ENDDO
       ENDDO
@@ -1546,6 +1594,7 @@
      e            julien, rmu0, 
      e            ok_veget, ocean, npas, nexca, ftsol,
-     e            paprs,pplay,radsol, fsnow,fqsol,fevap,falbe,
-     e            rain_fall, snow_fall, solsw, sollw, fder,
+     $            soil_model,ftsoil,
+     $            paprs,pplay,radsol, fsnow,fqsol,fevap,falbe,fluxlat,
+     e            rain_fall, snow_fall, solsw, sollw, sollwdown, fder,
      e            rlon, rlat, cufi, cvfi, frugs,
      e            debut, lafin, agesno,rugoro ,
@@ -1605,6 +1654,8 @@
       DO nsrf = 1, nbsrf
       DO i = 1, klon
-         ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
-         zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
+c$$$        IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN 
+            ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
+            zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
+c$$$        ENDIF 
       ENDDO
       ENDDO
@@ -1668,4 +1719,6 @@
      s            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
      s            kcbot, kctop, kdtop, pmflxr, pmflxs)
+      WHERE (rain_con < 0.) rain_con = 0.
+      WHERE (snow_con < 0.) snow_con = 0.
       DO i = 1, klon
          ibas_con(i) = klev+1 - kcbot(i)
@@ -1758,4 +1811,6 @@
      .           frac_impa, frac_nucl,
      .           prfl, psfl)
+      WHERE (rain_lsc < 0) rain_lsc = 0.
+      WHERE (snow_lsc < 0) snow_lsc = 0.
       DO k = 1, klev
       DO i = 1, klon
@@ -1861,4 +1916,5 @@
      s             heat,heat0,cool,cool0,radsol,albpla,
      s             topsw,toplw,solsw,sollw,
+     s             sollwdown,
      s             topsw0,toplw0,solsw0,sollw0)
       itaprad = 0
@@ -1893,12 +1949,12 @@
 c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
 c
-      DO nsrf = 1, nbsrf
-      DO i = 1, klon
-         IF (pctsrf(i,nsrf).LT.epsfra) THEN
-            fqsol(i,nsrf) = zxqsol(i)
-            fsnow(i,nsrf) = zxsnow(i)
-         ENDIF
-      ENDDO
-      ENDDO
+c$$$      DO nsrf = 1, nbsrf
+c$$$      DO i = 1, klon
+c$$$         IF (pctsrf(i,nsrf).LT.epsfra) THEN
+c$$$            fqsol(i,nsrf) = zxqsol(i)
+c$$$            fsnow(i,nsrf) = zxsnow(i)
+c$$$         ENDIF
+c$$$      ENDDO
+c$$$      ENDDO
 c
 c Calculer le bilan du sol et la derive de temperature (couplage)
@@ -2092,4 +2148,7 @@
       CALL histwrite(nid_day,"soll",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollwdown,zx_tmp_2d)
+      CALL histwrite(nid_day,"solldown",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
       CALL histwrite(nid_day,"bils",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
@@ -2108,24 +2167,39 @@
         CALL histwrite(nid_day,"pourc_"//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_day,"tsol_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
 C 
-        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
+        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)
+        zx_tmp_fi2d(1 : klon) = fluxlat( 1 : klon, 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)
+        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)
+        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
+        zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_day,"albe_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
+C
+        zx_tmp_fi2d(1 : klon) = frugs( 1 : klon, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_day,"rugs_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d) 
 C 
       END DO  
@@ -2258,4 +2332,7 @@
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
       CALL histwrite(nid_mth,"soll",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollwdown,zx_tmp_2d)
+      CALL histwrite(nid_mth,"solldown",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw0,zx_tmp_2d)
@@ -2300,24 +2377,39 @@
      $      zx_tmp_2d,iim*jjmp1,ndex2d) 
 C 
-        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
+        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_mth,"tsol_"//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)
+     $      zx_tmp_2d,iim*jjmp1,ndex2d)
+C
+        zx_tmp_fi2d(1 : klon) = fluxlat( 1 : klon, 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)
+        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)
+        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 
+        zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_mth,"albe_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d)
+C
+        zx_tmp_fi2d(1 : klon) = frugs( 1 : klon, nsrf)
+        CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
+        CALL histwrite(nid_mth,"rugs_"//clnsurf(nsrf),itap,
+     $      zx_tmp_2d,iim*jjmp1,ndex2d)
+
       END DO  
 c$$$      DO i = 1, klon
@@ -2559,4 +2651,7 @@
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
       CALL histwrite(nid_ins,"soll",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
+c
+      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollwdown,zx_tmp_2d)
+      CALL histwrite(nid_ins,"solldown",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
 c
       CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
@@ -2588,10 +2683,10 @@
      $      zx_tmp_2d,iim*jjmp1,ndex2d) 
 C 
-        zx_tmp_fi2d(1 : klon) = - fluxt( 1 : klon, 1, nsrf)
+        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)
+        zx_tmp_fi2d(1 : klon) = fluxlat( 1 : klon, nsrf)
         CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d , zx_tmp_2d)
         CALL histwrite(nid_ins,"lat_"//clnsurf(nsrf),itap,
@@ -2603,10 +2698,10 @@
      $      zx_tmp_2d,iim*jjmp1,ndex2d) 
 C
-        zx_tmp_fi2d(1 : klon) = - fluxu( 1 : klon, 1, nsrf)
+        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)
+        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,
@@ -2752,5 +2847,5 @@
      .      rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsol, fsnow,
      .      falbe, fevap, rain_fall, snow_fall,
-     .      solsw, sollw,fder,
+     .      solsw, sollwdown,fder,
      .      radsol,frugs,agesno,
      .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/phystokenc.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/phystokenc.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/phystokenc.F	(revision 177)
@@ -7,4 +7,5 @@
      O                   physid)
       USE ioipsl
+      USE histcom
 
       IMPLICIT none
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/phytrac.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/phytrac.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/phytrac.F	(revision 177)
@@ -10,4 +10,6 @@
      O                   tr_seri)
       USE ioipsl
+      USE histcom
+
       IMPLICIT none
 c======================================================================
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F	(revision 177)
@@ -4,4 +4,5 @@
      .                  heat,heat0,cool,cool0,radsol,albpla,
      .                  topsw,toplw,solsw,sollw,
+     .                  sollwdown,
      .                  topsw0,toplw0,solsw0,sollw0)
       IMPLICIT none
@@ -52,4 +53,5 @@
       real solsw(klon), sollw(klon), albpla(klon)
       real topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
+      real sollwdown(klon)
 c
       REAL*8 zx_alpha1, zx_alpha2
@@ -83,4 +85,5 @@
       REAL*8 ztopsw(kdlon), ztoplw(kdlon)
       REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
+      REAL*8 zsollwdown(kdlon)
       REAL*8 ztopsw0(kdlon), ztoplw0(kdlon)
       REAL*8 zsolsw0(kdlon), zsollw0(kdlon)
@@ -186,5 +189,6 @@
      .        PVIEW,
      .        zcool, zcool0,
-     .        ztoplw,zsollw,ztoplw0,zsollw0)
+     .        ztoplw,zsollw,ztoplw0,zsollw0,
+     .        zsollwdown)
       CALL SW(PSCT, RCO2, zrmu0, zfract,
      S        PPMB, PDP,
@@ -201,4 +205,5 @@
          solsw(iof+i) = zsolsw(i)
          sollw(iof+i) = zsollw(i)
+         sollwdown(iof+i) = zsollwdown(i)
          topsw0(iof+i) = ztopsw0(i)
          toplw0(iof+i) = ztoplw0(i)
@@ -2454,5 +2459,6 @@
      .              PVIEW,
      .              PCOLR, PCOLR0,
-     .              PTOPLW,PSOLLW,PTOPLW0,PSOLLW0)
+     .              PTOPLW,PSOLLW,PTOPLW0,PSOLLW0,
+     .              psollwdown)
       IMPLICIT none
 #include "dimensions.h"
@@ -2516,4 +2522,6 @@
       REAL*8 PTOPLW0(KDLON)       ! LONGWAVE FLUX AT T.O.A. (CLEAR-SKY)
       REAL*8 PSOLLW0(KDLON)       ! LONGWAVE FLUX AT SURFACE (CLEAR-SKY)
+c Rajout LF
+      real*8 psollwdown(kdlon)    ! LONGWAVE downwards flux at surface
 C
 C-------------------------------------------------------------------------
@@ -2588,4 +2596,5 @@
          PSOLLW0(i) = -ZFLUC(i,1,1)-ZFLUC(i,2,1)
          PTOPLW0(i) = ZFLUC(i,1,KFLEV+1) + ZFLUC(i,2,KFLEV+1)
+         psollwdown(i) = -ZFLUX(i,2,1)
       ENDDO
 C     ------------------------------------------------------------------
Index: LMDZ.3.3/branches/rel-LF/liste_des_sources_f90
===================================================================
--- LMDZ.3.3/branches/rel-LF/liste_des_sources_f90	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/liste_des_sources_f90	(revision 177)
@@ -11,4 +11,6 @@
       USE IOIPSL
 libf/dyn3d/abort_gcm.F
+      USE ioipsl
+libf/dyn3d/create_limit.F
       USE IOIPSL
 libf/dyn3d/dynredem.F
Index: LMDZ.3.3/branches/rel-LF/makegcm
===================================================================
--- LMDZ.3.3/branches/rel-LF/makegcm	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/makegcm	(revision 177)
@@ -168,8 +168,9 @@
 else if $SUN then
    set optim=" -fast "
+#   set optim=" -g "
    set optimbis=" "
    set optim90=" -fast -fixed "
-   set optimtru90=" -fast -free"
-   set opt_link="-lf77compat -L$IOIPSLDIR -lioipsl  -L$NCDFLIB -lnetcdf"
+   set optimtru90=" -fast -free "
+   set opt_link="-lf77compat -L$modipsl/lib -lsechiba -lparameters -lstomate -L$IOIPSLDIR -lioipsl -L$NCDFLIB -lnetcdf "
    set mod_loc_dir=$localdir
    set mod_suffix=mod
@@ -199,7 +200,7 @@
 else if $LINUX then
    set optim="-fast "
-   set optim90=" -fast -module $libo"
-   set optimtru90=" -fast -c -Mfree -module /d3/fairhead/sechiba/ioipsl"
-   set opt_link=" -Mfree -L/usr/local/pgi/linux86/lib -lpgf90 -lpgftnrtl -lpghpf -lpghpf2 -L$NCDFLIB -lnetcdf -L$IOIPSLDIR -lioipsl -Wl,-Bstatic "
+   set optim90=" -fast "
+   set optimtru90=" -fast -c -Mfree -module $IOIPSLDIR "
+   set opt_link=" -Mfree -L /usr/local/pgi/linux86/lib -lpgf90 -lpgftnrtl -lpghpf -lpghpf2 -L$modipsl/lib -lsechiba -lparameters -lstomate -L$NCDFLIB -lnetcdf -L$IOIPSLDIR -lioipsl -Wl,-Bstatic -L/usr/lib/gcc-lib/i386-linux/2.95.2/"
    set mod_loc_dir=$IOIPSLDIR
    set mod_suffix=mod
@@ -390,12 +391,15 @@
            setenv PARALLEL 2
            set optim=" -g -C "
-           set optim90=" -fixed -g "
-           set optimtru90=" -free -g "
+           set optim90=" -fixed -g -C "
+           set optimtru90=" -free -g -C "
+#           set optim=" -g "
+#           set optim90=" -fixed -g "
+#           set optimtru90=" -free -g "
         else if $CRAY then
            set optim="$optim"" -g "
            set optim90="$optim90"" -G1 "
         else if $LINUX then
-           set optim="$optim"" -g -Mbounds "
-           set optim90="$optim90"" -g -Mbounds "
+           set optim="$optim"" -g -Mbounds -C "
+           set optim90="$optim90"" -g -Mbounds -C "
         else 
            echo pas d option debug predefinie pour cette machine
@@ -715,10 +719,7 @@
  \cp $IOIPSLDIR/*.mod $libo
 else if $SUN then
- set optim90=" $optim90 -M$libo "
- set optimtru90=" $optimtru90 -M$libo "
+ set optim90=" $optim90 -M$libo -M$modipsl/lib "
+ set optimtru90=" $optimtru90 -M$libo -M$modipsl/lib "
  set optim="$optim90"
- \cp /d3/fairhead/sechiba_sun/parameters/*.mod $libo
- \cp /d3/fairhead/sechiba_sun/sechiba/*.mod $libo
- \cp /d3/fairhead/sechiba_sun/stomate/*.mod $libo
  \cp $IOIPSLDIR/*.mod $libo
 else if $NEC then
@@ -728,5 +729,5 @@
  set optimtru90=" $optimtru90 -I$libo "
 else if $LINUX then
- set optim90=" -fast -module $libo "
+ set optim90=" $optim90 -module $libo "
  set optim="$optim90"
  set mod_loc_dir=$libo
Index: LMDZ.3.3/branches/rel-LF/tmp
===================================================================
--- LMDZ.3.3/branches/rel-LF/tmp	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/tmp	(revision 177)
@@ -1,2 +1,1 @@
-g
-C
+fast
Index: LMDZ.3.3/branches/rel-LF/tmp90
===================================================================
--- LMDZ.3.3/branches/rel-LF/tmp90	(revision 176)
+++ LMDZ.3.3/branches/rel-LF/tmp90	(revision 177)
@@ -11,4 +11,6 @@
       USE IOIPSL
 libf/dyn3d/abort_gcm.F
+      USE ioipsl
+libf/dyn3d/create_limit.F
       USE IOIPSL
 libf/dyn3d/dynredem.F
