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,interfsurf_hq 

  INTERFACE interfsurf
    module procedure interfsurf_hq, interfsurf_vent
  END INTERFACE

  INTERFACE interfoce
    module procedure interfoce_cpl, interfoce_slab, interfoce_lim
  END INTERFACE

#include "YOMCST.inc"
#include "indicesol.inc"


! run_off      ruissellement total
  real, allocatable, dimension(:),save    :: run_off



  CONTAINS
!
!############################################################################
!
  SUBROUTINE interfsurf_hq(itime, dtime, jour, rmu0, &
      & klon, iim, jjm, nisurf, knon, knindex, pctsrf, rlon, rlat, &
      & debut, lafin, ok_veget, &
      & zlev,  u1_lay, v1_lay, temp_air, spechum, hum_air, ccanopy, & 
      & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
      & precip_rain, precip_snow, lwdown, swnet, swdown, &
      & fder, taux, tauy, &
      & albedo, snow, qsol, &
      & tsurf, p1lay, ps, radsol, &
      & ocean, npas, nexca, zmasq, &
      & evap, fluxsens, fluxlat, dflux_l, dflux_s, &              
      & tsol_rad, tsurf_new, alb_new, emis_new, z0_new, pctsrf_new, agesno)


! 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
!   iim, jjm     nbres de pts de grille
!   dtime        pas de temps de la physique (en s)
!   jour         jour dans l'annee en cours,
!   rmu0         cosinus de l'angle solaire zenithal
!   nexca        pas de temps couplage
!   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
!   pctsrf       tableau des pourcentages de surface de chaque maille
!   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
!   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
!   albedo       albedo de la surface
!   tsurf        temperature de surface
!   p1lay        pression 1er niveau (milieu de couche)
!   ps           pression au sol
!   radsol       rayonnement net aus sol (LW + SW)
!   ocean        type d'ocean utilise (force, slab, couple)
!   fder         derivee des flux (pour le couplage)
!   taux, tauy   tension de vents
!   zmasq        masque terre/ocean
!
! 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
!   pctsrf_new   nouvelle repartition des surfaces


! Parametres d'entree
  integer, intent(IN) :: itime
  integer, intent(IN) :: iim, jjm
  integer, intent(IN) :: klon
  real, intent(IN) :: dtime
  integer, intent(IN) :: jour
  real, intent(IN)    :: rmu0(klon)
  integer, intent(IN) :: nisurf
  integer, intent(IN) :: knon
  integer, dimension(klon), intent(in) :: knindex
  real, dimension(klon,nbsrf), intent(IN) :: pctsrf
  logical, intent(IN) :: debut, lafin, ok_veget
  real, dimension(klon), intent(IN) :: rlon, rlat
  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) :: hum_air, ccanopy
  real, dimension(klon), intent(IN) :: tq_cdrag, 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, albedo
  real, dimension(klon), intent(IN) :: tsurf, p1lay
  real, dimension(klon), intent(IN) :: radsol
  real, dimension(klon), intent(IN) :: zmasq
  real, dimension(klon), intent(IN) :: fder, taux, tauy
  character (len = 6)  :: ocean
  integer              :: npas, nexca ! nombre et pas de temps couplage
  real, dimension(klon), intent(INOUT) :: evap, snow, qsol

! Parametres de sortie
  real, dimension(klon), intent(OUT):: 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_l, dflux_s
  real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
  real, dimension(klon), intent(INOUT):: agesno

! Local
  character (len = 20) :: modname = 'interfsurf_hq'
  character (len = 80) :: abort_message 
  logical, save        :: first_call = .true.
  INTEGER              :: error, ii
  logical              :: check = .true.
  real, dimension(klon):: cal, beta, dif_grnd, capsol
  real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=1./86400.*5.
  real, parameter      :: calsno=1./(2.3867e+06*.15)
  real, dimension(klon):: alb_ice
  real, dimension(klon):: tsurf_temp
  real, dimension(klon):: alb_neig_grid, alb_eau
  real, dimension(klon):: alb_neig
  real, DIMENSION(klon):: zfra

  if (check) write(*,*) 'Entree ', modname
!
! 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 .and. klon > 1) 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
    if ( is_oce > is_sic ) then
      write(*,*)' *** Warning ***'
      write(*,*)' Pour des raisons de sequencement dans le code'
      write(*,*)' l''ocean doit etre traite avant la banquise'
      write(*,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
      abort_message='voir ci-dessus'
      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
!
! Calcul age de la neige
!

  CALL albsno(klon,agesno,alb_neig_grid)  
  
  
 
    if (.not. ok_veget) then
!
! calcul snow et qsol, hydrol adapté
!
      call calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
      cal = RCPD * capsol
      call calcul_fluxs( klon, knon, nisurf, dtime, &
     &   tsurf, 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)

!
! calcul albedo: lecture albedo fichier CL puis ajout albedo neige 
! 
       call interfsur_lim(itime, dtime, jour, &
     & klon, nisurf, knon, knindex, debut,  &
     & alb_new, z0_new)
!
!
       DO ii = 1, knon
         alb_neig(ii) = alb_neig_grid(knindex(ii))
       enddo
       zfra = MAX(0.0,MIN(1.0,snow/(snow+10.0)))
       alb_new = alb_neig*zfra + alb_new*(1.0-zfra)
    else
!
!  appel a sechiba
!
      call interfsol(itime, klon, dtime, nisurf, knon, &
     &  knindex, rlon, rlat, &
     &  debut, lafin, ok_veget, &
     &  zlev,  u1_lay, v1_lay, temp_air, spechum, hum_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)
    endif    
!
! Remplissage des pourcentages de surface
!
    pctsrf_new(:,nisurf) = pctsrf(:,nisurf)

  else if (nisurf == is_oce) then

    if (check) write(*,*)'ocean, nisurf = ',nisurf 


!
! Surface "ocean" appel a l'interface avec l'ocean
!
    if (ocean == 'couple') then
!     nexca = 0
      if (nexca == 0) then
        abort_message='nexca = 0 dans interfoce_cpl'
        call abort_gcm(modname,abort_message,1)
      endif

      call interfoce(itime, dtime, &
      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
      & ocean, npas, nexca, debut, lafin, &
      & swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
      & fder, albedo, taux, tauy, zmasq, &
      & tsurf_new, alb_new, alb_ice, pctsrf_new)

!    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)

    endif

    tsurf_temp = tsurf_new
    cal = 0.
    beta = 1.
    dif_grnd = 0.

    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)
!
! calcul albedo
!

    if ( minval(rmu0) == maxval(rmu0) .and. minval(rmu0) == -999.999 ) then
      CALL alboc(FLOAT(jour),rlat,alb_eau)
    else  ! cycle diurne
      CALL alboc_cd(rmu0,alb_eau)
    endif
    DO ii =1, knon
      alb_new(ii) = alb_eau(knindex(ii))
    enddo
!
  else if (nisurf == is_sic) then

    if (check) write(*,*)'sea ice, nisurf = ',nisurf 

!
! Surface "glace de mer" appel a l'interface avec l'ocean
!
!
    if (ocean == 'couple') then

      call interfoce(itime, dtime, &
      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
      & ocean, npas, nexca, debut, lafin, &
      & swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
      & fder, albedo, taux, tauy, zmasq, &
      & tsurf_new, alb_new, alb_ice, pctsrf_new)

      tsurf_temp = tsurf_new
      dif_grnd = 0.
      beta = 1.0

!    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
    endif

    cal = calice
    where (snow > 0.0) cal = calsno

    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)

!
! calcul albedo
!
       zfra = MAX(0.0,MIN(1.0,snow/(snow+10.0)))
       DO ii = 1, knon
         alb_neig = alb_neig_grid(knindex(ii))
       enddo
       alb_new = alb_neig*zfra + 0.6 * (1.0-zfra)

  else if (nisurf == is_lic) then

    if (check) write(*,*)'glacier, nisurf = ',nisurf 

!
! Surface "glacier continentaux" appel a l'interface avec le sol
!
!    call interfsol(nisurf)

    cal = calice
    where (snow > 0.0) cal = calsno
    beta = 1.0
    dif_grnd = 0.0

    call calcul_fluxs( klon, knon, nisurf, dtime, &
     &   tsurf, 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)

!
! calcul albedo
!
       zfra = MAX(0.0,MIN(1.0,snow/(snow+10.0)))
       DO ii =1, knon
         alb_neig = alb_neig_grid(knindex(ii))
       enddo
       alb_new = alb_neig*zfra + 0.6 * (1.0-zfra)
!
! Remplissage des pourcentages de surface
!
    pctsrf_new(:,nisurf) = pctsrf(:,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, ok_veget, &
     & zlev,  u1_lay, v1_lay, temp_air, spechum, hum_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)

! 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)
!   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
!   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
!   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
  integer, intent(IN) :: nisurf
  integer, intent(IN) :: knon
  integer, dimension(klon), intent(IN) :: knindex
  logical, intent(IN) :: debut, lafin, ok_veget
  real, dimension(klon), intent(IN) :: rlon, rlat
  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) :: hum_air, ccanopy
  real, dimension(klon), intent(IN) :: tq_cdrag, 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
  integer              :: error
  character (len = 20) :: modname = 'interfsol'
  character (len = 80) :: abort_message
  logical              :: 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
! 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

  real, dimension(klon):: snow, qsol

  if (check) write(*,*)'Entree ', modname
  if (check) write(*,*)'ok_veget = ',ok_veget

! 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,  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(itime, dtime, &
      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
      & ocean, npas, nexca, debut, lafin, &
      & swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
      & fder, albsol, taux, tauy, zmasq, &
      & tsurf_new, alb_new, alb_ice, pctsrf_new)

! Cette routine sert d'interface entre le modele atmospherique et un 
! coupleur avec un modele d'ocean 'complet' derriere
!
! Le modele de glace qu'il est prevu d'utiliser etant couple directement a 
! l'ocean presentement, on va passer deux fois dans cette routine par pas de 
! temps physique, une fois avec les points oceans et l'autre avec les points
! glace. A chaque pas de temps de couplage, la lecture des champs provenant
! du coupleur se fera "dans" l'ocean et l'ecriture des champs a envoyer
! au coupleur "dans" la glace. Il faut donc des tableaux de travail "tampons"
! dimensionnes sur toute la grille qui remplissent les champs sur les
! domaines ocean/glace quand il le faut. Il est aussi necessaire que l'index
! ocean soit traiter avant l'index glace (sinon tout intervertir)
!
!
! L. Fairhead 02/2000
!
! input:
!   itime        numero du pas de temps
!   iim, jjm     nbres de pts de grille
!   dtime        pas de tempsde la physique
!   klon         nombre total de points de grille
!   nisurf       index de la surface a traiter (1 = sol continental)
!   pctsrf       tableau des fractions de surface de chaque maille
!   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
!   ocean        type d'ocean
!   nexca        frequence de couplage
!   swdown       flux solaire entrant a la surface
!   lwdown       flux IR entrant a la surface
!   precip_rain  precipitation liquide
!   precip_snow  precipitation solide
!   evap         evaporation
!   tsurf        temperature de surface
!   fder         derivee dF/dT
!   albsol       albedo du sol (coherent avec swdown)
!   taux         tension de vent en x
!   tauy         tension de vent en y
!   nexca        frequence de couplage
!   zmasq        masque terre/ocean
!
!
! output:
!   tsurf_new    temperature au sol
!   alb_new      albedo
!   pctsrf_new   nouvelle repartition des surfaces
!   alb_ice      albedo de la glace
!


! Parametres d'entree
  integer, intent(IN) :: itime
  integer, intent(IN) :: iim, jjm
  real, intent(IN) :: dtime
  integer, intent(IN) :: klon
  integer, intent(IN) :: nisurf
  integer, intent(IN) :: knon
  real, dimension(klon,nbsrf), intent(IN) :: pctsrf
  integer, dimension(klon), intent(in) :: knindex
  logical, intent(IN) :: debut, lafin
  real, dimension(klon), intent(IN) :: rlon, rlat
  character (len = 6)  :: ocean
  real, dimension(klon), intent(IN) :: lwdown, swdown
  real, dimension(klon), intent(IN) :: precip_rain, precip_snow
  real, dimension(klon), intent(IN) :: tsurf, fder, albsol, taux, tauy
  INTEGER              :: nexca, npas, kstep
  real, dimension(klon), intent(IN) :: zmasq

  real, dimension(klon), intent(INOUT) :: evap

! Parametres de sortie
  real, dimension(klon), intent(OUT):: tsurf_new, alb_new, alb_ice
  real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new

! Variables locales
  integer                    :: j, error, sum_error, ig, cpl_index,i
  character (len = 20) :: modname = 'interfoce_cpl'
  character (len = 80) :: abort_message
  logical              :: check = .true.
! variables pour moyenner les variables de couplage
  real, allocatable, dimension(:,:),save :: cpl_sols, cpl_nsol, cpl_rain
  real, allocatable, dimension(:,:),save :: cpl_snow, cpl_evap, cpl_tsol
  real, allocatable, dimension(:,:),save :: cpl_fder, cpl_albe, cpl_taux
  real, allocatable, dimension(:,:),save :: cpl_tauy, cpl_rriv, cpl_rcoa
! variables tampons avant le passage au coupleur
  real, allocatable, dimension(:,:,:),save :: tmp_sols, tmp_nsol, tmp_rain
  real, allocatable, dimension(:,:,:),save :: tmp_snow, tmp_evap, tmp_tsol
  real, allocatable, dimension(:,:,:),save :: tmp_fder, tmp_albe, tmp_taux
  real, allocatable, dimension(:,:,:),save :: tmp_tauy, tmp_rriv, tmp_rcoa
! variables a passer au coupleur
  real, dimension(iim, jjm+1) :: wri_sol_ice, wri_sol_sea, wri_nsol_ice 
  real, dimension(iim, jjm+1) :: wri_nsol_sea, wri_fder_ice, wri_evap_ice
  real, dimension(iim, jjm+1) :: wri_evap_sea
  real, dimension(iim, jjm+1) :: wri_rain, wri_snow, wri_taux
  real, dimension(iim, jjm+1) :: wri_tauy, wri_rriv, wri_rcoa
! variables relues par le coupleur
! read_sic = fraction de glace
! read_sit = temperature de glace
  real, allocatable, dimension(:,:),save :: read_sst, read_sic, read_sit
  real, allocatable, dimension(:,:),save :: read_alb_sic
! variable tampon
  real, dimension(klon)       :: tamp
  real, dimension(klon)       :: tamp_sic
! sauvegarde des fractions de surface d'un pas de temps a l'autre apres 
! l'avoir lu
  real, allocatable,dimension(:,:),save :: pctsrf_sav
  real, dimension(iim, jjm+1, 2) :: tamp_srf
  integer, allocatable, dimension(:), save :: tamp_ind
  real, allocatable, dimension(:,:),save :: tamp_zmasq
  real, dimension(iim, jjm+1) :: deno
  integer                     :: idtime
  integer, allocatable, dimension(:,:) :: isst
  integer, allocatable,save,dimension(:) :: unity
! 
  logical, save    :: first_appel = .true.
  logical          :: print

!
! Initialisation
!
  if (check) write(*,*)'Entree ',modname,'nisurf = ',nisurf
 
  if (first_appel) then
    error = 0
    allocate(unity(klon), stat = error)
    if ( error  /=0) then
      abort_message='Pb allocation variable unity'
      call abort_gcm(modname,abort_message,1)
    endif
    allocate(pctsrf_sav(klon,nbsrf), stat = error)
    if ( error  /=0) then
      abort_message='Pb allocation variable pctsrf_sav'
      call abort_gcm(modname,abort_message,1)
    endif
    pctsrf_sav = 0.

    do ig = 1, klon
      unity(ig) = ig
    enddo
    sum_error = 0
    allocate(cpl_sols(klon,2), stat = error); sum_error = sum_error + error
    allocate(cpl_nsol(klon,2), stat = error); sum_error = sum_error + error
    allocate(cpl_rain(klon,2), stat = error); sum_error = sum_error + error
    allocate(cpl_snow(klon,2), stat = error); sum_error = sum_error + error
    allocate(cpl_evap(klon,2), stat = error); sum_error = sum_error + error
    allocate(cpl_tsol(klon,2), stat = error); sum_error = sum_error + error
    allocate(cpl_fder(klon,2), stat = error); sum_error = sum_error + error
    allocate(cpl_albe(klon,2), stat = error); sum_error = sum_error + error
    allocate(cpl_taux(klon,2), stat = error); sum_error = sum_error + error
    allocate(cpl_tauy(klon,2), stat = error); sum_error = sum_error + error
    allocate(cpl_rcoa(klon,2), stat = error); sum_error = sum_error + error
    allocate(cpl_rriv(klon,2), stat = error); sum_error = sum_error + error
    allocate(read_sst(iim, jjm+1), stat = error); sum_error = sum_error + error
    allocate(read_sic(iim, jjm+1), stat = error); sum_error = sum_error + error
    allocate(read_sit(iim, jjm+1), stat = error); sum_error = sum_error + error
    allocate(read_alb_sic(iim, jjm+1), stat = error); sum_error = sum_error + error

    if (sum_error /= 0) then
      abort_message='Pb allocation variables couplees'
      call abort_gcm(modname,abort_message,1)
    endif
    cpl_sols = 0.; cpl_nsol = 0.; cpl_rain = 0.; cpl_snow = 0.
    cpl_evap = 0.; cpl_tsol = 0.; cpl_fder = 0.; cpl_albe = 0.
    cpl_taux = 0.; cpl_tauy = 0.; cpl_rriv = 0.; cpl_rcoa = 0.

    sum_error = 0
    allocate(tamp_ind(klon), stat = error); sum_error = sum_error + error
    allocate(tamp_zmasq(iim, jjm+1), stat = error); sum_error = sum_error + error    
    do ig = 1, klon
      tamp_ind(ig) = ig
    enddo
    call gath2cpl(zmasq, tamp_zmasq, klon, klon, iim, jjm, tamp_ind)
!
! initialisation couplage
!
    idtime = int(dtime)
    call inicma(npas , nexca, idtime,(jjm+1)*iim)

    first_appel = .false.
  endif ! fin if (first_appel)

! fichier restart et fichiers histoires

! calcul des fluxs a passer

  cpl_index = 1
  if (nisurf == is_sic) cpl_index = 2
  do ig = 1, knon
    cpl_sols(ig,cpl_index) = cpl_sols(ig,cpl_index) &
     &                               + swdown(ig)      / FLOAT(nexca)
    cpl_nsol(ig,cpl_index) = cpl_nsol(ig,cpl_index) &
     &                               + lwdown(ig)      / FLOAT(nexca)
    cpl_rain(ig,cpl_index) = cpl_rain(ig,cpl_index) &
     &                               + precip_rain(ig) / FLOAT(nexca)
    cpl_snow(ig,cpl_index) = cpl_snow(ig,cpl_index) &
     &                               + precip_snow(ig) / FLOAT(nexca)
    cpl_evap(ig,cpl_index) = cpl_evap(ig,cpl_index) &
     &                               + evap(ig)        / FLOAT(nexca)
    cpl_tsol(ig,cpl_index) = cpl_tsol(ig,cpl_index) &
     &                               + tsurf(ig)       / FLOAT(nexca)
    cpl_fder(ig,cpl_index) = cpl_fder(ig,cpl_index) &
     &                               + fder(ig)        / FLOAT(nexca)
    cpl_albe(ig,cpl_index) = cpl_albe(ig,cpl_index) &
     &                               + albsol(ig)      / FLOAT(nexca)
    cpl_taux(ig,cpl_index) = cpl_taux(ig,cpl_index) &
     &                               + taux(ig)        / FLOAT(nexca)
    cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) &
     &                               + tauy(ig)        / FLOAT(nexca)
    cpl_rriv(ig,cpl_index) = cpl_rriv(ig,cpl_index) &
     &                               + 0.     / FLOAT(nexca)/dtime
    cpl_rcoa(ig,cpl_index) = cpl_rcoa(ig,cpl_index) &
     &                               + 0.     / FLOAT(nexca)/dtime
  enddo

  if (mod(itime, nexca) == 1) then
!
! Passage des champs au coupleur
!
! Si le domaine considere est l'ocean, on lit les champs venant du coupleur
!
    if (nisurf == is_oce) then
      if (check) write(*,*)'rentree fromcpl, itime-1 = ',itime-1
      call fromcpl(itime-1,(jjm+1)*iim,                                  &
     &        read_sst, read_sic, read_sit, read_alb_sic)
      do j = 1, jjm + 1
        do ig = 1, iim
          if (abs(1. - read_sic(ig,j)) < 0.00001) then
            read_sst(ig,j) = RTT - 1.8
            read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
            read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
          else if (abs(read_sic(ig,j)) < 0.00001) then
            read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
            read_sit(ig,j) = read_sst(ig,j)
            read_alb_sic(ig,j) =  0.6
          else
            read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
            read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
            read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
          endif
        enddo
      enddo
!
! transformer read_sic en pctsrf_sav
!
    call cpl2gath(read_sic, tamp_sic , klon, klon,iim,jjm, unity)
    do ig = 1, klon
      IF (pctsrf(ig,is_oce) > epsfra .OR.            &
     &             pctsrf(ig,is_sic) > epsfra) THEN
          pctsrf_sav(ig,is_sic) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
     &                               * tamp_sic(ig)
          pctsrf_sav(ig,is_oce) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
     &                        - pctsrf_sav(ig,is_sic)
        endif
      enddo
      if (check) then
        write(46,*)'pct_srf_sav_ice = '
        write(46,'(72f8.3)')pctsrf_sav(:,is_sic)
      endif
      if (minval(pctsrf_new(:,is_oce)) < 0.) then
        write(*,*)'Pb fraction ocean inferieure a 0'
        write(*,*)'au point ',minloc(pctsrf_new(:,is_oce))
        write(*,*)'valeur = ',minval(pctsrf_new(:,is_oce)) 
        abort_message = 'voir ci-dessus'
        call abort_gcm(modname,abort_message,1)
      endif
      if (minval(pctsrf_new(:,is_sic)) < 0.) then
        write(*,*)'Pb fraction glace inferieure a 0'
        write(*,*)'au point ',minloc(pctsrf_new(:,is_sic))
        write(*,*)'valeur = ',minval(pctsrf_new(:,is_sic)) 
        abort_message = 'voir ci-dessus'
        call abort_gcm(modname,abort_message,1)
      endif
    endif 
  endif                         ! fin mod(itime, nexca) == 1

  if (mod(itime, nexca) == 0) then
!
! allocation memoire
    if (nisurf == is_oce) then
      sum_error = 0
      allocate(tmp_sols(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      allocate(tmp_nsol(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      allocate(tmp_rain(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      allocate(tmp_snow(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      allocate(tmp_evap(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      allocate(tmp_tsol(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      allocate(tmp_fder(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      allocate(tmp_albe(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      allocate(tmp_taux(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      allocate(tmp_tauy(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      allocate(tmp_rriv(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      allocate(tmp_rcoa(iim,jjm+1,2), stat=error); sum_error = sum_error + error
      if (sum_error /= 0) then
        abort_message='Pb allocation variables couplees pour l''ecriture'
        call abort_gcm(modname,abort_message,1)
      endif
    endif

!
! Mise sur la bonne grille des champs a passer au coupleur
!
    cpl_index = 1
    if (nisurf == is_sic) cpl_index = 2
    call gath2cpl(cpl_sols(1,cpl_index), tmp_sols(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)
    call gath2cpl(cpl_nsol(1,cpl_index), tmp_nsol(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)
    call gath2cpl(cpl_rain(1,cpl_index), tmp_rain(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)
    call gath2cpl(cpl_snow(1,cpl_index), tmp_snow(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)
    call gath2cpl(cpl_evap(1,cpl_index), tmp_evap(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)
    call gath2cpl(cpl_tsol(1,cpl_index), tmp_tsol(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)
    call gath2cpl(cpl_fder(1,cpl_index), tmp_fder(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)
    call gath2cpl(cpl_albe(1,cpl_index), tmp_albe(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)
    call gath2cpl(cpl_taux(1,cpl_index), tmp_taux(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)
    call gath2cpl(cpl_tauy(1,cpl_index), tmp_tauy(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)
    call gath2cpl(cpl_rriv(1,cpl_index), tmp_rriv(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)
    call gath2cpl(cpl_rcoa(1,cpl_index), tmp_rcoa(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)

!
! Si le domaine considere est la banquise, on envoie les champs au coupleur
!
    if (nisurf == is_sic) then
      wri_rain = 0.; wri_snow = 0.; wri_rcoa = 0.; wri_rriv = 0. 
      wri_taux = 0.; wri_tauy = 0.
      call gath2cpl(pctsrf(1,is_oce), tamp_srf(1,1,1), klon, klon, iim, jjm, tamp_ind)
      call gath2cpl(pctsrf(1,is_sic), tamp_srf(1,1,2), klon, klon, iim, jjm, tamp_ind)

      wri_sol_ice = tmp_sols(:,:,2)
      wri_sol_sea = tmp_sols(:,:,1)
      wri_nsol_ice = tmp_nsol(:,:,2)
      wri_nsol_sea = tmp_nsol(:,:,1)
      wri_fder_ice = tmp_fder(:,:,2)
      wri_evap_ice = tmp_evap(:,:,2)
      wri_evap_sea = tmp_evap(:,:,1)
      where (tamp_zmasq /= 1.)
        deno =  tamp_srf(:,:,1) + tamp_srf(:,:,2)
        wri_rain = tmp_rain(:,:,1) * tamp_srf(:,:,1) / deno +    &
      &            tmp_rain(:,:,2) * tamp_srf(:,:,2) / deno
        wri_snow = tmp_snow(:,:,1) * tamp_srf(:,:,1) / deno +    &
      &            tmp_snow(:,:,2) * tamp_srf(:,:,2) / deno
        wri_rriv = tmp_rriv(:,:,1) * tamp_srf(:,:,1) / deno +    &
      &            tmp_rriv(:,:,2) * tamp_srf(:,:,2) / deno
        wri_rcoa = tmp_rcoa(:,:,1) * tamp_srf(:,:,1) / deno +    &
      &            tmp_rcoa(:,:,2) * tamp_srf(:,:,2) / deno
        wri_taux = tmp_taux(:,:,1) * tamp_srf(:,:,1) / deno +    &
      &            tmp_taux(:,:,2) * tamp_srf(:,:,2) / deno
        wri_tauy = tmp_tauy(:,:,1) * tamp_srf(:,:,1) / deno +    &
      &            tmp_tauy(:,:,2) * tamp_srf(:,:,2) / deno
      endwhere

      call intocpl(itime, (jjm+1)*iim, wri_sol_ice, wri_sol_sea, wri_nsol_ice,&
      & wri_nsol_sea, wri_fder_ice, wri_evap_ice, wri_evap_sea, wri_rain, &
      & wri_snow, wri_rcoa, wri_rriv, wri_taux, wri_tauy, wri_taux, wri_tauy, &
      & lafin ) 
      cpl_sols = 0.; cpl_nsol = 0.; cpl_rain = 0.; cpl_snow = 0.
      cpl_evap = 0.; cpl_tsol = 0.; cpl_fder = 0.; cpl_albe = 0.
      cpl_taux = 0.; cpl_tauy = 0.; cpl_rriv = 0.; cpl_rcoa = 0.
!
! deallocation memoire variables temporaires
!
      sum_error = 0
      deallocate(tmp_sols, stat=error); sum_error = sum_error + error
      deallocate(tmp_nsol, stat=error); sum_error = sum_error + error
      deallocate(tmp_rain, stat=error); sum_error = sum_error + error
      deallocate(tmp_snow, stat=error); sum_error = sum_error + error
      deallocate(tmp_evap, stat=error); sum_error = sum_error + error
      deallocate(tmp_fder, stat=error); sum_error = sum_error + error
      deallocate(tmp_tsol, stat=error); sum_error = sum_error + error
      deallocate(tmp_albe, stat=error); sum_error = sum_error + error
      deallocate(tmp_taux, stat=error); sum_error = sum_error + error
      deallocate(tmp_tauy, stat=error); sum_error = sum_error + error
      deallocate(tmp_rriv, stat=error); sum_error = sum_error + error
      deallocate(tmp_rcoa, stat=error); sum_error = sum_error + error
      if (sum_error /= 0) then
        abort_message='Pb deallocation variables couplees'
        call abort_gcm(modname,abort_message,1)
      endif

    endif

  endif            ! fin (mod(itime, nexca) == 0)
!
! on range les variables lues/sauvegardees dans les bonnes variables de sortie
!
  if (nisurf == is_oce) then
    call cpl2gath(read_sst, tsurf_new, klon, knon,iim,jjm, knindex)
  else if (nisurf == is_sic) then
    call cpl2gath(read_sit, tsurf_new, klon, knon,iim,jjm, knindex)
    call cpl2gath(read_alb_sic, alb_new, klon, knon,iim,jjm, knindex)
  endif
  pctsrf_new(:,nisurf) = pctsrf_sav(:,nisurf)
  
!  if (lafin) call quitcpl

  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, pctsrf_new)

! Cette routine sert d'interface entre le modele atmospherique et un fichier
! de conditions aux limites
!
! 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
!   pctsrf_new   sous-maille fractionnelle
!


! Parametres d'entree
  integer, intent(IN) :: itime
  real   , intent(IN) :: dtime
  integer, intent(IN) :: jour
  integer, intent(IN) :: nisurf
  integer, intent(IN) :: knon
  integer, intent(IN) :: klon
  integer, dimension(klon), intent(in) :: knindex
  logical, intent(IN) :: debut

! Parametres de sortie
  real, intent(out), dimension(klon) :: lmt_sst
  real, intent(out), dimension(klon,nbsrf) :: pctsrf_new

! Variables locales
  integer     :: ii
  INTEGER,save :: 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.nc'
  LOGICAL     :: newlmt = .TRUE.
  logical     :: check = .true.
! Champs lus dans le fichier de CL
  real, allocatable , save, dimension(:) :: sst_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 .and. .not. allocated(sst_lu)) then
    lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
    jour_lu = jour - 1
    allocate(sst_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 
  if (check) write(*,*)modname,' :: itime, lmt_pas ', itime, lmt_pas,dtime

! 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
!
    fich = trim(fich)
    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,pct_tmp(1,is_oce))
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(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,pct_tmp(1,is_sic))
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(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,pct_tmp(1,is_ter))
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(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,pct_tmp(1,is_lic))
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(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    

!
! 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))
  enddo
! je peux pas utiliser la ligne suivante a cause du compilo Sun
!  lmt_sst = sst_lu(knindex)
  pctsrf_new = pct_tmp

  END SUBROUTINE interfoce_lim

!
!#########################################################################
!
  SUBROUTINE interfsur_lim(itime, dtime, jour, &
     & klon, nisurf, knon, knindex, &
     & debut,  &
     & lmt_alb, lmt_rug)

! Cette routine sert d'interface entre le modele atmospherique et un fichier
! de conditions aux limites
!
! 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
!


! Parametres d'entree
  integer, intent(IN) :: itime
  real   , intent(IN) :: dtime
  integer, intent(IN) :: jour
  integer, intent(IN) :: nisurf
  integer, intent(IN) :: knon
  integer, intent(IN) :: klon
  integer, dimension(klon), intent(in) :: knindex
  logical, intent(IN) :: debut

! Parametres de sortie
  real, intent(out), dimension(klon) :: lmt_alb
  real, intent(out), dimension(klon) :: lmt_rug

! Variables locales
  integer     :: ii
  integer,save :: lmt_pas     ! frequence de lecture des conditions limites 
                             ! (en pas de physique)
  logical,save :: deja_lu_sur! pour indiquer que le jour a lire a deja
                             ! lu pour une surface precedente
  integer,save :: jour_lu_sur 
  integer      :: ierr
  character (len = 20) :: modname = 'interfsur_lim'
  character (len = 80) :: abort_message
  character (len = 20) :: fich ='limit.nc'
  logical     :: newlmt = .false.
  logical     :: check = .true.
! Champs lus dans le fichier de CL
  real, allocatable , save, dimension(:) :: alb_lu, rug_lu
!
! 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_sur = jour - 1
    allocate(alb_lu(klon))
    allocate(rug_lu(klon))
  endif

  if ((jour - jour_lu_sur) /= 0) deja_lu_sur = .false.
  
  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)

! Tester d'abord si c'est le moment de lire le fichier
  if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu_sur) then
!
! Ouverture du fichier
!
    fich = trim(fich)
    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
!
! 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_sur = .true.
    jour_lu_sur = jour
  endif
!
! Recopie des variables dans les champs de sortie
!
  DO ii = 1, knon
    lmt_alb(ii) = alb_lu(knindex(ii))
    lmt_rug(ii) = rug_lu(knindex(ii))
  enddo

  END SUBROUTINE interfsur_lim

!
!#########################################################################
!

  SUBROUTINE calcul_fluxs( klon, knon, nisurf, dtime, &
     & tsurf, p1lay, cal, beta, coef1lay, ps, &
     & precip_rain, precip_snow, snow, qsol, &
     & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
     & petAcoef, peqAcoef, petBcoef, peqBcoef, &
     & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)

! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
! une temperature de surface (au cas ou ok_veget = false)
!
! L. Fairhead 4/2000
!
! input:
!   knon         nombre de points a traiter
!   nisurf       surface a traiter
!   tsurf        temperature de surface
!   p1lay        pression 1er niveau (milieu de couche)
!   cal          capacite calorifique du sol
!   beta         evap reelle
!   coef1lay     coefficient d'echange
!   ps           pression au sol
!   precip_rain  precipitations liquides
!   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
!   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
!   radsol       rayonnement net aus sol (LW + SW)
!   dif_grnd     coeff. diffusion vers le sol profond
!
! output:
!   tsurf_new    temperature au sol
!   fluxsens     flux de chaleur sensible
!   fluxlat      flux de chaleur latente
!   dflux_s      derivee du flux de chaleur sensible / Ts
!   dflux_l      derivee du flux de chaleur latente  / Ts
!

#include "YOETHF.inc"
#include "FCTTRE.inc"

! Parametres d'entree
  integer, intent(IN) :: knon, nisurf, klon
  real   , intent(IN) :: dtime
  real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
  real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
  real, dimension(klon), intent(IN) :: ps, q1lay
  real, dimension(klon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
  real, dimension(klon), intent(IN) :: precip_rain, precip_snow
  real, dimension(klon), intent(IN) :: radsol, dif_grnd
  real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay
  real, dimension(klon), intent(INOUT) :: snow, qsol

! Parametres sorties
  real, dimension(klon), intent(OUT):: tsurf_new, evap, fluxsens, fluxlat
  real, dimension(klon), intent(OUT):: dflux_s, dflux_l

! Variables locales
  integer :: i
  real, dimension(klon) :: zx_mh, zx_nh, zx_oh
  real, dimension(klon) :: zx_mq, zx_nq, zx_oq
  real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
  real, dimension(klon) :: zx_sl, zx_k1,  zx_dq,  zx_cq,  zx_dh, zx_ch
  real, dimension(klon) :: zx_h_ts, zx_q_0 , d_ts
  real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
  real                  :: bilan_f, fq_fonte
  real, parameter :: t_grnd = 271.35, t_coup = 273.15
  logical         :: check = .false.
  character (len = 20)  :: modname = 'calcul_fluxs'
  logical         :: fonte_neige = .false.
  real            :: max_eau_sol = 150.0
  character (len = 80) :: abort_message 

  if (check) write(*,*)'Entree ', modname,' surface = ',nisurf

  if (size(run_off) /= knon .AND. nisurf == is_ter) then
    write(*,*)'Bizarre, le nombre de points continentaux'
    write(*,*)'a change entre deux appels. J''arrete ...'
    abort_message='Pb run_off'
    call abort_gcm(modname,abort_message,1)
  endif
!
! Traitement neige et humidite du sol
!
    if (nisurf == is_oce) then
      snow = 0.
      qsol = max_eau_sol
    else
      snow = max(0.0, snow + (precip_snow - evap) * dtime)
      qsol = qsol + (precip_rain - evap) * dtime
    endif 


! 
! Initialisation
!
!
! zx_qs = qsat en kg/kg
!
  DO i = 1, knon
    zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
  IF (thermcep) THEN
      zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
      zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
      zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
      zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
      zx_qs=MIN(0.5,zx_qs)
      zcor=1./(1.-retv*zx_qs)
      zx_qs=zx_qs*zcor
      zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
     &                 /RLVTT / zx_pkh(i)
    ELSE
      IF (tsurf(i).LT.t_coup) THEN
        zx_qs = qsats(tsurf(i)) / ps(i)
        zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
     &                    / zx_pkh(i)
      ELSE
        zx_qs = qsatl(tsurf(i)) / ps(i)
        zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
     &               / zx_pkh(i)
      ENDIF
    ENDIF
    zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
    zx_qsat(i) = zx_qs
    zx_coef(i) = coef1lay(i) &
     & * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
     & * p1lay(i)/(RD*t1lay(i))
    
  ENDDO


! === Calcul de la temperature de surface ===
! 
! zx_sl = chaleur latente d'evaporation ou de sublimation
!
  do i = 1, knon
    zx_sl(i) = RLVTT
    if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
    zx_k1(i) = zx_coef(i)
  enddo


  do i = 1, knon
! Q
    zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
    zx_mq(i) = beta(i) * zx_k1(i) * &
     &             (peqAcoef(i) - zx_qsat(i) &
     &                          + zx_dq_s_dt(i) * tsurf(i)) &
     &             / zx_oq(i)
    zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
     &                              / zx_oq(i)

! H
    zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
    zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
    zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)

! Tsurface
    tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
     &             (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) & 
     &                 + dif_grnd(i) * t_grnd * dtime)/ &
     &          ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
     &                       zx_nh(i) + zx_sl(i) * zx_nq(i)) &  
     &                     + dtime * dif_grnd(i))

!
! Y'a-t-il fonte de neige?
!
    fonte_neige = (nisurf /= is_oce) .AND. &
     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
     & .AND. (tsurf_new(i) >= RTT)
    if (fonte_neige) tsurf_new(i) = RTT  
    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
    d_ts(i) = tsurf_new(i) - tsurf(i)
    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
    evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i) 
    fluxlat(i) = - evap(i) * zx_sl(i)
    fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
! Derives des flux dF/dTs (W m-2 K-1):
    dflux_s(i) = zx_nh(i)
    dflux_l(i) = (zx_sl(i) * zx_nq(i))
!
! en cas de fonte de neige
!
    if (fonte_neige) then
      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
      bilan_f = max(0., bilan_f)
      fq_fonte = bilan_f / zx_sl(i)
      snow(i) = max(0., snow(i) - fq_fonte * dtime)
      qsol(i) = qsol(i) + (fq_fonte * dtime)
    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) 
  ENDDO

  END SUBROUTINE calcul_fluxs
!
!#########################################################################
!
  SUBROUTINE gath2cpl(champ_in, champ_out, klon, knon, iim, jjm, knindex)

! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
! au coupleur.
!
! 
! input:         
!   champ_in     champ sur la grille gathere        
!   knon         nombre de points dans le domaine a traiter
!   knindex      index des points de la surface a traiter
!   klon         taille de la grille
!   iim,jjm      dimension de la grille 2D
!
! output:
!   champ_out    champ sur la grille 2D
!
! input
  integer                   :: klon, knon, iim, jjm
  real, dimension(klon)     :: champ_in
  integer, dimension(klon)  :: knindex
! output
  real, dimension(iim,jjm+1)  :: champ_out
! local
  integer                   :: i, ig, j
  real, dimension(klon)     :: tamp

  tamp = 0.
  do i = 1, knon
    ig = knindex(i)
    tamp(ig) = champ_in(i)
  enddo    
  ig = 1
  champ_out(:,1) = tamp(ig)
  do j = 2, jjm
    do i = 1, iim
      ig = ig + 1
      champ_out(i,j) = tamp(ig)
    enddo
  enddo
  ig = ig + 1
  champ_out(:,jjm+1) = tamp(ig)

  END SUBROUTINE gath2cpl
!
!#########################################################################
!
  SUBROUTINE cpl2gath(champ_in, champ_out, klon, knon, iim, jjm, knindex)

! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
! au coupleur.
!
! 
! input:         
!   champ_in     champ sur la grille gathere        
!   knon         nombre de points dans le domaine a traiter
!   knindex      index des points de la surface a traiter
!   klon         taille de la grille
!   iim,jjm      dimension de la grille 2D
!
! output:
!   champ_out    champ sur la grille 2D
!
! input
  integer                   :: klon, knon, iim, jjm
  real, dimension(iim,jjm+1)     :: champ_in
  integer, dimension(klon)  :: knindex
! output
  real, dimension(klon)  :: champ_out
! local
  integer                   :: i, ig, j
  real, dimension(klon)     :: tamp
  logical                   :: check = .false.

  ig = 1
  tamp(ig) = champ_in(1,1)
  do j = 2, jjm
    do i = 1, iim
      ig = ig + 1
      tamp(ig) = champ_in(i,j)
    enddo
  enddo
  ig = ig + 1
  tamp(ig) = champ_in(1,jjm+1)

  do i = 1, knon
    ig = knindex(i)
    champ_out(i) = tamp(ig)
  enddo    

  END SUBROUTINE cpl2gath
!
!#########################################################################
!
  SUBROUTINE albsno(klon, agesno,alb_neig_grid)
  IMPLICIT none
 
  integer :: klon
  INTEGER, PARAMETER :: nvm = 8
  REAL, dimension(klon,nvm) :: veget
  REAL, DIMENSION(klon) :: alb_neig_grid, agesno
 
  INTEGER :: i, nv
 
  REAL, DIMENSION(nvm),SAVE :: init, decay
  REAL :: as
  DATA init /0.55, 0.14, 0.18, 0.29, 0.15, 0.15, 0.14, 0./
  DATA decay/0.30, 0.67, 0.63, 0.45, 0.40, 0.14, 0.06, 1./
 
  veget = 0.
  veget(:,1) = 1.     ! desert partout
  DO i = 1, klon
    alb_neig_grid(i) = 0.0
  ENDDO
  DO nv = 1, nvm
    DO i = 1, klon
      as = init(nv)+decay(nv)*EXP(-agesno(i)/5.)
      alb_neig_grid(i) = alb_neig_grid(i) + veget(i,nv)*as
    ENDDO
  ENDDO
 
  END SUBROUTINE albsno
!
!#########################################################################
!
  END MODULE interface_surf