source: LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90 @ 457

! $Header$

  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

  IMPLICIT none

  PRIVATE
  PUBLIC :: interfsurf,interfsurf_hq, gath2cpl 

  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
  real, allocatable, dimension(:),save    :: coastalflow, riverflow
!!$PB
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: tmp_rriv, tmp_rcoa
!! pour simuler la fonte des glaciers antarctiques
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: coeff_iceberg
  real, save                              :: surf_maille 
  real, save                              :: cte_flux_iceberg = 6.3e7
  integer, save                           :: num_antarctic = 1
!!$
  CONTAINS
!
!############################################################################
!
  SUBROUTINE interfsurf_hq(itime, dtime, date0, jour, rmu0, &
      & klon, iim, jjm, nisurf, knon, knindex, pctsrf, &
      & rlon, rlat, cufi, cvfi,&
      & debut, lafin, ok_veget, soil_model, nsoilmx, tsoil, qsol,&
      & zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, & 
      & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
      & precip_rain, precip_snow, sollw, sollwdown, swnet, swdown, &
      & fder, taux, tauy, rugos, rugoro, &
      & albedo, snow, qsurf, &
      & tsurf, p1lay, ps, radsol, &
      & ocean, npas, nexca, zmasq, &
      & evap, fluxsens, fluxlat, dflux_l, dflux_s, &              
      & tsol_rad, tsurf_new, alb_new, alblw, emis_new, &
!IM cf. JLD      & z0_new, pctsrf_new, agesno)
      & z0_new, pctsrf_new, agesno,fqcalving,ffonte)


! 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)
!   date0        jour initial 
!   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
!   cufi,cvfi    resolution des mailles en x et y (m)
!   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
!   epot_air     temp potentielle 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
!   sollw        flux IR net a la surface
!   sollwdown    flux IR descendant 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
!   rugos        rugosite
!   zmasq        masque terre/ocean
!   rugoro       rugosite orographique
!
! 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
  real, intent(IN) :: date0
  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) :: cufi, cvfi
  real, dimension(klon), intent(INOUT) :: tq_cdrag
  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(IN) :: petAcoef, peqAcoef
  real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
  real, dimension(klon), intent(IN) :: precip_rain, precip_snow
  real, dimension(klon), intent(IN) :: sollw, sollwdown, swnet, swdown
  real, dimension(klon), intent(IN) :: ps, albedo
!IM cf LF
! real, dimension(klon), intent(INOUT) :: tsurf
! real, dimension(klon), intent(IN) :: p1lay
  real, dimension(klon), intent(IN) :: tsurf, p1lay
  REAL, DIMENSION(klon), INTENT(INOUT) :: radsol,fder
  real, dimension(klon), intent(IN) :: zmasq
  real, dimension(klon), intent(IN) :: taux, tauy, rugos, rugoro
  character (len = 6)  :: ocean
  integer              :: npas, nexca ! nombre et pas de temps couplage
  real, dimension(klon), intent(INOUT) :: evap, snow, qsurf
!! PB ajout pour soil
  logical          :: soil_model
  integer          :: nsoilmx
  REAL, DIMENSION(klon, nsoilmx) :: tsoil
  REAL, dimension(klon), intent(INOUT) :: qsol
  REAL, dimension(klon)          :: soilcap
  REAL, dimension(klon)          :: soilflux
! 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):: alblw
  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

! Flux thermique utiliser pour fondre la neige
!jld a rajouter   real, dimension(klon), intent(INOUT):: ffonte
  real, dimension(klon), intent(INOUT):: ffonte
! Flux d'eau "perdue" par la surface et nécessaire pour que limiter la
! hauteur de neige, en kg/m2/s
!jld a rajouter   real, dimension(klon), intent(INOUT):: fqcalving
  real, dimension(klon), intent(INOUT):: fqcalving

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

  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
!    allocate(alb_neig_grid(klon), stat = error)
!    if (error /= 0) then
!      abort_message='Pb allocation alb_neig_grid'
!      call abort_gcm(modname,abort_message,1)
!    endif
  endif
  first_call = .false.
  
! Initialisations diverses
!
!!$  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.
  alblw = 999999.
! 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(coastalflow)) then
      allocate(coastalflow(knon), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation coastalflow'
        call abort_gcm(modname,abort_message,1)
      endif
      allocate(riverflow(knon), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation riverflow'
        call abort_gcm(modname,abort_message,1)
      endif
      allocate(run_off(knon), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation run_off'
        call abort_gcm(modname,abort_message,1)
      endif
!!$PB
      ALLOCATE (tmp_rriv(iim,jjm+1), stat=error)
      if (error /= 0) then
        abort_message='Pb allocation tmp_rriv'
        call abort_gcm(modname,abort_message,1)
      endif
      ALLOCATE (tmp_rcoa(iim,jjm+1), stat=error)
      if (error /= 0) then
        abort_message='Pb allocation tmp_rcoa'
        call abort_gcm(modname,abort_message,1)
      endif
!!$
    else if (size(coastalflow) /= knon) then
      write(*,*)'Bizarre, le nombre de points continentaux'
      write(*,*)'a change entre deux appels. J''arrete ...'
      abort_message='voir ci-dessus'
      call abort_gcm(modname,abort_message,1)
    endif
    coastalflow = 0.
    riverflow = 0.
!
! Calcul age de la neige
!
!!$ PB ATTENTION changement ordre des appels
!!$    CALL albsno(klon,agesno,alb_neig_grid)  

    if (.not. ok_veget) then
!
! 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)
!  
! calcul snow et qsurf, hydrol adapté
!
       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)
       ELSE 
           cal = RCPD * capsol
!!$      cal = capsol
       ENDIF
       CALL calcul_fluxs( klon, knon, nisurf, dtime, &
     &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &
     &   precip_rain, precip_snow, snow, qsurf,  &
     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)

       CALL fonte_neige( 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, &
!IM cf JLD
     &   fqcalving,ffonte)


     call albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
     where (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
     zfra(1:knon) = max(0.0,min(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
     alb_new(1 : knon)  = alb_neig(1 : knon) *zfra(1:knon) + &
    &                     alb_new(1 : knon)*(1.0-zfra(1:knon))
     z0_new = sqrt(z0_new**2+rugoro**2)
     alblw(1 : knon) = alb_new(1 : knon)

    else
!!      CALL albsno(klon,agesno,alb_neig_grid)  
!
!  appel a sechiba
!
      call interfsol(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, sollwdown, swnet, swdown, &
     &  tsurf, p1lay/100., ps/100., radsol, &
     &  evap, fluxsens, fluxlat, &              
     &  tsol_rad, tsurf_new, alb_new, alblw, &
     &  emis_new, z0_new, dflux_l, dflux_s, qsurf_new)

!  
! ajout de la contribution du relief
!  
      z0_new = SQRT(z0_new**2+rugoro**2)
!
! mise a jour de l'humidite saturante calculee par ORCHIDEE
      qsurf(1:knon) = qsurf_new(1:knon)

    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
      if (nexca == 0) then
        abort_message='nexca = 0 dans interfoce_cpl'
        call abort_gcm(modname,abort_message,1)
      endif

      cumul = .false.

      call interfoce(itime, dtime, cumul, &
      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
      & ocean, npas, nexca, debut, lafin, &
      & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
      & fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
      & tsurf_new, alb_new, 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.
    alb_neig(:) = 0.
    agesno(:) = 0.

    call calcul_fluxs( klon, knon, nisurf, dtime, &
     &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
     &   precip_rain, precip_snow, snow, qsurf,  &
     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
    
    fder = fder + dflux_s + dflux_l

!
! 2eme appel a interfoce pour le cumul des champs (en particulier
! fluxsens et fluxlat calcules dans calcul_fluxs)
!
    if (ocean == 'couple') then

      cumul = .true.

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

!    else if (ocean == 'slab  ') then
!      call interfoce(nisurf)

    endif

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

    z0_new = sqrt(rugos**2 + rugoro**2)
    alblw(1:knon) = alb_new(1:knon)

!
  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

      cumul =.false.

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

      tsurf_temp = tsurf_new
      cal = 0.
      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)

!IM cf LF
      DO i = 1, knon
       IF (pctsrf_new(i,nisurf) < EPSFRA) then
          snow(i) = 0.0
!IM cf LF/JLD         tsurf(i) = RTT - 1.8
          tsurf_new(i) = RTT - 1.8
          IF (soil_model) tsoil(i,:) = RTT -1.8
        endif
      enddo

      CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
      
      IF (soil_model) THEN 
!IM cf LF/JLD        CALL soil(dtime, nisurf, knon,snow, tsurf, tsoil,soilcap, soilflux)
         CALL soil(dtime, nisurf, knon,snow, tsurf_new, tsoil,soilcap, soilflux)
         cal(1:knon) = RCPD / soilcap(1:knon)
         radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
         dif_grnd = 0.
      ELSE 
         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, qsurf,  &
         &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
         &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
         &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)

!IM cf JP 12.02.03
!     DO i = 1, knon
!      IF (pctsrf_new(i,nisurf) < EPSFRA) then
!         snow(i) = 0.0      
!         tsurf_new(i) = RTT - 1.8
!         IF (soil_model) tsoil(i,:) = RTT -1.8
!       endif
!     enddo

    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, &
!IM cf JLD
             &   fqcalving,ffonte)

!     calcul albedo

      CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
      WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
      zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
      alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + & 
     &                    0.6 * (1.0-zfra(1:knon))
!!      alb_new(1 : knon) = 0.6
    ENDIF

    fder = fder + dflux_s + dflux_l

!
! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
!
    if (ocean == 'couple') then

      cumul =.true.

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

!    else if (ocean == 'slab  ') then
!      call interfoce(nisurf)

    endif

     
    z0_new = 0.001
    z0_new = SQRT(z0_new**2+rugoro**2)
    alblw(1:knon) = alb_new(1:knon)

  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)
    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)
    ELSE 
        cal = RCPD * calice
        WHERE (snow > 0.0) cal = RCPD * calsno
    ENDIF 
    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, qsurf,  &
     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)

    call fonte_neige( 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, &
!IM cf JLD
     &   fqcalving,ffonte)

!
! calcul albedo
!
     CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
     WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
     zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
     alb_new(1 : knon)  = alb_neig(1 : knon)*zfra(1:knon) + &
    &                     0.6 * (1.0-zfra(1:knon))
!!     alb_new(1 : knon)  = 0.6
!
! Rugosite
!
    z0_new = rugoro
!
! Remplissage des pourcentages de surface
!
    pctsrf_new(:,nisurf) = pctsrf(:,nisurf)

    alblw(1:knon) = alb_new(1:knon)
  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, date0, nisurf, knon, &
     & knindex, rlon, rlat, cufi, cvfi, iim, jjm, pctsrf, &
     & debut, lafin, ok_veget, &
     & plev,  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, alblw, &
     & emis_new, z0_new, dflux_l, dflux_s, qsurf)

  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)
!   plev         hauteur de la premiere couche (Pa)      
!   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
!   qsurf        air moisture at surface

! 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) :: plev
  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
!IM cf. JP +++
  real, dimension(klon) :: swdown_vrai
!IM cf. JP ---
  real, dimension(klon), intent(IN) :: tsurf, p1lay
  real, dimension(klon), intent(IN) :: radsol
! Parametres de sortie
  real, dimension(klon), intent(OUT):: evap, fluxsens, fluxlat, qsurf
  real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new, alblw
  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, iglob
  integer              :: error
  character (len = 20) :: modname = 'interfsol'
  character (len = 80) :: abort_message
  logical,save              :: check = .FALSE.
  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
! Preserved albedo
!IM cf. JP +++
  real, allocatable, dimension(:), save :: albedo_keep, zlev
!IM cf. JP ---
! 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):: snow
  real, dimension(knon,2) :: albedo_out
! Pb de nomenclature
  real, dimension(klon) :: petA_orc, peqA_orc
  real, dimension(klon) :: petB_orc, peqB_orc
! Pb de correspondances de grilles
  integer, dimension(:), save, allocatable :: ig, jg
  integer :: indi, indj
  integer, dimension(klon) :: ktindex
  REAL, dimension(klon) :: bidule
! Essai cdrag
  real, dimension(klon) :: cdrag

#include "temps.inc"
#include "YOMCST.inc"

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

  ktindex(:) = knindex(:) + iim - 1

! initialisation
  if (debut) then

  IF ( .NOT. allocated(albedo_keep)) THEN
     ALLOCATE(albedo_keep(klon))
     ALLOCATE(zlev(klon))
  ENDIF
! Pb de correspondances de grilles
   allocate(ig(klon))
   allocate(jg(klon))
   ig(1) = 1
   jg(1) = 1
   indi = 0
   indj = 2
   do igrid = 2, klon - 1
     indi = indi + 1
     if ( indi > iim) then
       indi = 1
       indj = indj + 1
     endif
     ig(igrid) = indi
     jg(igrid) = indj
   enddo
   ig(klon) = 1
   jg(klon) = jjm + 1
!
!  Initialisation des offset    
!
! offset bord ouest
   off_ini(1,1) = - iim  ; off_ini(2,1) = - iim + 1; off_ini(3,1) = 1
   off_ini(4,1) = iim + 1; off_ini(5,1) = iim      ; off_ini(6,1) = 2 * iim - 1
   off_ini(7,1) = iim -1 ; off_ini(8,1) = - 1
! offset point normal
   off_ini(1,2) = - iim  ; off_ini(2,2) = - iim + 1; off_ini(3,2) = 1
   off_ini(4,2) = iim + 1; off_ini(5,2) = iim      ; off_ini(6,2) = iim - 1
   off_ini(7,2) = -1     ; off_ini(8,2) = - iim - 1
! offset bord   est
   off_ini(1,3) = - iim; off_ini(2,3) = - 2 * iim + 1; off_ini(3,3) = - iim + 1
   off_ini(4,3) =  1   ; off_ini(5,3) = iim          ; off_ini(6,3) = iim - 1
   off_ini(7,3) = -1   ; off_ini(8,3) = - iim - 1
!
! Initialisation des correspondances point -> indices i,j
!
    if (( .not. allocated(correspond))) then
      allocate(correspond(iim,jjm+1), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation correspond'
        call abort_gcm(modname,abort_message,1)
      endif      
    endif
!
! Attention aux poles
!
    do igrid = 1, knon
      index = ktindex(igrid)
          jj = int((index - 1)/iim) + 1
          ij = index - (jj - 1) * iim
      correspond(ij,jj) = igrid
    enddo

! 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
    if ((.not. allocated(lon_scat))) then
      allocate(lon_scat(iim,jjm+1), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation lon_scat'
        call abort_gcm(modname,abort_message,1)
      endif      
    endif
    if ((.not. allocated(lat_scat))) then
      allocate(lat_scat(iim,jjm+1), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation lat_scat'
        call abort_gcm(modname,abort_message,1)
      endif      
    endif
    lon_scat = 0.
    lat_scat = 0.
    do igrid = 1, knon
      index = knindex(igrid)
      lalo(igrid,2) = rlon(index)
      lalo(igrid,1) = rlat(index)
      ij = index - int((index-1)/iim)*iim - 1
      jj = 2 + int((index-1)/iim)
      if (mod(index,iim) == 1 ) then
        jj = 1 + int((index-1)/iim)
        ij = iim
      endif
!      lon_scat(ij,jj) = rlon(index)
!      lat_scat(ij,jj) = rlat(index)
    enddo
    index = 1
    do jj = 2, jjm
      do ij = 1, iim
        index = index + 1
        lon_scat(ij,jj) = rlon(index)
        lat_scat(ij,jj) = rlat(index)
      enddo
    enddo
    lon_scat(:,1) = lon_scat(:,2)
    lat_scat(:,1) = rlat(1)
    lon_scat(:,jjm+1) = lon_scat(:,2)
    lat_scat(:,jjm+1) = rlat(klon)
! Pb de correspondances de grilles!
!    do igrid = 1, knon
!      index = ktindex(igrid)
!      ij = ig(index)
!      jj = jg(index)
!      lon_scat(ij,jj) = rlon(index)
!      lat_scat(ij,jj) = rlat(index)
!    enddo

!
! Allouer et initialiser le tableau des voisins et des fraction de continents
!
    if ( (.not.allocated(neighbours))) THEN
      allocate(neighbours(knon,8), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation neighbours'
        call abort_gcm(modname,abort_message,1)
      endif
    endif
    neighbours = -1.
    if (( .not. allocated(contfrac))) then
      allocate(contfrac(knon), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation contfrac'
        call abort_gcm(modname,abort_message,1)
      endif      
    endif

    do igrid = 1, knon
      ireal = knindex(igrid)
      contfrac(igrid) = pctsrf(ireal,is_ter)
    enddo

    do igrid = 1, knon
      iglob = ktindex(igrid)
      if (mod(iglob, iim) == 1) then
        offset = off_ini(:,1)
      else if(mod(iglob, iim) == 0) then
        offset = off_ini(:,3)
      else
        offset = off_ini(:,2)
      endif 
      do i = 1, 8
        index = iglob + offset(i)
        ireal = (min(max(1, index - iim + 1), klon))
        if (pctsrf(ireal, is_ter) > EPSFRA) then
          jj = int((index - 1)/iim) + 1
          ij = index - (jj - 1) * iim
            neighbours(igrid, i) = correspond(ij, jj)
        endif
      enddo
    enddo

!
!  Allocation et calcul resolutions
    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
    do igrid = 1, knon
      ij = knindex(igrid)
      resolution(igrid,1) = cufi(ij)
      resolution(igrid,2) = cvfi(ij)
    enddo  

  endif                          ! (fin debut) 

! 
! Appel a la routine sols continentaux
!
  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

  cdrag = 0.
  cdrag(1:knon) = tq_cdrag(1:knon)

!IM cf. JP +++
! zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/287.05*temp_air(1:knon))*9.80665)
  zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/RD*temp_air(1:knon))*RG)
!IM cf. JP ---

   where(cdrag > 0.01) 
     cdrag = 0.01
   endwhere
!  write(*,*)'Cdrag = ',minval(cdrag),maxval(cdrag)

!
! Init Orchidee
!
  if (debut) then
    call intersurf_main (itime+itau_phy-1, iim, jjm+1, knon, ktindex, dtime, &
     & lrestart_read, lrestart_write, lalo, &
     & contfrac, neighbours, resolution, date0, &
     & zlev,  u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
     & cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
     & precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
     & evap, fluxsens, fluxlat, coastalflow, riverflow, &
     & tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
     & lon_scat, lat_scat)

!IM cf. JP +++
    albedo_keep(:) = (albedo_out(:,1)+albedo_out(:,2))/2.
!IM cf. JP ---

  endif

!IM cf. JP +++
  swdown_vrai(:) = swnet(:)/(1. - albedo_keep(:))
!IM cf. JP ---

  call intersurf_main (itime+itau_phy, iim, jjm+1, knon, ktindex, dtime, &
     & lrestart_read, lrestart_write, lalo, &
     & contfrac, neighbours, resolution, date0, &
     & zlev,  u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
     & cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
!IM cf. JP +++
     & precip_rain, precip_snow, lwdown, swnet, swdown_vrai, ps, &
!IM cf. JP ---
     & evap, fluxsens, fluxlat, coastalflow, riverflow, &
     & tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
     & lon_scat, lat_scat)

!IM cf. JP +++
    albedo_keep(:) = (albedo_out(:,1)+albedo_out(:,2))/2.
!IM cf. JP ---

    bidule=0.
    bidule(1:knon)=riverflow(1:knon)
    call gath2cpl(bidule, tmp_rriv, klon, knon,iim,jjm,knindex)
    bidule=0.
    bidule(1:knon)=coastalflow(1:knon)
    call gath2cpl(bidule, tmp_rcoa, klon, knon,iim,jjm,knindex)
    alb_new(1:knon) = albedo_out(1:knon,1) 
    alblw(1:knon) = albedo_out(1:knon,2)


! Convention orchidee: positif vers le haut
  fluxsens(1:knon) = -1. * fluxsens(1:knon)
  fluxlat(1:knon)  = -1. * fluxlat(1:knon)

!  evap     = -1. * evap

  if (debut) lrestart_read = .false.

  END SUBROUTINE interfsol
!
!#########################################################################
!
  SUBROUTINE interfoce_cpl(itime, dtime, cumul, &
      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
      & ocean, npas, nexca, debut, lafin, &
      & swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
      & fluxlat, fluxsens, fder, albsol, taux, tauy, zmasq, &
      & tsurf_new, alb_new, 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 temps de 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 net 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(IN) :: fluxlat, fluxsens
  logical, intent(IN)               :: cumul
  real, dimension(klon), intent(INOUT) :: evap

! Parametres de sortie
  real, dimension(klon), intent(OUT):: tsurf_new, alb_new
  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,save              :: check = .FALSE.
! 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
!!$PB  real, allocatable, dimension(:,:),save :: cpl_tauy, cpl_rriv, cpl_rcoa
  real, allocatable, dimension(:,:),save :: cpl_tauy
  real, allocatable, dimension(:,:),save :: 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
  REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: tmp_tauy
! 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, wri_rcoa, wri_rriv
  REAL, DIMENSION(iim, jjm+1) :: wri_rain, wri_snow, wri_taux, wri_tauy
  REAL, DIMENSION(iim, jjm+1) :: wri_calv
  REAL, DIMENSION(iim, jjm+1) :: wri_tauxx, wri_tauyy, wri_tauzz
  REAL, DIMENSION(iim, jjm+1) :: tmp_lon, tmp_lat
! 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_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(:),save :: unity
! 
  logical, save    :: first_appel = .true.
  logical,save          :: print
!maf
! variables pour avoir une sortie IOIPSL des champs echanges
  CHARACTER*80,SAVE :: clintocplnam, clfromcplnam
  INTEGER, SAVE :: jf,nhoridct,nidct
  INTEGER, SAVE :: nhoridcs,nidcs
  INTEGER :: ndexct(iim*(jjm+1)),ndexcs(iim*(jjm+1))
  REAL :: zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian
  integer :: idayref, itau_w
#include "param_cou.h"
#include "inc_cpl.h"
#include "temps.inc"
!
! 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
!!$PB
!!$    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(cpl_rriv(iim,jjm+1), stat=error); sum_error = sum_error + error
    ALLOCATE(cpl_rcoa(iim,jjm+1), 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)

!
! initialisation sorties netcdf
!
    idayref = day_ini
    CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
    CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
    DO i = 1, iim
      zx_lon(i,1) = rlon(i+1)
      zx_lon(i,jjm+1) = rlon(i+1)
    ENDDO
    CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
    clintocplnam="cpl_atm_tauflx"
    CALL histbeg(clintocplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),1,iim,1,jjm+1, &
       & itau_phy,zjulian,dtime,nhoridct,nidct) 
! no vertical axis
    CALL histdef(nidct, 'tauxe','tauxe', &
         & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
    CALL histdef(nidct, 'tauyn','tauyn', &
         & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
    CALL histdef(nidct, 'tmp_lon','tmp_lon', &
         & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
    CALL histdef(nidct, 'tmp_lat','tmp_lat', &
         & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
    DO jf=1,jpflda2o1 + jpflda2o2
      CALL histdef(nidct, cl_writ(jf),cl_writ(jf), &
         & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
    END DO
    CALL histend(nidct)
    CALL histsync(nidct)

    clfromcplnam="cpl_atm_sst"
    CALL histbeg(clfromcplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),1,iim,1,jjm+1, &
       & 0,zjulian,dtime,nhoridcs,nidcs) 
! no vertical axis
    DO jf=1,jpfldo2a
      CALL histdef(nidcs, cl_read(jf),cl_read(jf), &
         & "-",iim, jjm+1, nhoridcs, 1, 1, 1, -99, 32, "inst", dtime,dtime)
    END DO
    CALL histend(nidcs)
    CALL histsync(nidcs)

! pour simuler la fonte des glaciers antarctiques
! 
    surf_maille = (4. * rpi * ra**2) / (iim * (jjm +1))
    ALLOCATE(coeff_iceberg(iim,jjm+1), stat=error)
    if (error /= 0) then
      abort_message='Pb allocation variable coeff_iceberg'
      call abort_gcm(modname,abort_message,1)
    endif
    open (12,file='flux_iceberg',form='formatted',status='old')
    read (12,*) coeff_iceberg
    close (12)
    num_antarctic = max(1, count(coeff_iceberg > 0))
    
    first_appel = .false.
  endif ! fin if (first_appel)

! Initialisations

! calcul des fluxs a passer

  cpl_index = 1
  if (nisurf == is_sic) cpl_index = 2
  if (cumul) then
    if (check) write(*,*) modname, 'cumul des champs'
    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) + fluxlat(ig) +fluxsens(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) &
!!$       &                          + riverflow(ig)   / FLOAT(nexca)/dtime
!!$      cpl_rcoa(ig,cpl_index) = cpl_rcoa(ig,cpl_index) &
!!$       &                          + coastalflow(ig) / FLOAT(nexca)/dtime
    enddo
    IF (cpl_index .EQ. 1) THEN 
        cpl_rriv(:,:) = cpl_rriv(:,:) + tmp_rriv(:,:) / FLOAT(nexca)
        cpl_rcoa(:,:) = cpl_rcoa(:,:) + tmp_rcoa(:,:) / FLOAT(nexca)
    ENDIF
  endif

  if (mod(itime, nexca) == 1) then
!
! Demande des champs au coupleur
!
! Si le domaine considere est l'ocean, on lit les champs venant du coupleur
!
    if (nisurf == is_oce .and. .not. cumul) 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)
!
! sorties NETCDF des champs recus
!
      ndexcs(:)=0
      itau_w = itau_phy + itime
      CALL histwrite(nidcs,cl_read(1),itau_w,read_sst,iim*(jjm+1),ndexcs)
      CALL histwrite(nidcs,cl_read(2),itau_w,read_sic,iim*(jjm+1),ndexcs)
      CALL histwrite(nidcs,cl_read(3),itau_w,read_alb_sic,iim*(jjm+1),ndexcs)
      CALL histwrite(nidcs,cl_read(4),itau_w,read_sit,iim*(jjm+1),ndexcs)
      CALL histsync(nidcs)
! pas utile      IF (npas-itime.LT.nexca )CALL histclo(nidcs)

      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
!
! Pour rattraper des erreurs d'arrondis
!
      where (abs(pctsrf_sav(:,is_sic)) .le. 2.*epsilon(pctsrf_sav(1,is_sic)))
        pctsrf_sav(:,is_sic) = 0.
        pctsrf_sav(:,is_oce) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
      endwhere
      where (abs(pctsrf_sav(:,is_oce)) .le. 2.*epsilon(pctsrf_sav(1,is_oce)))
        pctsrf_sav(:,is_sic) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
        pctsrf_sav(:,is_oce) = 0.
      endwhere
      if (minval(pctsrf_sav(:,is_oce)) < 0.) then
        write(*,*)'Pb fraction ocean inferieure a 0'
        write(*,*)'au point ',minloc(pctsrf_sav(:,is_oce))
        write(*,*)'valeur = ',minval(pctsrf_sav(:,is_oce)) 
        abort_message = 'voir ci-dessus'
        call abort_gcm(modname,abort_message,1)
      endif
      if (minval(pctsrf_sav(:,is_sic)) < 0.) then
        write(*,*)'Pb fraction glace inferieure a 0'
        write(*,*)'au point ',minloc(pctsrf_sav(:,is_sic))
        write(*,*)'valeur = ',minval(pctsrf_sav(:,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 .and. (.not. cumul) ) 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 .and. cumul) 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)
!!$PB
      wri_rriv = cpl_rriv(:,:)
      wri_rcoa = cpl_rcoa(:,:)

      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
!!$PB
!!$        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
!
! pour simuler la fonte des glaciers antarctiques
!
!$$$        wri_rain = wri_rain      &
!$$$      &     + coeff_iceberg * cte_flux_iceberg / (num_antarctic * surf_maille)
      wri_calv = coeff_iceberg * cte_flux_iceberg / (num_antarctic * surf_maille)
!
! on passe les coordonnées de la grille
!

      CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,tmp_lon)
      CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,tmp_lat)

      DO i = 1, iim
        tmp_lon(i,1) = rlon(i+1)
        tmp_lon(i,jjm + 1) = rlon(i+1)
      ENDDO
!
! sortie netcdf des champs pour le changement de repere
!
      ndexct(:)=0
      CALL histwrite(nidct,'tauxe',itau_w,wri_taux,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,'tauyn',itau_w,wri_tauy,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,'tmp_lon',itau_w,tmp_lon,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,'tmp_lat',itau_w,tmp_lat,iim*(jjm+1),ndexct)

!
! calcul 3 coordonnées du vent
!
      CALL atm2geo (iim , jjm + 1, wri_taux, wri_tauy, tmp_lon, tmp_lat, &
         & wri_tauxx, wri_tauyy, wri_tauzz )
!
! sortie netcdf des champs apres changement de repere et juste avant
! envoi au coupleur
!
      CALL histwrite(nidct,cl_writ(1),itau_w,wri_sol_ice,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(2),itau_w,wri_sol_sea,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(3),itau_w,wri_nsol_ice,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(4),itau_w,wri_nsol_sea,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(5),itau_w,wri_fder_ice,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(6),itau_w,wri_evap_ice,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(7),itau_w,wri_evap_sea,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(8),itau_w,wri_rain,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(9),itau_w,wri_snow,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(10),itau_w,wri_rcoa,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(11),itau_w,wri_rriv,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(12),itau_w,wri_calv,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(13),itau_w,wri_tauxx,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(14),itau_w,wri_tauyy,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(15),itau_w,wri_tauzz,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(16),itau_w,wri_tauxx,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(17),itau_w,wri_tauyy,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(18),itau_w,wri_tauzz,iim*(jjm+1),ndexct)
      CALL histsync(nidct)
! pas utile      IF (lafin) CALL histclo(nidct)

      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_calv, wri_tauxx, wri_tauyy,     &
      & wri_tauzz, wri_tauxx, wri_tauyy, wri_tauzz,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
!!$PB
!!$      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),save :: fich ='limit.nc'
  logical, save     :: newlmt = .TRUE.
  logical, save     :: check = .FALSE.
! 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, jour_lu, deja_lu', jour, 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
!
  lmt_sst = 999999999.
  do ii = 1, knon
    lmt_sst(ii) = sst_lu(knindex(ii))
  enddo

  pctsrf_new(:,is_oce) = pct_tmp(:,is_oce)
  pctsrf_new(:,is_sic) = pct_tmp(:,is_sic)

  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),save :: fich ='limit.nc'
  logical,save     :: newlmt = .false.
  logical,save     :: 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 ,save             :: nid, nvarid
  integer, dimension(2),save :: 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
  if (check) 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)
    IF (check) WRITE(*,*)modname,' ouverture fichier ',fich
    if (check) CALL flush(6)
    ierr = NF_OPEN (fich, NF_NOWRITE,nid)
    if (ierr.NE.NF_NOERR) then
      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
!
!!$  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))
    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, qsurf, &
     & 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
!   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
!   qsurf        humidite de l'air au dessus du 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"
#include "indicesol.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, qsurf

! 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
  real, dimension(klon) :: zx_q_0 , d_ts
  real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
  real                  :: bilan_f, fq_fonte
  REAL                  :: subli, fsno
  REAL                  :: qsat_new, q1_new
  real, parameter :: t_grnd = 271.35, t_coup = 273.15
!! PB temporaire en attendant mieux pour le modele de neige
  REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
!
  logical, save         :: check = .true.
  character (len = 20)  :: modname = 'calcul_fluxs'
  logical, save         :: fonte_neige = .false.
  real, save            :: max_eau_sol = 150.0
  character (len = 80) :: abort_message 
  logical,save         :: first = .true.,second=.false.

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

  IF (check) THEN
      WRITE(*,*)' radsol (min, max)' &
         &     , MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
      CALL flush(6)
  ENDIF

  if (size(coastalflow) /= knon .AND. nisurf == is_ter) then
    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
!
!!$  WRITE(*,*)'test calcul_flux, surface ', nisurf
!!PB test
!!$    if (nisurf == is_oce) then
!!$      snow = 0.
!!$      qsol = max_eau_sol
!!$    else
!!$      where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
!!$      where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
!!$!      snow = max(0.0, snow + (precip_snow - evap) * dtime)
!!$      where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
!!$    endif 
!!$    IF (nisurf /= is_ter) qsol = max_eau_sol


! 
! Initialisation
!
  evap = 0.
  fluxsens=0.
  fluxlat=0.
  dflux_s = 0.
  dflux_l = 0.  
!
! 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  
    d_ts(i) = tsurf_new(i) - tsurf(i)
!    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(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))
! Nouvelle valeure de l'humidite au dessus du sol
    qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
    q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
    qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new
!
! en cas de fonte de neige
!
!    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 ,save                  :: 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, knon,dtime,agesno,alb_neig_grid, precip_snow)
  IMPLICIT none
 
  INTEGER :: klon, knon
  INTEGER, PARAMETER :: nvm = 8
  REAL   :: dtime
  REAL, dimension(klon,nvm) :: veget
  REAL, DIMENSION(klon) :: alb_neig_grid, agesno, precip_snow
 
  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, knon
    alb_neig_grid(i) = 0.0
  ENDDO
  DO nv = 1, nvm
    DO i = 1, knon
      as = init(nv)+decay(nv)*EXP(-agesno(i)/5.)
      alb_neig_grid(i) = alb_neig_grid(i) + veget(i,nv)*as
    ENDDO
  ENDDO
!
!! modilation en fonction de l'age de la neige
!
  DO i = 1, knon
    agesno(i)  = (agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
            &             * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/0.3)
    agesno(i) =  MAX(agesno(i),0.0)
  ENDDO
 
  END SUBROUTINE albsno
!
!#########################################################################
!

  SUBROUTINE fonte_neige( 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, &
!IM cf JLD
     & fqcalving,ffonte)

! Routine de traitement de la fonte de la neige dans le cas du traitement
! de sol simplifié
!
! LF 03/2001
! 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         hauteur d'eau contenu dans le sol
!   runoff       runoff en cas de trop plein
!   petAcoef     coeff. A de la resolution de la CL pour t
!   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"
#include "indicesol.inc"
!IM cf JLD
!#include "YOMCST.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(INOUT):: tsurf_new, evap, fluxsens, fluxlat
  real, dimension(klon), intent(INOUT):: dflux_s, dflux_l
! Flux thermique utiliser pour fondre la neige
  real, dimension(klon), intent(INOUT):: ffonte
! Flux d'eau "perdue" par la surface et necessaire pour que limiter la
! hauteur de neige, en kg/m2/s
  real, dimension(klon), intent(INOUT):: fqcalving

! Variables locales
! Masse maximum de neige (kg/m2). Au dessus de ce seuil, la neige
! en exces "s'ecoule" (calving)
  real, parameter :: snow_max=1.
  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
  real, dimension(klon) :: zx_q_0 , d_ts
  real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
  real                  :: bilan_f, fq_fonte
  REAL                  :: subli, fsno
  real, dimension(klon) :: bil_eau_s
  real, parameter :: t_grnd = 271.35, t_coup = 273.15
!! PB temporaire en attendant mieux pour le modele de neige
! REAL, parameter :: chasno = RLMLT/(2.3867E+06*0.15)
  REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
!
  logical, save         :: check = .FALSE.
  character (len = 20)  :: modname = 'fonte_neige'
  logical, save         :: neige_fond = .false.
  real, save            :: max_eau_sol = 150.0
  character (len = 80) :: abort_message 
  logical,save         :: first = .true.,second=.false.

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

! Initialisations
  bil_eau_s(:) = 0.
  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)
  enddo


  WHERE (precip_snow > 0.) snow = snow + (precip_snow * dtime)
  WHERE (evap > 0 ) snow = MAX(0.0, snow - (evap * dtime))
  bil_eau_s = bil_eau_s + (precip_rain - evap) * dtime
!
! Y'a-t-il fonte de neige?
!
  ffonte=0.
  do i = 1, knon
    neige_fond = ((snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
     & .AND. tsurf_new(i) >= RTT)
    if (neige_fond) then
      fq_fonte = MIN( MAX((tsurf_new(i)-RTT )/chasno,0.0),snow(i))
      ffonte(i) = fq_fonte * RLMLT/dtime
      snow(i) = max(0., snow(i) - fq_fonte)
      bil_eau_s(i) = bil_eau_s(i) + fq_fonte 
      tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno  
!IM cf JLD OK     
      IF (nisurf == is_sic .OR. nisurf == is_lic ) tsurf_new(i) = RTT
      d_ts(i) = tsurf_new(i) - tsurf(i)
!      zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(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))
!!$      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)
    endif
!
!   s'il y a une hauteur trop importante de neige, elle s'coule
    fqcalving(i) = max(0., snow(i) - snow_max)/dtime
    snow(i)=min(snow(i),snow_max)
!
    IF (nisurf == is_ter) then
      qsol(i) = qsol(i) + bil_eau_s(i)
      run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0)
      qsol(i) = MIN(qsol(i), max_eau_sol) 
    else
      run_off(i) = run_off(i) + MAX(bil_eau_s(i), 0.0)
    endif
  enddo

  END SUBROUTINE fonte_neige
!
!#########################################################################
!
  END MODULE interface_surf