Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 234)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 235)
@@ -282,5 +282,5 @@
       y_flux_u = 0.0
       y_flux_v = 0.0
-      ytsoil = 0.0
+      ytsoil = 999999.
 
       DO nsrf = 1, nbsrf
@@ -330,5 +330,5 @@
 
       DO 99999 nsrf = 1, nbsrf
-c$$$   PB   totalflu = radsol
+       totalflu = radsol
 
 c chercher les indices:
@@ -384,5 +384,5 @@
         yu1(j) = u1lay(i)
         yv1(j) = v1lay(i)
-c$$$    PB    yrads(j) = totalflu(i)
+c$$$        yrads(j) = totalflu(i)
         yrads(j) = (1 - albe(i,nsrf)) 
      $      /(1 - pctsrf(i,is_ter) * albe(i,is_ter) 
@@ -493,4 +493,9 @@
       evap(:,nsrf) = - flux_q(:,1,nsrf)
 c
+      albe(:, nsrf) = 0.
+      snow(:, nsrf) = 0.
+      qsol(:, nsrf) = 0.
+      rugos(:, nsrf) = 0.
+      fluxlat(:,nsrf) = 0.
       DO j = 1, knon
          i = ni(j)
@@ -502,5 +507,8 @@
          fluxlat(i,nsrf) = yfluxlat(j)
 c$$$ pb         rugmer(i) = yrugm(j)
-         IF (nsrf .EQ. is_oce) rugmer(i) = yrugm(j)
+         IF (nsrf .EQ. is_oce) then 
+	       rugmer(i) = yrugm(j)
+               rugos(i,nsrf) = yrugm(i)
+          endif	
          cdragh(i) = cdragh(i) + ycoefh(j,1)
          cdragm(i) = cdragm(i) + ycoefm(j,1)
@@ -511,4 +519,5 @@
       END DO
 c$$$ PB ajout pour soil
+      ftsoil(:,:,nsrf) = 0.
       DO k = 1, nsoilmx
         DO j = 1, knon
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 234)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 235)
@@ -3,14 +3,8 @@
 !
 
-!module conf_phys
-!
-!  use IOIPSL
-!  implicit none
-!
-!  public conf_phys
-!
-!contains
+  subroutine conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan)
 
-  subroutine conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan)
+   use IOIPSL
+   implicit none
 
 !
@@ -21,6 +15,4 @@
 !
 
-  use IOIPSL     
-  implicit none
 !
 ! ocean:      type d'ocean (force, slab, couple)
@@ -104,3 +96,2 @@
   end subroutine conf_phys
 
-!end module conf_phys
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 234)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90	(revision 235)
@@ -159,7 +159,7 @@
   real, dimension(klon), intent(IN) :: ps, albedo
   real, dimension(klon), intent(IN) :: tsurf, p1lay
-  real, dimension(klon), intent(INOUT) :: radsol
+  REAL, DIMENSION(klon), INTENT(INOUT) :: radsol,fder
   real, dimension(klon), intent(IN) :: zmasq
-  real, dimension(klon), intent(IN) :: fder, taux, tauy, rugos, rugoro
+  real, dimension(klon), intent(IN) :: taux, tauy, rugos, rugoro
   character (len = 6)  :: ocean
   integer              :: npas, nexca ! nombre et pas de temps couplage
@@ -196,5 +196,4 @@
   real, DIMENSION(klon):: zfra
   logical              :: cumul = .false.
-  logical,save         :: scatter = .false.
 
   if (check) write(*,*) 'Entree ', modname
@@ -235,7 +234,12 @@
 ! Initialisations diverses
 !
-  cal=0.; beta=1.; dif_grnd=0.; capsol=0.
-  alb_new = 0.; z0_new = 0.; alb_neig = 0.0
-
+!!$  cal=0.; beta=1.; dif_grnd=0.; capsol=0.
+!!$  alb_new = 0.; z0_new = 0.; alb_neig = 0.0
+!!$! PB
+!!$  tsurf_new = 0.
+
+  cal = 999999. ; beta = 999999. ; dif_grnd = 999999. ; capsol = 999999.
+  alb_new = 999999. ; z0_new = 999999. ; alb_neig = 999999.
+  tsurf_new = 999999.
 ! Aiguillage vers les differents schemas de surface
 
@@ -290,7 +294,6 @@
 !
 !!$ PB ATTENTION changement ordre des appels
-!!$   CALL albsno(klon,agesno,alb_neig_grid)  
- 
- 
+    CALL albsno(klon,agesno,alb_neig_grid)  
+
     if (.not. ok_veget) then
 !
@@ -410,4 +413,6 @@
      &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
      &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+    
+    fder = fder + dflux_s + dflux_l
 
 !
@@ -472,32 +477,49 @@
 !    else if (ocean == 'slab  ') then
 !      call interfoce(nisurf)
-    else                              ! lecture conditions limites
-      call interfoce(itime, dtime, jour, & 
-     &  klon, nisurf, knon, knindex, &
-     &  debut, &
-     &  tsurf_new, pctsrf_new)
-
-      tsurf_temp = tsurf
-           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, &
-     &   tsurf_temp, 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)
-
-    if (ocean /= 'couple') then
-      call fonte_neige( klon, knon, nisurf, dtime, &
-     &   tsurf_temp, 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)
-    endif
+      ELSE
+!                              ! lecture conditions limites
+          CALL interfoce(itime, dtime, jour, & 
+             &  klon, nisurf, knon, knindex, &
+             &  debut, &
+             &  tsurf_new, pctsrf_new)
+
+          CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
+      
+          IF (soil_model) THEN 
+              CALL soil(dtime, nisurf, knon,snow, tsurf, tsoil,soilcap, soilflux)
+              cal(1:knon) = RCPD / soilcap(1:knon)
+              radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
+              dif_grnd = 0.
+!!$           WRITE(*,*) 'radsol'
+!!$           WRITE(*,*) radsol(1 : knon)
+!!$           WRITE(*,*) 'soilflux'
+!!$           WRITE(*,*) soilflux(1 : knon)
+          ELSE 
+!      if (check) write(*,*)'Sortie calbeta'
+!      if (check) write(*,*)'RCPD = ',RCPD,' capsol = '
+!      if (check) write(*,*)capsol
+              dif_grnd = 1.0 / tau_gl
+              cal = RCPD * calice
+              WHERE (snow > 0.0) cal = RCPD * calsno 
+          ENDIF
+          tsurf_temp = tsurf
+          beta = 1.0
+      ENDIF
+
+      CALL calcul_fluxs( klon, knon, nisurf, dtime, &
+         &   tsurf_temp, 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)
+
+      IF (ocean /= 'couple') THEN
+          CALL fonte_neige( klon, knon, nisurf, dtime, &
+             &   tsurf_temp, 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)
+      ENDIF
 !
 ! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
@@ -539,7 +561,11 @@
 !    call interfsol(nisurf)
     IF (soil_model) THEN 
-        CALL soil(dtime, nisurf, snow, tsurf, tsoil,soilcap, soilflux)
-        cal = RCPD / soilcap
-        radsol = radsol + soilflux
+        CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil,soilcap, soilflux)
+        cal(1:knon) = RCPD / soilcap(1:knon)
+        radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
+!!$           WRITE(*,*) 'radsol'
+!!$           WRITE(*,'(16f17.4)') radsol(1 : knon)
+!!$           WRITE(*,*) 'soilflux'
+!!$           WRITE(*,'(16f17.4)')soilflux(1:knon)
     ELSE 
         cal = RCPD * calice
@@ -962,282 +988,4 @@
 
   END SUBROUTINE interfsol
-!
-!#########################################################################
-!
-  SUBROUTINE interfsol_scat(itime, klon, dtime, date0, nisurf, knon, &
-     & knindex, rlon, rlat, cufi, cvfi, iim, jjm, pctsrf, &
-     & debut, lafin, ok_veget, &
-     & 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, &
-     & 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 
-! modele de sol continental. Appel a sechiba
-!
-! L. Fairhead 02/2000
-!
-! input:
-!   itime        numero du pas de temps
-!   klon         nombre total de points de grille
-!   dtime        pas de temps de la physique (en s)
-!   nisurf       index de la surface a traiter (1 = sol continental)
-!   knon         nombre de points de la surface a traiter
-!   knindex      index des points de la surface a traiter
-!   rlon         longitudes de la grille entiere
-!   rlat         latitudes de la grille entiere
-!   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)
-!   ok_veget     logical: appel ou non au schema de surface continental
-!                     (si false calcul simplifie des fluxs sur les continents)
-!   zlev         hauteur de la premiere couche       
-!   u1_lay       vitesse u 1ere couche
-!   v1_lay       vitesse v 1ere couche
-!   temp_air     temperature de l'air 1ere couche
-!   spechum      humidite specifique 1ere couche
-!   epot_air     temp pot de l'air
-!   ccanopy      concentration CO2 canopee
-!   tq_cdrag     cdrag
-!   petAcoef     coeff. A de la resolution de la CL pour t
-!   peqAcoef     coeff. A de la resolution de la CL pour q
-!   petBcoef     coeff. B de la resolution de la CL pour t
-!   peqBcoef     coeff. B de la resolution de la CL pour q
-!   precip_rain  precipitation liquide
-!   precip_snow  precipitation solide
-!   lwdown       flux IR descendant a la surface
-!   swnet        flux solaire net
-!   swdown       flux solaire entrant a la surface
-!   tsurf        temperature de surface
-!   p1lay        pression 1er niveau (milieu de couche)
-!   ps           pression au sol
-!   radsol       rayonnement net aus sol (LW + SW)
-!   
-!
-! input/output
-!   run_off      ruissellement total
-!
-! output:
-!   evap         evaporation totale
-!   fluxsens     flux de chaleur sensible
-!   fluxlat      flux de chaleur latente
-!   tsol_rad     
-!   tsurf_new    temperature au sol
-!   alb_new      albedo
-!   emis_new     emissivite
-!   z0_new       surface roughness
-
-
-! Parametres d'entree
-  integer, intent(IN) :: itime
-  integer, intent(IN) :: klon
-  real, intent(IN)    :: dtime
-  real, intent(IN)    :: date0
-  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) :: 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
-  real, dimension(klon), intent(IN) :: lwdown, swnet, swdown, ps
-  real, dimension(klon), intent(IN) :: tsurf, p1lay
-  real, dimension(klon), intent(IN) :: radsol
-! Parametres de sortie
-  real, dimension(klon), intent(OUT):: evap, fluxsens, fluxlat
-  real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new
-  real, dimension(klon), intent(OUT):: emis_new, z0_new
-  real, dimension(klon), intent(OUT):: dflux_s, dflux_l
-
-! Local
-!
-  integer              :: ii, ij, jj, igrid, ireal, i, index
-  integer              :: error
-  character (len = 20) :: modname = 'interfsol_scat'
-  character (len = 80) :: abort_message
-  logical,save              :: check = .TRUE.
-  real, dimension(klon) :: cal, beta, dif_grnd, capsol
-! type de couplage dans sechiba
-!  character (len=10)   :: coupling = 'implicit' 
-! drapeaux controlant les appels dans SECHIBA
-!  type(control_type), save   :: control_in
-! coordonnees geographiques
-  real, allocatable, dimension(:,:), save :: lalo
-! 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, save           :: lrestart_read = .true. , lrestart_write = .false.
-
-  real, dimension(klon):: qsurf
-  real, dimension(klon):: snow, qsol
-  real, dimension(knon,2) :: albedo_out
-! Pb de nomenclature
-  real, dimension(klon)      :: petA_orc, petB_orc, peqA_orc, peqB_orc
-! champs a passer a ORCHIDEE
-  real, dimension(:,:), allocatable, save :: lon_sc, lat_sc, contfrac_sc
-  real, dimension(iim,jjm+1) :: zlev_sc
-  real, dimension(iim,jjm+1) :: u1_lay_sc, v1_lay_sc, spechum_sc, temp_air_sc
-  real, dimension(iim,jjm+1) :: epot_air_sc, ccanopy_sc, tq_cdrag_sc
-  real, dimension(iim,jjm+1) :: petA_or_sc, peqA_or_sc, petB_or_sc, peqB_or_sc
-  real, dimension(iim,jjm+1) :: precip_rain_sc, precip_snow_sc, lwdown_sc
-  real, dimension(iim,jjm+1) :: swnet_sc, swdown_sc, p1lay_sc, evap_sc
-  real, dimension(iim,jjm+1) :: fluxsens_sc, fluxlat_sc, coastalflow_sc
-  real, dimension(iim,jjm+1) :: riverflow_sc, tsol_rad_sc, tsurf_new_sc
-  real, dimension(iim,jjm+1) :: qsurf_sc, emis_new_sc, z0_new_sc
-  real, dimension(iim,jjm+1,2) :: albedo_out_sc
-
-  if (check) write(*,*)'Entree ', modname
-  if (check) write(*,*)'ok_veget = ',ok_veget
-
-! initialisation
-  if (debut) then
-    if ((.not. allocated(lon_sc))) then
-      allocate(lon_sc(iim,jjm+1), stat = error)
-      if (error /= 0) then
-        abort_message='Pb allocation lon_sc'
-        call abort_gcm(modname,abort_message,1)
-      endif      
-    endif
-    if ((.not. allocated(lat_sc))) then
-      allocate(lat_sc(iim,jjm+1), stat = error)
-      if (error /= 0) then
-        abort_message='Pb allocation lat_sc'
-        call abort_gcm(modname,abort_message,1)
-      endif      
-    endif
-    index = 1
-    do jj = 2, jjm
-      do ij = 1, iim
-        index = index + 1
-        lon_sc(ij,jj) = rlon(index)
-        lat_sc(ij,jj) = rlat(index)
-      enddo
-    enddo
-    lon_sc(:,1) = lon_sc(:,2)
-    lat_sc(:,1) = rlat(1)
-    lon_sc(:,jjm+1) = lon_sc(:,jjm)
-    lat_sc(:,jjm+1) = rlat(klon)
-    if (( .not. allocated(contfrac_sc))) then
-      allocate(contfrac_sc(iim, jjm+1), stat = error)
-      if (error /= 0) then
-        abort_message='Pb allocation contfrac'
-        call abort_gcm(modname,abort_message,1)
-      endif      
-    endif
-    contfrac_sc = 0.
-    call gath2cpl(pctsrf(:,is_ter), contfrac_sc, klon, knon,iim,jjm, knindex)
-  endif                          ! (fin debut) 
-
-! 
-! Appel a la routine sols continentaux
-!
-! petit pb de nomenclature
-!
-  if (lafin) lrestart_write = .true.
-  if (check) write(*,*)'lafin ',lafin,lrestart_write
-
-  petA_orc = petBcoef * dtime
-  petB_orc = petAcoef
-  peqA_orc = peqBcoef * dtime
-  peqB_orc = peqAcoef
-
-!
-! Passage sur la grille 2D avant envoi a ORCHIDEE
-!
-  call gath2cpl(zlev, zlev_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(u1_lay, u1_lay_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(v1_lay, v1_lay_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(spechum, spechum_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(temp_air, temp_air_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(epot_air, epot_air_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(ccanopy, ccanopy_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(tq_cdrag, tq_cdrag_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(petA_orc, petA_or_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(peqA_orc, peqA_or_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(petB_orc, petB_or_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(peqB_orc, peqB_or_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(precip_rain, precip_rain_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(precip_snow, precip_snow_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(lwdown, lwdown_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(swnet, swnet_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(swdown, swdown_sc, klon, knon,iim,jjm, knindex)
-  call gath2cpl(p1lay, p1lay_sc, klon, knon,iim,jjm, knindex)
-!
-! Init Orchidee
-!
-  if (debut) then
-    if (check) write(*,*) 'debut orchidee itime - 1', itime-1,date0
-    call intersurf_main (itime-1, iim, jjm+1 , knon, knindex, dtime, &
-     & lrestart_read, lrestart_write, lon_sc, lat_sc, &
-     & contfrac_sc, date0, zlev_sc,  u1_lay_sc, v1_lay_sc, &
-     & spechum_sc, temp_air_sc, epot_air_sc, ccanopy_sc, &
-     & tq_cdrag_sc, petA_or_sc, peqA_or_sc, petB_or_sc, peqB_or_sc, &
-     & precip_rain_sc, precip_snow_sc, lwdown_sc, swnet_sc, &
-     & swdown_sc, p1lay_sc, & 
-     & evap_sc, fluxsens_sc, fluxlat_sc, &
-     & coastalflow_sc, riverflow_sc, tsol_rad_sc, tsurf_new_sc, &
-     & qsurf_sc, albedo_out_sc, emis_new_sc, z0_new_sc)
-  endif
-
-  call intersurf_main (itime, iim, jjm+1 , knon, knindex, dtime, &
-     & lrestart_read, lrestart_write, lon_sc, lat_sc, &
-     & contfrac_sc, date0, zlev_sc,  u1_lay_sc, v1_lay_sc, &
-     & spechum_sc, temp_air_sc, epot_air_sc, ccanopy_sc, &
-     & tq_cdrag_sc, petA_or_sc, peqA_or_sc, petB_or_sc, peqB_or_sc, &
-     & precip_rain_sc, precip_snow_sc, lwdown_sc, swnet_sc, &
-     & swdown_sc, p1lay_sc, & 
-     & evap_sc, fluxsens_sc, fluxlat_sc, &
-     & coastalflow_sc, riverflow_sc, tsol_rad_sc, tsurf_new_sc, &
-     & qsurf_sc, albedo_out_sc, emis_new_sc, z0_new_sc)
-! 
-! sorties mises sur la grille physique
-!
-
-  call cpl2gath(evap_sc, evap, klon, knon,iim,jjm, knindex)
-  call cpl2gath(fluxsens_sc, fluxsens, klon, knon,iim,jjm, knindex)
-  call cpl2gath(fluxlat_sc, fluxlat, klon, knon,iim,jjm, knindex)
-  call cpl2gath(coastalflow_sc, coastalflow, klon, knon,iim,jjm, knindex)
-  call cpl2gath(riverflow_sc, riverflow, klon, knon,iim,jjm, knindex)
-  call cpl2gath(tsol_rad_sc, tsol_rad, klon, knon,iim,jjm, knindex)
-  call cpl2gath(tsurf_new_sc, tsurf_new, klon, knon,iim,jjm, knindex)
-  call cpl2gath(qsurf_sc, qsurf, klon, knon,iim,jjm, knindex)
-  call cpl2gath(albedo_out_sc(:,:,1), alb_new, klon, knon,iim,jjm, knindex)
-  call cpl2gath(emis_new_sc, emis_new, klon, knon,iim,jjm, knindex)
-  call cpl2gath(z0_new_sc, z0_new, klon, knon,iim,jjm, knindex)
-
-! LF essai sensible
-  fluxsens = -1. * fluxsens
-  fluxlat = -1. * fluxlat
-
-  if (debut) lrestart_read = .false.
-
-  END SUBROUTINE interfsol_scat
 !
 !#########################################################################
@@ -2021,4 +1769,5 @@
 ! Recopie des variables dans les champs de sortie
 !
+  lmt_sst = 999999999.
   do ii = 1, knon
     lmt_sst(ii) = sst_lu(knindex(ii))
@@ -2176,6 +1925,8 @@
 ! Recopie des variables dans les champs de sortie
 !
-  lmt_alb(:) = 0.0
-  lmt_rug(:) = 0.0
+!!$  lmt_alb(:) = 0.0
+!!$  lmt_rug(:) = 0.0
+  lmt_alb(:) = 999999.
+  lmt_rug(:) = 999999.
   DO ii = 1, knon
     lmt_alb(ii) = alb_lu(knindex(ii))
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 234)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F	(revision 235)
@@ -447,4 +447,5 @@
 c
       CHARACTER*2 str2
+      CHARACTER*2 iqn
 c
       REAL qcheck
@@ -549,4 +550,5 @@
 c
       REAL tr_seri(klon,klev,nbtr)
+      REAL d_tr(klon,klev,nbtr)
 
       REAL zx_rh(klon,klev)
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/phytrac.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/phytrac.F	(revision 234)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/phytrac.F	(revision 235)
@@ -242,7 +242,7 @@
             zx_lon(i,jjm+1) = xlon(i+1)
          ENDDO
-         DO ll=1,klev
-            znivsig(ll)=float(ll)
-         ENDDO
+c         DO ll=1,klev
+c            znivsig(ll)=float(ll)
+c         ENDDO
          CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
          CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/raddim.h
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/raddim.h	(revision 234)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/raddim.h	(revision 235)
@@ -1,16 +1,2 @@
       INTEGER kdlon, kflev
-c
-ccc      PARAMETER (kdlon=klon,kflev=klev)
-c
-c resolution 72 45:
-      PARAMETER (kdlon=317,kflev=klev)
-c resolution 64 32:
-ccc      PARAMETER (kdlon=331,kflev=klev)
-c resolution 96 49:
-ccc      PARAMETER (kdlon=461,kflev=klev)
-c resolution 144 73:
-ccc      PARAMETER (kdlon=610,kflev=klev)
-c resolution 96 72:
-c      PARAMETER (kdlon=487,kflev=klev)
-c resolution 128 64:
-ccc       PARAMETER (kdlon=4033,kflev=klev)
+      PARAMETER (kdlon=klon,kflev=klev)
Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/soil.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/soil.F	(revision 234)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/soil.F	(revision 235)
@@ -1,3 +1,3 @@
-      SUBROUTINE soil(ptimestep, indice, snow, ptsrf, ptsoil,
+      SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil,
      s          pcapcal, pfluxgrd)
       IMPLICIT NONE
@@ -47,4 +47,5 @@
 
 #include "dimensions.h"
+#include "YOMCST.h"
 #include "dimphy.h"
 #include "dimsoil.h"
@@ -56,5 +57,5 @@
 
       REAL ptimestep
-      INTEGER indice
+      INTEGER indice, knon
       REAL ptsrf(klon),ptsoil(klon,nsoilmx),snow(klon)
       REAL pcapcal(klon),pfluxgrd(klon)
@@ -65,5 +66,6 @@
 
       INTEGER ig,jk
-      REAL zdz2(nsoilmx),z1(klon)
+c$$$      REAL zdz2(nsoilmx),z1(klon)
+      REAL zdz2(nsoilmx),z1(klon,nbsrf)
       REAL min_period,dalph_soil
       REAL ztherm_i(klon)
@@ -72,13 +74,15 @@
 c   ----------------------
       REAL dz1(nsoilmx),dz2(nsoilmx)
-      REAL zc(klon,nsoilmx),zd(klon,nsoilmx)
+c$$$          REAL zc(klon,nsoilmx),zd(klon,nsoilmx)
+      REAL zc(klon,nsoilmx,nbsrf),zd(klon,nsoilmx,nbsrf)
       REAL lambda
       SAVE dz1,dz2,zc,zd,lambda
-      LOGICAL firstcall
-      SAVE firstcall
+      LOGICAL firstcall, firstsurf(nbsrf)
+      SAVE firstcall, firstsurf
       REAL isol,isno,iice
       SAVE isol,isno,iice
 
       DATA firstcall/.true./
+      DATA firstsurf/.TRUE.,.TRUE.,.TRUE.,.TRUE./
 
       DATA isol,isno,iice/2000.,2000.,2000./
@@ -90,5 +94,5 @@
       REAL fz,rk,fz1,rk1,rk2
       fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
-
+      pfluxgrd(:) = 0.
 c   calcul de l'inertie thermique a partir de la variable rnat.
 c   on initialise a iice meme au-dessus d'un point de mer au cas 
@@ -97,20 +101,20 @@
 c
       IF (indice.EQ.is_sic) THEN
-         DO ig = 1, klon
+         DO ig = 1, knon
             ztherm_i(ig)   = iice
             IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno
          ENDDO
       ELSE IF (indice.EQ.is_lic) THEN
-         DO ig = 1, klon
+         DO ig = 1, knon
             ztherm_i(ig)   = iice
             IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno
          ENDDO
       ELSE IF (indice.EQ.is_ter) THEN
-         DO ig = 1, klon
+         DO ig = 1, knon
             ztherm_i(ig)   = isol
             IF (snow(ig).GT.0.0) ztherm_i(ig)   = isno
          ENDDO
       ELSE IF (indice.EQ.is_oce) THEN
-         DO ig = 1, klon
+         DO ig = 1, knon
             ztherm_i(ig)   = iice
          ENDDO
@@ -121,5 +125,6 @@
 
 
-      IF (firstcall) THEN
+c$$$      IF (firstcall) THEN
+      IF (firstsurf(indice)) THEN 
 
 c-----------------------------------------------------------------------
@@ -162,6 +167,7 @@
      .               fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
          ENDDO
-
-         firstcall =.false.
+C PB
+         firstsurf(indice) = .FALSE. 
+c$$$         firstcall =.false.
 
 c   Initialisations:
@@ -175,13 +181,14 @@
 
 c    surface temperature
-         DO ig=1,klon
-            ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/
-     s      (lambda*(1.-zd(ig,1))+1.)
+         DO ig=1,knon
+            ptsoil(ig,1)=(lambda*zc(ig,1,indice)+ptsrf(ig))/
+     s      (lambda*(1.-zd(ig,1,indice))+1.)
          ENDDO
 
 c   other temperatures
          DO jk=1,nsoilmx-1
-            DO ig=1,klon
-               ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk)
+            DO ig=1,knon
+               ptsoil(ig,jk+1)=zc(ig,jk,indice)+zd(ig,jk,indice)
+     $            *ptsoil(ig,jk)
             ENDDO
          ENDDO
@@ -192,20 +199,30 @@
 c   ---------------------------------------------------------------
 
+c$$$  PB ajout pour cas glace de mer
+      IF (indice .EQ. is_sic) THEN
+          DO ig = 1 , knon
+            ptsoil(ig,nsoilmx) = RTT - 1.8
+          END DO 
+      ENDIF 
+
       DO jk=1,nsoilmx
          zdz2(jk)=dz2(jk)/ptimestep
       ENDDO
 
-      DO ig=1,klon
-         z1(ig)=zdz2(nsoilmx)+dz1(nsoilmx-1)
-         zc(ig,nsoilmx-1)=zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig)
-         zd(ig,nsoilmx-1)=dz1(nsoilmx-1)/z1(ig)
+      DO ig=1,knon
+         z1(ig,indice)=zdz2(nsoilmx)+dz1(nsoilmx-1)
+         zc(ig,nsoilmx-1,indice)=
+     $       zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig,indice)
+         zd(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1(ig,indice)
       ENDDO
 
       DO jk=nsoilmx-1,2,-1
-         DO ig=1,klon
-            z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk)))
-            zc(ig,jk-1)=
-     s      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig)
-            zd(ig,jk-1)=dz1(jk-1)*z1(ig)
+         DO ig=1,knon
+            z1(ig,indice)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)
+     $         *(1.-zd(ig,jk,indice)))
+            zc(ig,jk-1,indice)=
+     s      (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice))
+     $          *z1(ig,indice)
+            zd(ig,jk-1,indice)=dz1(jk-1)*z1(ig,indice)
          ENDDO
       ENDDO
@@ -216,13 +233,15 @@
 c   ---------------------------------
 
-      DO ig=1,klon
+      DO ig=1,knon
          pfluxgrd(ig)=ztherm_i(ig)*dz1(1)*
-     s   (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1))
+     s   (zc(ig,1,indice)+(zd(ig,1,indice)-1.)*ptsoil(ig,1))
          pcapcal(ig)=ztherm_i(ig)*
-     s   (dz2(1)+ptimestep*(1.-zd(ig,1))*dz1(1))
-         z1(ig)=lambda*(1.-zd(ig,1))+1.
-         pcapcal(ig)=pcapcal(ig)/z1(ig)
-         pfluxgrd(ig)=pfluxgrd(ig)
-     s   +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig))
+     s   (dz2(1)+ptimestep*(1.-zd(ig,1,indice))*dz1(1))
+         z1(ig,indice)=lambda*(1.-zd(ig,1,indice))+1.
+         pcapcal(ig)=pcapcal(ig)/z1(ig,indice)
+         pfluxgrd(ig) = pfluxgrd(ig)
+     s   + pcapcal(ig) * (ptsoil(ig,1) * z1(ig,indice)
+     $       - lambda * zc(ig,1,indice)
+     $       - ptsrf(ig))
      s   /ptimestep
       ENDDO
