Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/interfsurf.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/interfsurf.F	(revision 81)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/interfsurf.F	(revision 81)
@@ -0,0 +1,758 @@
+
+  MODULE interface_surf
+
+! Ce module regroupe toutes les routines gerant l'interface entre le modele 
+! atmospherique et les modeles de surface (sols continentaux, oceans, glaces)
+! Les routines sont les suivantes:
+!
+!   interfsurf_*: routines d'aiguillage vers les interfaces avec les 
+!                 differents modeles de surface
+!   interfsol\
+!             > routines d'interface proprement dite
+!   interfoce/
+!
+!   interfstart: routine d'initialisation et de lecture de l'etat initial
+!                "interface"
+!   interffin  : routine d'ecriture de l'etat de redemmarage de l'interface
+!
+! 
+! L. Fairhead, LMD, 02/2000
+
+  USE ioipsl
+  USE constantes
+
+  IMPLICIT none
+
+  PRIVATE
+  PUBLIC :: interfsurf
+
+  INTERFACE interfsurf
+    module procedure interfsurf_hq, interfsurf_vent
+  END INTERFACE
+
+  INTERFACE interfoce
+    module procedure interfoce_cpl, interfoce_slab, interfoce_lim
+  END INTERFACE
+
+
+! run_off      ruissellement total
+  real, allocatable, dimension(:),save    :: run_off
+
+
+  CONTAINS
+!
+!############################################################################
+!
+  SUBROUTINE interfsurf_hq(itime, dtime, jour, &
+     & klon, nisurf, knon, knindex, rlon, rlat, &
+     & debut, lafin, ok_veget, &
+     & zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, & 
+     & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
+     & precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
+     & ocean, &
+     & evap, fluxsens, fluxlat, &              
+     & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new)
+
+! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general 
+! (sols continentaux, oceans, glaces) pour les fluxs de chaleur et d'humidite.
+! En pratique l'interface se fait entre la couche limite du modele 
+! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
+!
+! 
+! 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)
+!   jour         jour dans l'annee en cours
+!   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
+!   rlat         latitudes
+!   debut        logical: 1er appel a la physique
+!   lafin        logical: dernier appel a la physique
+!   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
+!   zlflu        epaisseur de la premier 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
+!   hum_air      humidite 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 entrant a la surface
+!   swnet        flux solaire net
+!   swdown       flux solaire entrant a la surface
+!   ps           pression au sol
+!   ocean        type d'ocean utilise (force, slab, couple)
+!
+! 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
+
+  include 'indicesol.h'
+
+! Parametres d'entree
+  integer, intent(IN) :: itime
+  integer, intent(IN) :: klon
+  integer, intent(IN) :: dtime
+  integer, intent(IN) :: jour
+  integer, intent(IN) :: nisurf
+  integer, intent(IN) :: knon
+  integer, dimension(knon), intent(in) :: knindex
+  logical, intent(IN) :: debut, lafin, ok_veget
+  real, dimension(klon), intent(IN) :: rlon, rlat
+  real, dimension(knon), intent(IN) :: zlev, zlflu
+  real, dimension(knon), intent(IN) :: u1_lay, v1_lay
+  real, dimension(knon), intent(IN) :: temp_air, spechum
+  real, dimension(knon), intent(IN) :: hum_air, ccanopy
+  real, dimension(knon), intent(IN) :: tq_cdrag, petAcoef, peqAcoef
+  real, dimension(knon), intent(IN) :: petBcoef, peqBcoef
+  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
+  real, dimension(knon), intent(IN) :: lwdown, swnet, swdown, ps
+  character (len = 6)  :: ocean
+
+! Parametres de sortie
+  real, dimension(knon), intent(OUT):: evap, fluxsens, fluxlat
+  real, dimension(knon), intent(OUT):: tsol_rad, tsurf_new, alb_new
+  real, dimension(knon), intent(OUT):: emis_new, z0_new
+  real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
+
+! Local
+  character (len = 20) :: modname = 'interfsurf_hq'
+  character (len = 80) :: abort_message 
+  logical, save        :: first_call = .true.
+  integer              :: error
+
+
+!
+! On doit commencer par appeler les schemas de surfaces continentales
+! car l'ocean a besoin du ruissellement qui est y calcule
+!
+  if (first_call) then
+    if (nisurf /= is_ter) then 
+      write(*,*)' *** Warning ***'
+      write(*,*)' nisurf = ',nisurf,' /= is_ter = ',is_ter
+      write(*,*)'or on doit commencer par les surfaces continentales'
+      abort_message='voir ci-dessus'
+      call abort_gcm(modname,abort_message,1)
+    endif
+    if (ocean /= 'slab  ' .and. ocean /= 'force ' .and. ocean /= 'couple') then
+      write(*,*)' *** Warning ***'
+      write(*,*)'Option couplage pour l''ocean = ', ocean
+      abort_message='option pour l''ocean non valable'
+      call abort_gcm(modname,abort_message,1)
+    endif
+  endif
+  first_call = .false.
+!
+! Aiguillage vers les differents schemas de surface
+!
+  if (nisurf == is_ter) then
+!
+! Surface "terre" appel a l'interface avec les sols continentaux
+!
+! allocation du run-off
+    if (.not. allocated(run_off)) then
+      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
+      write(*,*)'Bizarre, le nombre de points continentaux'
+      write(*,*)'a change entre deux appels. Je continue ...'
+      deallocate(run_off, stat = error)
+      allocate(run_off(knon), stat = error)
+      if (error /= 0) then
+        abort_message='Pb allocation run_off'
+        call abort_gcm(modname,abort_message,1)
+      endif
+    endif
+!
+    call interfsol(itime, klon, dtime, nisurf, knon, &
+     & knindex, rlon, rlat, debut, lafin, &
+     & zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, 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)
+!
+  else if (nisurf == is_oce) then
+!
+! Surface "ocean" appel a l'interface avec l'ocean
+!
+    if (ocean == 'couple') then
+      call interfoce(nisurf, ocean)
+    else if (ocean == 'slab  ') then
+      call interfoce(nisurf)
+    else                              ! lecture conditions limites
+      call interfoce(itime, dtime, jour, & 
+     &  klon, nisurf, knon, knindex, &
+     &  debut, &
+     &  tsurf_new, alb_new, z0_new, pctsrf_new)
+    endif
+!
+  else if (nisurf == is_sic) then
+!
+! Surface "glace de mer" appel a l'interface avec l'ocean
+!
+    call interfoce(nisurf, ocean)
+!
+  else if (nisurf == is_lic) then
+!
+! Surface "glacier continentaux" appel a l'interface avec le sol
+!
+!    call interfsol(nisurf)
+  else
+    write(*,*)'Index surface = ',nisurf
+    abort_message = 'Index surface non valable'
+    call abort_gcm(modname,abort_message,1)
+  endif
+
+  END SUBROUTINE interfsurf_hq
+
+!
+!#########################################################################
+!
+  SUBROUTINE interfsurf_vent(nisurf, knon   &         
+  &                     )
+!
+! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general 
+! (sols continentaux, oceans, glaces) pour les tensions de vents.
+! En pratique l'interface se fait entre la couche limite du modele 
+! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
+!
+! 
+! L.Fairhead 02/2000
+!
+! input:
+!   nisurf       index de la surface a traiter (1 = sol continental)
+!   knon         nombre de points de la surface a traiter
+
+! Parametres d'entree
+  integer, intent(IN) :: nisurf
+  integer, intent(IN) :: knon
+
+
+  return
+  END SUBROUTINE interfsurf_vent
+!
+!#########################################################################
+!
+  SUBROUTINE interfsol(itime, klon, dtime, nisurf, knon, &
+     & knindex, rlon, rlat, debut, lafin, &
+     & zlev, zlflu, u1_lay, v1_lay, temp_air, spechum, 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)
+
+! 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
+!   debut        logical: 1er appel a la physique (lire les restart)
+!   lafin        logical: dernier appel a la physique (ecrire les restart)
+!   zlev         hauteur de la premiere couche
+!   zlflu        
+!   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
+!   hum_air      humidite 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 entrant a la surface
+!   swnet        flux solaire net
+!   swdown       flux solaire entrant a la surface
+!   ps           pression au sol
+!
+! 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
+  integer, intent(IN) :: dtime
+  integer, intent(IN) :: nisurf
+  integer, intent(IN) :: knon
+  integer, dimension(knon), intent(IN) :: knindex
+  logical, intent(IN) :: debut, lafin
+  real, dimension(klon), intent(IN) :: rlon, rlat
+  real, dimension(knon), intent(IN) :: zlev, zlflu
+  real, dimension(knon), intent(IN) :: u1_lay, v1_lay
+  real, dimension(knon), intent(IN) :: temp_air, spechum
+  real, dimension(knon), intent(IN) :: hum_air, ccanopy
+  real, dimension(knon), intent(IN) :: tq_cdrag, petAcoef, peqAcoef
+  real, dimension(knon), intent(IN) :: petBcoef, peqBcoef
+  real, dimension(knon), intent(IN) :: precip_rain, precip_snow
+  real, dimension(knon), intent(IN) :: lwdown, swnet, swdown, ps
+
+! Parametres de sortie
+  real, dimension(knon), intent(OUT):: evap, fluxsens, fluxlat
+  real, dimension(knon), intent(OUT):: tsol_rad, tsurf_new, alb_new
+  real, dimension(knon), intent(OUT):: emis_new, z0_new
+
+! Local
+!
+  integer              :: ii
+  integer              :: error
+  character (len = 20) :: modname = 'interfsol'
+  character (len = 80) :: abort_message
+! 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
+! resolution de la grille
+  real, allocatable, dimension (:,:), save :: resolution
+! Identifieurs des fichiers restart et histoire
+  integer, save          :: rest_id, hist_id 
+  integer, save          :: rest_id_stom, hist_id_stom
+
+
+! 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) 
+
+! 
+! Appel a la routine sols continentaux
+!
+
+  call sechiba_main(itime, klon, knon, knindex, dtime, &
+     & debut, lafin, coupling, control_in, &
+     & lalo, neighbours, resolution,&
+     & zlev, zlflu, 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
+!
+
+  END SUBROUTINE interfsol
+!
+!#########################################################################
+!
+  SUBROUTINE interfoce_cpl(nisurf, ocean)
+
+! Cette routine sert d'interface entre le modele atmospherique et un 
+! coupleur avec un modele d'ocean 'complet'
+!
+! L. Fairhead 02/2000
+!
+! input:
+!   nisurf       index de la surface a traiter (1 = sol continental)
+!   ocean        type d'ocean
+!
+! output:
+!
+
+! Parametres d'entree
+  integer, intent(IN) :: nisurf
+  character (len = 6)  :: ocean
+
+! Parametres de sortie
+
+! Variables locales
+
+
+! Initialisation
+! fichier restart et fichiers histoires
+
+! calcul des fluxs a passer
+
+  END SUBROUTINE interfoce_cpl
+!
+!#########################################################################
+!
+  SUBROUTINE interfoce_slab(nisurf)
+
+! Cette routine sert d'interface entre le modele atmospherique et un 
+! modele de 'slab' ocean
+!
+! L. Fairhead 02/2000
+!
+! input:
+!   nisurf       index de la surface a traiter (1 = sol continental)
+!
+! output:
+!
+
+! Parametres d'entree
+  integer, intent(IN) :: nisurf
+
+  END SUBROUTINE interfoce_slab
+!
+!#########################################################################
+!
+  SUBROUTINE interfoce_lim(itime, dtime, jour, &
+     & klon, nisurf, knon, knindex, &
+     & debut,  &
+     & lmt_sst, lmt_alb, lmt_rug, pctsrf_new)
+
+! Cette routine sert d'interface entre le modele atmospherique et un 
+!
+! L. Fairhead 02/2000
+!
+! input:
+!   itime        numero du pas de temps courant
+!   dtime        pas de temps de la physique (en s)
+!   jour         jour a lire dans l'annee
+!   nisurf       index de la surface a traiter (1 = sol continental)
+!   knon         nombre de points dans le domaine a traiter
+!   knindex      index des points de la surface a traiter
+!   klon         taille de la grille
+!   debut        logical: 1er appel a la physique (initialisation)
+!
+! output:
+!   lmt_sst      SST lues dans le fichier de CL
+!   lmt_alb      Albedo lu 
+!   lmt_rug      longueur de rugosité lue
+!   pctsrf_new   sous-maille fractionnelle
+!
+
+#include "indicesol.h"
+
+! Parametres d'entree
+  integer, intent(IN) :: itime
+  integer, intent(IN) :: dtime
+  integer, intent(IN) :: jour
+  integer, intent(IN) :: nisurf
+  integer, intent(IN) :: knon
+  integer, intent(IN) :: klon
+  integer, dimension(knon), intent(in) :: knindex
+  logical, intent(IN) :: debut
+
+! Parametres de sortie
+  real, intent(out), dimension(knon) :: lmt_sst
+  real, intent(out), dimension(knon) :: lmt_alb
+  real, intent(out), dimension(knon) :: lmt_rug
+  real, intent(out), dimension(klon,nbsrf) :: pctsrf_new
+
+! Variables locales
+  integer     :: ii
+  integer     :: lmt_pas     ! frequence de lecture des conditions limites 
+                             ! (en pas de physique)
+  logical,save :: deja_lu    ! pour indiquer que le jour a lire a deja
+                             ! lu pour une surface precedente
+  integer,save :: jour_lu 
+  integer      :: ierr
+  character (len = 20) :: modname = 'interfoce_lim'
+  character (len = 80) :: abort_message
+  character (len = 20) :: fich ='limit'
+  logical     :: newlmt = .false.
+  logical     :: check = .true.
+! Champs lus dans le fichier de CL
+  real, allocatable , save, dimension(:) :: sst_lu, alb_lu, rug_lu, nat_lu
+  real, allocatable , save, dimension(:,:) :: pct_tmp
+!
+! quelques variables pour netcdf
+!
+#include "netcdf.inc"
+  integer              :: nid, nvarid
+  integer, dimension(2) :: start, epais
+!
+! Fin déclaration
+!
+    
+  if (debut) then
+    lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
+    jour_lu = jour - 1
+    allocate(sst_lu(klon))
+    allocate(alb_lu(klon))
+    allocate(rug_lu(klon))
+    allocate(nat_lu(klon))
+    allocate(pct_tmp(klon,nbsrf))
+  endif
+
+  if ((jour - jour_lu) /= 0) deja_lu = .false.
+  
+  if (check) write(*,*)modname,' :: jour_lu, deja_lu', jour_lu, deja_lu 
+
+! Tester d'abord si c'est le moment de lire le fichier
+  if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu) then
+!
+! Ouverture du fichier
+!
+    ierr = NF_OPEN (fich, NF_NOWRITE,nid)
+    if (ierr.NE.NF_NOERR) then
+      abort_message = 'Pb d''ouverture du fichier de conditions aux limites'
+      call abort_gcm(modname,abort_message,1)
+    endif
+! 
+! La tranche de donnees a lire:
+!
+    start(1) = 1
+    start(2) = jour + 1
+    epais(1) = klon
+    epais(2) = 1
+!
+    if (newlmt) then
+!
+! Fraction "ocean" 
+!
+      ierr = NF_INQ_VARID(nid, 'FOCE', nvarid)
+      if (ierr /= NF_NOERR) then
+        abort_message = 'Le champ <FOCE> est absent'
+        call abort_gcm(modname,abort_message,1)
+      endif
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pctsrf_new(1,is_oce))
+#else
+      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pctsrf_new(1,is_oce))
+#endif
+      if (ierr /= NF_NOERR) then
+        abort_message = 'Lecture echouee pour <FOCE>'
+        call abort_gcm(modname,abort_message,1)
+      endif
+!
+! Fraction "glace de mer" 
+!
+      ierr = NF_INQ_VARID(nid, 'FSIC', nvarid)
+      if (ierr /= NF_NOERR) then
+        abort_message = 'Le champ <FSIC> est absent'
+        call abort_gcm(modname,abort_message,1)
+      endif
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pctsrf_new(1,is_sic))
+#else
+      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pctsrf_new(1,is_sic))
+#endif
+      if (ierr /= NF_NOERR) then
+        abort_message = 'Lecture echouee pour <FSIC>'
+        call abort_gcm(modname,abort_message,1)
+      endif
+!
+! Fraction "terre" 
+!
+      ierr = NF_INQ_VARID(nid, 'FTER', nvarid)
+      if (ierr /= NF_NOERR) then
+        abort_message = 'Le champ <FTER> est absent'
+        call abort_gcm(modname,abort_message,1)
+      endif
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pctsrf_new(1,is_ter))
+#else
+      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pctsrf_new(1,is_ter))
+#endif
+      if (ierr /= NF_NOERR) then
+        abort_message = 'Lecture echouee pour <FTER>'
+        call abort_gcm(modname,abort_message,1)
+      endif
+!
+! Fraction "glacier terre" 
+!
+      ierr = NF_INQ_VARID(nid, 'FLIC', nvarid)
+      if (ierr /= NF_NOERR) then
+        abort_message = 'Le champ <FLIC> est absent'
+        call abort_gcm(modname,abort_message,1)
+      endif
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pctsrf_new(1,is_lic))
+#else
+      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pctsrf_new(1,is_lic))
+#endif
+      if (ierr /= NF_NOERR) then
+        abort_message = 'Lecture echouee pour <FLIC>'
+        call abort_gcm(modname,abort_message,1)
+      endif
+!
+    else  ! on en est toujours a rnatur
+! 
+      ierr = NF_INQ_VARID(nid, 'NAT', nvarid)
+      if (ierr /= NF_NOERR) then
+        abort_message = 'Le champ <NAT> est absent'
+        call abort_gcm(modname,abort_message,1)
+      endif
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, nat_lu)
+#else
+      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, nat_lu)
+#endif
+      if (ierr /= NF_NOERR) then
+        abort_message = 'Lecture echouee pour <NAT>'
+        call abort_gcm(modname,abort_message,1)
+      endif
+!
+! Remplissage des fractions de surface
+! nat = 0, 1, 2, 3 pour ocean, terre, glacier, seaice
+! 
+      pct_tmp = 0.0
+      do ii = 1, klon
+        pct_tmp(ii,nint(nat_lu(ii)) + 1) = 1.
+      enddo
+
+!
+!  On se retrouve avec ocean en 1 et terre en 2 alors qu'on veut le contraire
+!
+      pctsrf_new = pct_tmp
+      pctsrf_new (:,2)= pct_tmp (:,1)
+      pctsrf_new (:,1)= pct_tmp (:,2)
+      pct_tmp = pctsrf_new 
+    endif ! fin test sur newlmt
+!
+! Lecture SST
+!
+    ierr = NF_INQ_VARID(nid, 'SST', nvarid)
+    if (ierr /= NF_NOERR) then
+      abort_message = 'Le champ <SST> est absent'
+      call abort_gcm(modname,abort_message,1)
+    endif
+#ifdef NC_DOUBLE
+    ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, sst_lu)
+#else
+    ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, sst_lu)
+#endif
+    if (ierr /= NF_NOERR) then
+      abort_message = 'Lecture echouee pour <SST>'
+      call abort_gcm(modname,abort_message,1)
+    endif    
+!
+! Lecture Albedo
+!
+    ierr = NF_INQ_VARID(nid, 'ALB', nvarid)
+    if (ierr /= NF_NOERR) then
+      abort_message = 'Le champ <ALB> est absent'
+      call abort_gcm(modname,abort_message,1)
+    endif
+#ifdef NC_DOUBLE
+    ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, alb_lu)
+#else
+    ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, alb_lu)
+#endif
+    if (ierr /= NF_NOERR) then
+      abort_message = 'Lecture echouee pour <ALB>'
+      call abort_gcm(modname,abort_message,1)
+    endif
+!
+! Lecture rugosité
+!
+    ierr = NF_INQ_VARID(nid, 'RUG', nvarid)
+    if (ierr /= NF_NOERR) then
+      abort_message = 'Le champ <RUG> est absent'
+      call abort_gcm(modname,abort_message,1)
+    endif
+#ifdef NC_DOUBLE
+    ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, rug_lu)
+#else
+    ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, rug_lu)
+#endif
+    if (ierr /= NF_NOERR) then
+      abort_message = 'Lecture echouee pour <RUG>'
+      call abort_gcm(modname,abort_message,1)
+    endif
+
+!
+! Fin de lecture
+!
+    ierr = NF_CLOSE(nid)
+    deja_lu = .true.
+    jour_lu = jour
+  endif
+!
+! Recopie des variables dans les champs de sortie
+!
+  do ii = 1, knon
+    lmt_sst(ii) = sst_lu(knindex(ii))
+    lmt_alb(ii) = alb_lu(knindex(ii))
+    lmt_rug(ii) = rug_lu(knindex(ii))
+  enddo
+  pctsrf_new = pct_tmp
+
+  END SUBROUTINE interfoce_lim
+
+!
+!#########################################################################
+!
+  END MODULE interface_surf
