Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 88)
+++ 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
 
