Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 440)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 441)
@@ -57,10 +57,10 @@
       & klon, iim, jjm, nisurf, knon, knindex, pctsrf, &
       & rlon, rlat, cufi, cvfi,&
-      & debut, lafin, ok_veget, soil_model, nsoilmx, tsoil,&
+      & debut, lafin, ok_veget, soil_model, nsoilmx, tsoil, qsol,&
       & zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, & 
       & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
       & precip_rain, precip_snow, sollw, sollwdown, swnet, swdown, &
       & fder, taux, tauy, rugos, rugoro, &
-      & albedo, snow, qsol, &
+      & albedo, snow, qsurf, &
       & tsurf, p1lay, ps, radsol, &
       & ocean, npas, nexca, zmasq, &
@@ -171,9 +171,10 @@
   character (len = 6)  :: ocean
   integer              :: npas, nexca ! nombre et pas de temps couplage
-  real, dimension(klon), intent(INOUT) :: evap, snow, qsol
+  real, dimension(klon), intent(INOUT) :: evap, snow, qsurf
 !! PB ajout pour soil
   logical          :: soil_model
   integer          :: nsoilmx
   REAL, DIMENSION(klon, nsoilmx) :: tsoil
+  REAL, dimension(klon), intent(INOUT) :: qsol
   REAL, dimension(klon)          :: soilcap
   REAL, dimension(klon)          :: soilflux
@@ -193,5 +194,5 @@
   integer, save        :: error
   integer              :: ii, index
-  logical,save              :: check = .false.
+  logical,save              :: check = .true.
   real, dimension(klon):: cal, beta, dif_grnd, capsol
 !!$PB  real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5.
@@ -200,5 +201,5 @@
   real, dimension(klon):: alb_ice
   real, dimension(klon):: tsurf_temp
-  real, dimension(klon):: qsol_new
+  real, dimension(klon):: qsurf_new
 !!  real, allocatable, dimension(:), save :: alb_neig_grid
   real, dimension(klon):: alb_neig, alb_eau
@@ -309,5 +310,5 @@
      & alb_new, z0_new)
 !  
-! calcul snow et qsol, hydrol adapté
+! calcul snow et qsurf, hydrol adapté
 !
        CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
@@ -323,5 +324,5 @@
        CALL calcul_fluxs( klon, knon, nisurf, dtime, &
      &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &
-     &   precip_rain, precip_snow, snow, qsol,  &
+     &   precip_rain, precip_snow, snow, qsurf,  &
      &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
      &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
@@ -358,5 +359,5 @@
      &  evap, fluxsens, fluxlat, &              
      &  tsol_rad, tsurf_new, alb_new, alblw, &
-     &  emis_new, z0_new, dflux_l, dflux_s, qsol_new)
+     &  emis_new, z0_new, dflux_l, dflux_s, qsurf_new)
 
 !  
@@ -366,5 +367,5 @@
 !
 ! mise a jour de l'humidite saturante calculee par ORCHIDEE
-      qsol(1:knon) = qsol_new(1:knon)
+      qsurf(1:knon) = qsurf_new(1:knon)
 
     endif    
@@ -416,5 +417,5 @@
     call calcul_fluxs( klon, knon, nisurf, dtime, &
      &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
-     &   precip_rain, precip_snow, snow, qsol,  &
+     &   precip_rain, precip_snow, snow, qsurf,  &
      &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
      &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
@@ -511,5 +512,5 @@
     CALL calcul_fluxs( klon, knon, nisurf, dtime, &
          &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
-         &   precip_rain, precip_snow, snow, qsol,  &
+         &   precip_rain, precip_snow, snow, qsurf,  &
          &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
          &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
@@ -581,5 +582,5 @@
     call calcul_fluxs( klon, knon, nisurf, dtime, &
      &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &
-     &   precip_rain, precip_snow, snow, qsol,  &
+     &   precip_rain, precip_snow, snow, qsurf,  &
      &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
      &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
@@ -651,5 +652,5 @@
      & knindex, rlon, rlat, cufi, cvfi, iim, jjm, pctsrf, &
      & debut, lafin, ok_veget, &
-     & zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, & 
+     & plev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, & 
      & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
      & precip_rain, precip_snow, lwdown, swnet, swdown, &
@@ -680,5 +681,5 @@
 !   ok_veget     logical: appel ou non au schema de surface continental
 !                     (si false calcul simplifie des fluxs sur les continents)
-!   zlev         hauteur de la premiere couche       
+!   plev         hauteur de la premiere couche (Pa)      
 !   u1_lay       vitesse u 1ere couche
 !   v1_lay       vitesse v 1ere couche
@@ -715,5 +716,5 @@
 !   emis_new     emissivite
 !   z0_new       surface roughness
-!   qsurf        saturated air moisture at surface
+!   qsurf        air moisture at surface
 
 ! Parametres d'entree
@@ -730,5 +731,5 @@
   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) :: plev
   real, dimension(klon), intent(IN) :: u1_lay, v1_lay
   real, dimension(klon), intent(IN) :: temp_air, spechum
@@ -739,4 +740,7 @@
   real, dimension(klon), intent(IN) :: precip_rain, precip_snow
   real, dimension(klon), intent(IN) :: lwdown, swnet, swdown, ps
+!IM cf. JP +++
+  real, dimension(klon) :: swdown_vrai
+!IM cf. JP ---
   real, dimension(klon), intent(IN) :: tsurf, p1lay
   real, dimension(klon), intent(IN) :: radsol
@@ -759,4 +763,8 @@
 ! drapeaux controlant les appels dans SECHIBA
 !  type(control_type), save   :: control_in
+! Preserved albedo
+!IM cf. JP +++
+  real, allocatable, dimension(:), save :: albedo_keep, zlev
+!IM cf. JP ---
 ! coordonnees geographiques
   real, allocatable, dimension(:,:), save :: lalo
@@ -794,4 +802,5 @@
 
 #include "temps.inc"
+#include "YOMCST.inc"
 
   if (check) write(*,*)'Entree ', modname
@@ -803,4 +812,8 @@
   if (debut) then
 
+  IF ( .NOT. allocated(albedo_keep)) THEN
+     ALLOCATE(albedo_keep(klon))
+     ALLOCATE(zlev(klon))
+  ENDIF
 ! Pb de correspondances de grilles
    allocate(ig(klon))
@@ -990,4 +1003,9 @@
   cdrag(1:knon) = tq_cdrag(1:knon)
 
+!IM cf. JP +++
+! zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/287.05*temp_air(1:knon))*9.80665)
+  zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/RD*temp_air(1:knon))*RG)
+!IM cf. JP ---
+
    where(cdrag > 0.01) 
      cdrag = 0.01
@@ -1008,5 +1026,14 @@
      & tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
      & lon_scat, lat_scat)
+
+!IM cf. JP +++
+    albedo_keep(:) = (albedo_out(:,1)+albedo_out(:,2))/2.
+!IM cf. JP ---
+
   endif
+
+!IM cf. JP +++
+  swdown_vrai(:) = swnet(:)/(1. - albedo_keep(:))
+!IM cf. JP ---
 
   call intersurf_main (itime+itau_phy, iim, jjm+1, knon, ktindex, dtime, &
@@ -1015,8 +1042,14 @@
      & zlev,  u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
      & cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
-     & precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
+!IM cf. JP +++
+     & precip_rain, precip_snow, lwdown, swnet, swdown_vrai, ps, &
+!IM cf. JP ---
      & evap, fluxsens, fluxlat, coastalflow, riverflow, &
      & tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
      & lon_scat, lat_scat)
+
+!IM cf. JP +++
+    albedo_keep(:) = (albedo_out(:,1)+albedo_out(:,2))/2.
+!IM cf. JP ---
 
     bidule=0.
@@ -1927,5 +1960,5 @@
   character (len = 20),save :: fich ='limit.nc'
   logical,save     :: newlmt = .false.
-  logical,save     :: check = .FALSE.
+  logical,save     :: check = .true.
 ! Champs lus dans le fichier de CL
   real, allocatable , save, dimension(:) :: alb_lu, rug_lu
@@ -1951,5 +1984,5 @@
   if (check) write(*,*)modname,':: jour_lu_sur, deja_lu_sur', jour_lu_sur, deja_lu_sur 
   if (check) write(*,*)modname,':: itime, lmt_pas', itime, lmt_pas
-  call flush(6)
+  if (check) call flush(6)
 
 ! Tester d'abord si c'est le moment de lire le fichier
@@ -1959,4 +1992,6 @@
 !
     fich = trim(fich)
+    IF (check) WRITE(*,*)modname,' ouverture fichier ',fich
+    if (check) CALL flush(6)
     ierr = NF_OPEN (fich, NF_NOWRITE,nid)
     if (ierr.NE.NF_NOERR) then
@@ -2033,5 +2068,5 @@
   SUBROUTINE calcul_fluxs( klon, knon, nisurf, dtime, &
      & tsurf, p1lay, cal, beta, coef1lay, ps, &
-     & precip_rain, precip_snow, snow, qsol, &
+     & precip_rain, precip_snow, snow, qsurf, &
      & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
      & petAcoef, peqAcoef, petBcoef, peqBcoef, &
@@ -2055,5 +2090,4 @@
 !   precip_snow  precipitations solides
 !   snow         champs hauteur de neige
-!   qsol         humidite du sol
 !   runoff       runoff en cas de trop plein
 !   petAcoef     coeff. A de la resolution de la CL pour t
@@ -2066,4 +2100,5 @@
 ! output:
 !   tsurf_new    temperature au sol
+!   qsurf        humidite de l'air au dessus du sol
 !   fluxsens     flux de chaleur sensible
 !   fluxlat      flux de chaleur latente
@@ -2086,5 +2121,5 @@
   real, dimension(klon), intent(IN) :: radsol, dif_grnd
   real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay
-  real, dimension(klon), intent(INOUT) :: snow, qsol
+  real, dimension(klon), intent(INOUT) :: snow, qsurf
 
 ! Parametres sorties
@@ -2102,9 +2137,10 @@
   real                  :: bilan_f, fq_fonte
   REAL                  :: subli, fsno
+  REAL                  :: qsat_new, q1_new
   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, save         :: check = .FALSE.
+  logical, save         :: check = .true.
   character (len = 20)  :: modname = 'calcul_fluxs'
   logical, save         :: fonte_neige = .false.
@@ -2114,4 +2150,10 @@
 
   if (check) write(*,*)'Entree ', modname,' surface = ',nisurf
+
+  IF (check) THEN
+      WRITE(*,*)' radsol (min, max)' &
+         &     , MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
+      CALL flush(6)
+  ENDIF
 
   if (size(coastalflow) /= knon .AND. nisurf == is_ter) then
@@ -2135,5 +2177,5 @@
 !!$      where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
 !!$    endif 
-    IF (nisurf /= is_ter) qsol = max_eau_sol
+!!$    IF (nisurf /= is_ter) qsol = max_eau_sol
 
 
@@ -2233,5 +2275,8 @@
     dflux_s(i) = zx_nh(i)
     dflux_l(i) = (zx_sl(i) * zx_nq(i))
-
+! Nouvelle valeure de l'humidite au dessus du sol
+    qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
+    q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
+    qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new
 !
 ! en cas de fonte de neige
@@ -2412,5 +2457,5 @@
 !   precip_snow  precipitations solides
 !   snow         champs hauteur de neige
-!   qsol         humidite du sol
+!   qsol         hauteur d'eau contenu dans le sol
 !   runoff       runoff en cas de trop plein
 !   petAcoef     coeff. A de la resolution de la CL pour t
@@ -2459,4 +2504,5 @@
   real                  :: bilan_f, fq_fonte
   REAL                  :: subli, fsno
+  real, dimension(klon) :: bil_eau_s
   real, parameter :: t_grnd = 271.35, t_coup = 273.15
 !! PB temporaire en attendant mieux pour le modele de neige
@@ -2473,4 +2519,5 @@
 
 ! Initialisations
+  bil_eau_s(:) = 0.
   DO i = 1, knon
     zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
@@ -2534,5 +2581,5 @@
   WHERE (precip_snow > 0.) snow = snow + (precip_snow * dtime)
   WHERE (evap > 0 ) snow = MAX(0.0, snow - (evap * dtime))
-  qsol = qsol + (precip_rain - evap) * dtime
+  bil_eau_s = bil_eau_s + (precip_rain - evap) * dtime
 !
 ! Y'a-t-il fonte de neige?
@@ -2544,5 +2591,5 @@
       fq_fonte = MIN( MAX((tsurf_new(i)-RTT )/chasno,0.0),snow(i))
       snow(i) = max(0., snow(i) - fq_fonte)
-      qsol(i) = qsol(i) + fq_fonte 
+      bil_eau_s(i) = bil_eau_s(i) + fq_fonte 
       tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno  
       IF (nisurf == is_sic .OR. nisurf == is_lic ) tsurf_new(i) = RTT -1.8
@@ -2564,7 +2611,12 @@
 !!$      fq_fonte = bilan_f / zx_sl(i)
     endif
-    IF (nisurf == is_ter)  &
-       &  run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0)
-    qsol(i) = MIN(qsol(i), max_eau_sol) 
+
+    IF (nisurf == is_ter) then
+      qsol(i) = qsol(i) + bil_eau_s(i)
+      run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0)
+      qsol(i) = MIN(qsol(i), max_eau_sol) 
+    else
+      run_off(i) = run_off(i) + MAX(bil_eau_s(i), 0.0)
+    endif
   enddo
 
