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