source: LMDZ4/trunk/libf/phylmd/interface_surf.F90 @ 776

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

!ym  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"
!IM
#include "clesphys.inc"

! run_off      ruissellement total
  REAL, ALLOCATABLE, DIMENSION(:),SAVE    :: run_off, run_off_lic
!$OMP THREADPRIVATE(run_off, run_off_lic)
  real, allocatable, dimension(:),save    :: coastalflow, riverflow
!$OMP THREADPRIVATE(coastalflow, riverflow)
!!$PB
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: tmp_rriv, tmp_rcoa,tmp_rlic
!$OMP THREADPRIVATE(tmp_rriv, tmp_rcoa,tmp_rlic)
!! pour simuler la fonte des glaciers antarctiques
  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: coeff_iceberg
!$OMP THREADPRIVATE(coeff_iceberg)
  real, save                              :: surf_maille 
!$OMP THREADPRIVATE(surf_maille)
  real, save                              :: cte_flux_iceberg = 6.3e7
!$OMP THREADPRIVATE(cte_flux_iceberg)  
  integer, save                           :: num_antarctic = 1
!$OMP THREADPRIVATE(num_antarctic)
  REAL, save                              :: tau_calv
!$OMP THREADPRIVATE(tau_calv)
!!$
  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, &
! -- LOOP
      & windsp, &
! -- LOOP
      & 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, &
      & z0_new, pctsrf_new, agesno,fqcalving,fqfonte,ffonte, run_off_lic_0,&
!IM "slab" ocean
      & flux_o, flux_g, tslab, seaice)


   USE mod_phys_lmdz_para
! 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
!   tslab        temperature slab ocean
!   pctsrf_slab  pourcentages (0-1) des sous-surfaces dans le slab
!   tmp_pctsrf_slab = pctsrf_slab
!   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
! -- LOOP
!   windsp       module du vent a 10m
! -- LOOP
!   rugos        rugosite
!   zmasq        masque terre/ocean
!   rugoro       rugosite orographique
!   run_off_lic_0 runoff glacier du pas de temps precedent
!
! 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

#include "iniprint.h"


! 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
  real, dimension(klon), intent(IN) :: tsurf, p1lay
!IM: "slab" ocean
  real :: amn, amx 
  real, dimension(klon), intent(INOUT) :: tslab
  real, allocatable, dimension(:), save :: tmp_tslab
!$OMP THREADPRIVATE(tmp_tslab)
  real, dimension(klon), intent(OUT) :: flux_o, flux_g
  real, dimension(klon), intent(INOUT)        :: seaice ! glace de mer (kg/m2)
  real, dimension(klon)                       :: siceh  ! hauteur glace de mer (m)
  REAL, DIMENSION(klon), INTENT(INOUT) :: radsol,fder
  
!  real, dimension(klon), intent(IN) :: zmasq
  real, dimension(klon), intent(IN) :: zmasq
  real, dimension(klon), intent(IN) :: taux, tauy, rugos, rugoro
! -- LOOP
  real, dimension(klon), intent(IN) :: windsp
! -- LOOP
  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
  real, dimension(klon), intent(INOUT):: run_off_lic_0

! 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 necessaire pour que limiter la
! hauteur de neige, en kg/m2/s. Et quantite d'eau de fonte de la calotte.
!jld a rajouter   real, dimension(klon), intent(INOUT):: fqcalving, fqfonte
  REAL, DIMENSION(klon), INTENT(INOUT):: fqcalving, fqfonte
!IM: "slab" ocean
  real, dimension(klon) :: new_dif_grnd 
!IM: "slab" ocean - Local
  real, parameter :: t_grnd=271.35
  real, dimension(klon) :: zx_sl
  integer i
  real, allocatable, dimension(:), save :: tmp_flux_o, tmp_flux_g
!$OMP THREADPRIVATE(tmp_flux_o, tmp_flux_g)  
  real, allocatable, dimension(:), save :: tmp_radsol
!$OMP THREADPRIVATE(tmp_radsol)
  real, allocatable, dimension(:,:), save :: tmp_pctsrf_slab
!$OMP THREADPRIVATE(tmp_pctsrf_slab)
  real, allocatable, dimension(:), save :: tmp_seaice
!$OMP THREADPRIVATE(tmp_seaice)
! Local
  character (len = 20),save :: modname = 'interfsurf_hq'
!$OMP THREADPRIVATE(modname)
  character (len = 80) :: abort_message 
  logical, save        :: first_call = .true.
!$OMP THREADPRIVATE(first_call)  
  integer, save        :: error
!$OMP THREADPRIVATE(error)  
  integer              :: ii, index
  logical,save              :: check = .true.
!$OMP THREADPRIVATE(check)  
  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.
  INTEGER,dimension(1) :: iloc
  INTEGER                 :: isize
  real, dimension(klon):: fder_prev
  REAL, dimension(klon) :: bidule
  real, dimension(klon) :: ps_tmp,p1lay_tmp
  INTEGER :: j
!
!IM ?? quelques variables pour netcdf
#include "netcdf.inc"

  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
    call conf_interface(tau_calv)
    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.

!IM cf JLD
  ffonte(1:knon)=0.
  fqcalving(1:knon)=0.
  fqfonte  (1:knon)=0.


  cal = 999999. ; beta = 999999. ; dif_grnd = 999999. ; capsol = 999999.
  alb_new = 999999. ; z0_new = 999999. ; alb_neig = 999999.
  tsurf_new = 999999.
  alblw = 999999.

!IM: "slab" ocean; initialisations
  flux_o = 0.
  flux_g = 0.
!
    if (.not. allocated(tmp_flux_o)) then
      allocate(tmp_flux_o(klon), stat = error)
      DO i=1, knon
       tmp_flux_o(knindex(i))=flux_o(i)
      ENDDO
      if (error /= 0) then
        abort_message='Pb allocation tmp_flux_o'
        call abort_gcm(modname,abort_message,1)
      endif
    endif
    if (.not. allocated(tmp_flux_g)) then  
      allocate(tmp_flux_g(klon), stat = error)
      DO i=1, knon
       tmp_flux_g(knindex(i))=flux_g(i)
      ENDDO
      if (error /= 0) then
        abort_message='Pb allocation tmp_flux_g'
        call abort_gcm(modname,abort_message,1)
      endif
    endif
    if (.not. allocated(tmp_radsol)) then  
      allocate(tmp_radsol(klon), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation tmp_radsol'
        call abort_gcm(modname,abort_message,1)
      endif
    endif
    if (.not. allocated(tmp_pctsrf_slab)) then
      allocate(tmp_pctsrf_slab(klon,nbsrf), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation tmp_pctsrf_slab'
        call abort_gcm(modname,abort_message,1)
      endif
    DO i=1, klon
     tmp_pctsrf_slab(i,1:nbsrf)=pctsrf(i,1:nbsrf)
    ENDDO
    endif
!
    if (.not. allocated(tmp_seaice)) then
      allocate(tmp_seaice(klon), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation tmp_seaice'
        call abort_gcm(modname,abort_message,1)
      endif
    IF(check) THEN
     PRINT*,'allocation tmp_seaice nisurf itime',nisurf, itime
    ENDIF
    endif
!
    if (.not. allocated(tmp_tslab)) then
      allocate(tmp_tslab(klon), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation tmp_tslab'
        call abort_gcm(modname,abort_message,1)
      endif
    endif
    DO i=1, klon
     tmp_tslab(i)=tslab(i)
    ENDDO
!
! 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
!cym      
      run_off=0.0
!cym

!!$PB
      ALLOCATE (tmp_rriv(iim,jj_nb), stat=error)
      if (error /= 0) then
        abort_message='Pb allocation tmp_rriv'
        call abort_gcm(modname,abort_message,1)
      endif
      tmp_rriv=0.
      ALLOCATE (tmp_rcoa(iim,jj_nb), stat=error)
      if (error /= 0) then
        abort_message='Pb allocation tmp_rcoa'
        call abort_gcm(modname,abort_message,1)
      endif
      tmp_rcoa=0.
!ym      ALLOCATE (tmp_rlic(iim,jjm+1), stat=error)
      ALLOCATE (tmp_rlic(iim,jj_nb), stat=error)
      if (error /= 0) then
        abort_message='Pb allocation tmp_rlic'
        call abort_gcm(modname,abort_message,1)
      endif
      tmp_rlic=0.
!!$
    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 adapte !
       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, &
     &   fqcalving,fqfonte,ffonte, run_off_lic_0)


     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
!
#ifdef CPP_VEGET
      p1lay_tmp(1:knon)=p1lay(1:knon)/100.
      ps_tmp(1:knon)=ps(1:knon)/100.
 
       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_tmp, ps_tmp, 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)
#else
      abort_message='Pb de coherence: ok_veget = .t. mais CPP_VEGET = .f.'
      call abort_gcm(modname,abort_message,1)
#endif

    endif    
!
! Remplissage des pourcentages de surface
!
    pctsrf_new(:,nisurf) = pctsrf(:,nisurf)

  else if (nisurf == is_oce) then

    if (check) write(*,*)'ocean, nisurf = ',nisurf,'knon=',knon 
!
! 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.

      iloc = maxloc(fder(1:klon))
      if (check) then
        if (fder(iloc(1))> 0.) then
          WRITE(*,*)'**** Debug fder ****'
          WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))
        endif
      endif
!!$
!!$      where(fder.gt.0.) 
!!$        fder = 0.
!!$      endwhere

      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, &
! -- LOOP
      & windsp, &
! -- LOOP
      & zmasq, &
      & tsurf_new, alb_new, &
      & pctsrf_new)

!IM: "slab" ocean
     else if (ocean == 'slab  ') then
      DO i=1, knon
       tsurf_new(i) = tmp_tslab(knindex(i))
      ENDDO
      pctsrf_new = tmp_pctsrf_slab
!
    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_prev = fder    
    fder = fder_prev + dflux_s + dflux_l

      iloc = maxloc(fder(1:klon))
        if (check.and.fder(iloc(1))> 0.) then
          WRITE(*,*)'**** Debug fder****'
          WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))
          WRITE(*,*)'fder_prev, dflux_s, dflux_l',fder_prev(iloc(1)), &
     &                        dflux_s(iloc(1)), dflux_l(iloc(1))
        endif
!!$
!!$      where(fder.gt.0.) 
!!$        fder = 0.
!!$      endwhere

!IM: flux ocean-atmosphere utile pour le "slab" ocean
    DO i=1, knon
     zx_sl(i) = RLVTT
     if (tsurf_new(i) .LT. RTT) zx_sl(i) = RLSTT
!IM     flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i)
     flux_o(i) = fluxsens(i) + fluxlat(i)
     tmp_flux_o(knindex(i)) = flux_o(i)
     tmp_radsol(knindex(i))=radsol(i)
    ENDDO
!
! 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, &
! -- LOOP
      & windsp, &
! -- LOOP
      & zmasq, &
      & tsurf_new, alb_new, &
      & pctsrf_new)

    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,'knon=',knon 
!
! Surface "glace de mer" appel a l'interface avec l'ocean
!
!
    if (ocean == 'couple') then

      cumul =.false.

      iloc = maxloc(fder(1:klon))
      if (check.and.fder(iloc(1))> 0.) then
        WRITE(*,*)'**** Debug fder ****'
        WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))
      endif
!!$
!!$      where(fder.gt.0.) 
!!$        fder = 0.
!!$      endwhere

      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, &
! -- LOOP
      & windsp, &
! -- LOOP
      & zmasq, &
      & tsurf_new, alb_new, &
      & pctsrf_new)

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

!IM: "slab" ocean
    else if (ocean == 'slab  ') then
!IM ajout sicOBSERVE BEG
   IF ( ok_slab_sicOBS) THEN 
!                              ! lecture conditions limites
      CALL interfoce(itime, dtime, jour, & 
             &  klon, nisurf, knon, knindex, &
             &  debut, &
             &  tsurf_new, pctsrf_new)
!
    tmp_pctsrf_slab=pctsrf_new
    print*,'jour lecture pctsrf_new sic =',jour
!
   ELSE !ok_slab_sicOBS 
     pctsrf_new=tmp_pctsrf_slab
   ENDIF
!IM ajout sicOBSERVE END
!
     DO ii = 1, knon
      tsurf_new(ii) = tsurf(ii)
      IF (pctsrf_new(knindex(ii),nisurf) < EPSFRA) then
       snow(ii) = 0.0
       tsurf_new(ii) = RTT - 1.8
       IF (soil_model) tsoil(ii,:) = RTT -1.8
      ENDIF
     ENDDO

     CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
      
     IF (soil_model) THEN 
      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)
     ELSE 
      dif_grnd = 1.0 / tau_gl
      cal = RCPD * calice
      WHERE (snow > 0.0) cal = RCPD * calsno 
     ENDIF
     tsurf_temp = tsurf_new
     beta = 1.0
!
    ELSE
!                              ! lecture conditions limites
      CALL interfoce(itime, dtime, jour, & 
             &  klon, nisurf, knon, knindex, &
             &  debut, &
             &  tsurf_new, pctsrf_new)

!IM cf LF
      DO ii = 1, knon
       tsurf_new(ii) = tsurf(ii)
       IF (pctsrf_new(knindex(ii),nisurf) < EPSFRA) then
          snow(ii) = 0.0
!IM cf LF/JLD         tsurf(ii) = RTT - 1.8
          tsurf_new(ii) = RTT - 1.8
          IF (soil_model) tsoil(ii,:) = 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 = 1.0 / tau_gl
      ELSE 
         dif_grnd = 1.0 / tau_gl
         cal = RCPD * calice
         WHERE (snow > 0.0) cal = RCPD * calsno 
      ENDIF
!IMbadtsurf_temp = tsurf
      tsurf_temp = tsurf_new
      beta = 1.0
    ENDIF !ocean ==

    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)
!
    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, &
             &   fqcalving,fqfonte,ffonte, run_off_lic_0)

!     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
!
!IM: flux entre l'ocean et la glace de mer pour le "slab" ocean
!
    DO i = 1, knon
!
!IM: faire dependre le coefficient de conduction de la glace de mer
!    de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
!    actuel correspond a 3m de glace de mer, cf. L.Li
!
       IF(1.EQ.0) THEN
        IF(siceh(i).GT.0.) THEN
         new_dif_grnd(i) = dif_grnd(i)*3./siceh(i)
        ELSE
         new_dif_grnd(i) = 0.
        ENDIF
       ENDIF !(1.EQ.0) THEN
!
     IF (cal(i).GT.1.0e-15) THEN
       flux_g(i)=(tsurf_new(i)-t_grnd) &
    &                          * dif_grnd(i) *RCPD/cal(i)
!   &                          * new_dif_grnd(i) *RCPD/cal(i)
     ENDIF
     tmp_flux_g(knindex(i))=flux_g(i)
!
!IM: Attention: ne pas initialiser le tmp_radsol puisque c'est deja fait sur is_oce;
!IM:            tmp_radsol doit etre le flux solaire qui arrive sur l'ocean 
!IM:            et non pas celui qui arrive sur la glace de mer
!
    ENDDO

    fder_prev = fder    
    fder = fder_prev + dflux_s + dflux_l

      iloc = maxloc(fder(1:klon))
      if (check.and.fder(iloc(1))> 0.) then
        WRITE(*,*)'**** Debug fder ****'
        WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))
        WRITE(*,*)'fder_prev, dflux_s, dflux_l',fder_prev(iloc(1)), &
     &                        dflux_s(iloc(1)), dflux_l(iloc(1))
      endif
!!$      where(fder.gt.0.) 
!!$        fder = 0.
!!$      endwhere

!
! 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, &
! -- LOOP
      & windsp, &
! -- LOOP
      & zmasq, &
      & tsurf_new, alb_new, &
      & pctsrf_new)

!IM: "slab" ocean 
    else if (ocean == 'slab  ') then
!
   IF (check) THEN
    amn=MIN(tmp_tslab(1),1000.)
    amx=MAX(tmp_tslab(1),-1000.)
    DO i=2, klon
     amn=MIN(tmp_tslab(i),amn)
     amx=MAX(tmp_tslab(i),amx)
    ENDDO
!
    PRINT*,' debut avant interfoce_slab min max tmp_tslab',amn,amx
   ENDIF !(check) THEN
!
       cumul = .true.
       tslab = tmp_tslab   
       call interfoce(klon, debut, itime, dtime, jour, &
     & tmp_radsol, tmp_flux_o, tmp_flux_g, tmp_pctsrf_slab, &
     & tslab, seaice, pctsrf_new)
!
       tmp_seaice=seaice
       tmp_pctsrf_slab=pctsrf_new
       DO i=1, knon
         tmp_tslab(knindex(i))=tslab(knindex(i))
       ENDDO !i
!

    endif

     
    z0_new = 0.002
    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 

    if (.not. allocated(run_off_lic)) then
      allocate(run_off_lic(klon), stat = error)
      if (error /= 0) then
        abort_message='Pb allocation run_off_lic'
        call abort_gcm(modname,abort_message,1)
      endif
      run_off_lic = 0.
    endif
!
! 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, &
     &   fqcalving,fqfonte,ffonte, run_off_lic_0)

! passage du run-off des glaciers calcule dans fonte_neige au coupleur
    bidule=0.
    bidule(1:knon)= run_off_lic(1:knon)    
    call gath2cpl(bidule, tmp_rlic, klon, knon,iim,jjm,knindex)
!
! 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))
!
!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
!       alb_new(1 : knon)  = 0.6 !IM cf FH/GK 
!       alb_new(1 : knon)  = 0.82
!       alb_new(1 : knon)  = 0.77 !211003 Ksta0.77
!       alb_new(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
!IM: KstaTER0.77 & LMD_ARMIP6    
        alb_new(1 : knon)  = 0.77

!
! 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
!
!#########################################################################
!
#ifdef CPP_VEGET
  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
  USE mod_grid_phy_lmdz
  USE mod_phys_lmdz_para
  IMPLICIT NONE
! 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
 ! 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, save, allocatable,dimension(:) :: ktindex
  REAL, dimension(klon) :: bidule
! Essai cdrag
  real, dimension(klon) :: cdrag
  integer :: jjb,jje,ijb,ije
  INTEGER,SAVE :: offset
  REAL, dimension(klon_glo) :: rlon_g,rlat_g
  INTEGER, SAVE          :: orch_comm
#include "temps.inc"
#include "YOMCST.inc"
#include "iniprint.h"

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


  
! initialisation
  
  if (debut) then
    ALLOCATE(ktindex(klon))
  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

    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)

    enddo
   

    
    CALL Gather(rlon,rlon_g)
    CALL Gather(rlat,rlat_g)

    IF (is_mpi_root) THEN
      index = 1
      do jj = 2, jjm
        do ij = 1, iim
          index = index + 1
          lon_scat(ij,jj) = rlon_g(index)
          lat_scat(ij,jj) = rlat_g(index)
        enddo
      enddo
     lon_scat(:,1) = lon_scat(:,2)
     lat_scat(:,1) = rlat_g(1)
     lon_scat(:,jjm+1) = lon_scat(:,2)
     lat_scat(:,jjm+1) = rlat_g(klon_glo)
   ENDIF
   

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


   CALL Init_neighbours(iim,jjm,knon,neighbours,knindex,pctsrf(:,is_ter))

!
!  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(lunout,*)'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 ---


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

!
! Init Orchidee
!
!  if (pole_nord) then 
!    offset=0
!    ktindex(:)=ktindex(:)+iim-1
!  else
!    offset = klon_begin-1+iim-1
!    ktindex(:)=ktindex(:)+MOD(offset,iim)
!    offset=offset-MOD(offset,iim)
!  endif
  
   
  if (debut) then
    CALL Get_orchidee_communicator(knon,orch_comm)
    IF (knon /=0) THEN
      CALL Init_orchidee_index(iim,knon,orch_comm,knindex,offset,ktindex)
    
      call intersurf_main (itime+itau_phy-1, iim, jjm+1,offset, knon, ktindex, & 
       & orch_comm,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)

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

  endif

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

!IM cf. JP ---
    IF (knon /=0) THEN
    
      call intersurf_main (itime+itau_phy, iim, jjm+1,offset, knon, ktindex, & 
       & orch_comm,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_vrai, ps, &
       & evap, fluxsens, fluxlat, coastalflow, riverflow, &
       & tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
       & lon_scat, lat_scat)

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

    bidule=0.
    bidule(1:knon)=riverflow(1:knon)
    call gath2cpl(bidule, tmp_rriv, klon, knon,iim,jj_nb,knindex)
    bidule=0.
    bidule(1:knon)=coastalflow(1:knon)
    call gath2cpl(bidule, tmp_rcoa, klon, knon,iim,jj_nb,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 Init_orchidee_index(iim,knon,orch_comm,knindex,offset,ktindex)
  USE mod_phys_lmdz_para, ONLY : is_parallel
  
  IMPLICIT NONE
    INTEGER,INTENT(IN)  :: iim
    INTEGER,INTENT(IN)  :: knon
    INTEGER,INTENT(IN)  :: orch_comm
    INTEGER,INTENT(IN)  :: knindex(knon)
    INTEGER,INTENT(OUT) :: offset
    INTEGER,INTENT(OUT) :: ktindex(knon)

#ifdef CPP_PARA
    INCLUDE 'mpif.h'
    INTEGER :: status(MPI_STATUS_SIZE)
#endif
    INTEGER :: MyLastPoint
    INTEGER :: LastPoint
    INTEGER :: mpi_rank
    INTEGER :: mpi_size
    INTEGER :: ierr   
    
    MyLastPoint=klon_begin-1+knindex(knon)+iim-1

    IF (is_parallel) THEN
#ifdef CPP_PARA    
      call MPI_COMM_SIZE(orch_comm,mpi_size,ierr)
      call MPI_COMM_RANK(orch_comm,mpi_rank,ierr)
#endif
    ELSE
      mpi_rank=0
      mpi_size=1
    ENDIF
    
    IF (is_parallel) THEN
      IF (mpi_rank /= 0) then
#ifdef CPP_PARA
        CALL MPI_RECV(LastPoint,1,MPI_INTEGER,mpi_rank-1,1234,orch_comm,status,ierr)
#endif
      ENDIF
        
      IF (mpi_rank /= mpi_size-1) THEN
#ifdef CPP_PARA
        CALL MPI_SEND(MyLastPoint,1,MPI_INTEGER,mpi_rank+1,1234,orch_comm,ierr)  
#endif
      ENDIF
    ENDIF
    
    IF (mpi_rank==0) THEN 
      offset=0
    ELSE
     offset=LastPoint-MOD(LastPoint,iim)
    ENDIF
     
    ktindex(:)=knindex(:)+(klon_begin+iim-1)-offset-1   
    

   END SUBROUTINE  Init_orchidee_index 

  
  SUBROUTINE Get_orchidee_communicator(knon,orch_comm)
  USE mod_phys_lmdz_para, ONLY: mpi_rank, COMM_LMDZ_PHY

  IMPLICIT NONE
#ifdef CPP_PARA
    include 'mpif.h'
#endif    
    INTEGER,INTENT(IN)  :: knon
    INTEGER,INTENT(OUT) :: orch_comm
    
    INTEGER :: color
    INTEGER :: ierr
    
    IF (knon==0) THEN 
      color = 0
    ELSE 
      color = 1
    ENDIF

#ifdef CPP_PARA    
    CALL MPI_COMM_SPLIT(COMM_LMDZ_PHY,color,mpi_rank,orch_comm,ierr)
#endif
    
  END SUBROUTINE Get_orchidee_communicator 
    
    
  SUBROUTINE Init_neighbours(iim,jjm,knon,neighbours,ktindex,pctsrf)
  USE mod_phys_lmdz_para, not_use_mpi_root=>mpi_root
  USE dimphy
  IMPLICIT NONE

#ifdef CPP_PARA
  include 'mpif.h'
#endif
  INTEGER :: iim,jjm
  INTEGER :: knon
  INTEGER :: neighbours(knon,8)
  INTEGER :: ktindex(knon)
  REAL :: pctsrf(klon)
  
  INTEGER :: knon_nb(0:mpi_size-1)
  INTEGER,DIMENSION(0:mpi_size-1) :: displs,sendcount
  INTEGER,ALLOCATABLE :: ktindex_g(:)
  REAL*8  :: pctsrf_g(klon2)
  INTEGER,ALLOCATABLE ::neighbours_g(:,:)
  INTEGER :: knon_g
  REAL*8 :: correspond(iim,jjm+1)
  INTEGER :: i,igrid,jj,ij,iglob,ierr,ireal,index
  integer, dimension(8,3) :: off_ini
  integer, dimension(8)   :: offset  
  INTEGER :: ktindex_p(knon)

  IF (is_sequential) THEN
    knon_nb(:)=knon
  ELSE  

#ifdef CPP_PARA  
    CALL MPI_GATHER(knon,1,MPI_INTEGER,knon_nb,1,MPI_INTEGER,0,COMM_LMDZ_PHY,ierr)
#endif
  
  ENDIF
    
  IF (is_mpi_root) THEN
    knon_g=sum(knon_nb(:))
    ALLOCATE(ktindex_g(knon_g))
    ALLOCATE(neighbours_g(knon_g,8))
    neighbours_g(:,:)=-1
    displs(0)=0
    DO i=1,mpi_size-1
      displs(i)=displs(i-1)+knon_nb(i-1)
    ENDDO  
  ENDIF
  
  ktindex_p(:)=ktindex(:)+klon_begin-1+iim-1

  IF (is_sequential) THEN
    ktindex_g(:)=ktindex_p(:)
  ELSE

#ifdef CPP_PARA  
    CALL MPI_GATHERV(ktindex_p,knon,MPI_INTEGER,ktindex_g,knon_nb,displs,MPI_INTEGER,0,COMM_LMDZ_PHY,ierr)
#endif

  ENDIF
     
  CALL Gather(pctsrf,pctsrf_g)
  
  IF (is_mpi_root) THEN
!  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
!
!
! Attention aux poles
!
    do igrid = 1, knon_g
      index = ktindex_g(igrid)
          jj = int((index - 1)/iim) + 1
          ij = index - (jj - 1) * iim
      correspond(ij,jj) = igrid
    enddo

    do igrid = 1, knon_g
      iglob = ktindex_g(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), klon2))
        if (pctsrf_g(ireal) > EPSFRA) then
          jj = int((index - 1)/iim) + 1
          ij = index - (jj - 1) * iim
            neighbours_g(igrid, i) = correspond(ij, jj)
        endif
      enddo
    enddo

!    DO i=0,phy_size-1
!      displs(i)=displs(i)*8
!      sendcount(i)=knon_nb(i)*8
!    ENDDO
  
  ENDIF
  
  DO i=1,8
    IF (is_sequential) THEN
      neighbours(:,i)=neighbours_g(:,i)
    ELSE
#ifdef CPP_PARA
    CALL MPI_SCATTERV(neighbours_g(:,i),knon_nb,displs,MPI_INTEGER,neighbours(:,i),knon,MPI_INTEGER,0,COMM_LMDZ_PHY,ierr) 
#endif
    ENDIF
  ENDDO
  
  END SUBROUTINE Init_neighbours
#endif
!
!#########################################################################
!
  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, &
! -- LOOP
      & windsp, &
! -- LOOP
      & zmasq, &
      & tsurf_new, alb_new, &
      & pctsrf_new)

   USE ioipsl
   USE mod_phys_lmdz_para, not_use_mpi_root=>mpi_root
   
   USE iophy

#ifdef CPP_PSMILE
   USE oasis
#endif
   USE write_field_phy
   implicit none
#include "indicesol.inc"
#include "YOMCST.inc"
! 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
! -- LOOP
!    windsp       module du vent a 10m
! -- LOOP
!   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
! -- LOOP
   real, dimension(klon), intent(IN) :: windsp
! -- LOOP
  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
! -- LOOP
  INTEGER :: nsrf
! -- LOOP
  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
! -- LOOP
  real, allocatable, dimension(:,:),save :: cpl_windsp
! -- LOOP
  real, allocatable, dimension(:,:),save :: cpl_tauy
  REAL, ALLOCATABLE, DIMENSION(:,:),SAVE :: cpl_rriv, cpl_rcoa, cpl_rlic
!!$
! 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
! -- LOOP
  real, allocatable, dimension(:,:,:),save :: tmp_windsp
! -- LOOP
!!$  real, allocatable, dimension(:,:,:),save :: tmp_tauy, tmp_rriv, tmp_rcoa
  REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: tmp_tauy
! variables a passer au coupleur
!ym  real, dimension(iim, jjm+1) :: wri_sol_ice, wri_sol_sea, wri_nsol_ice 
!ym  real, dimension(iim, jjm+1) :: wri_nsol_sea, wri_fder_ice, wri_evap_ice
!ym  REAL, DIMENSION(iim, jjm+1) :: wri_evap_sea, wri_rcoa, wri_rriv
!ym  REAL, DIMENSION(iim, jjm+1) :: wri_rain, wri_snow, wri_taux, wri_tauy
!ym  REAL, DIMENSION(iim, jjm+1) :: wri_calv
!ym  REAL, DIMENSION(iim, jjm+1) :: wri_tauxx, wri_tauyy, wri_tauzz
!ym  REAL, DIMENSION(iim, jjm+1) :: tmp_lon, tmp_lat

  real, dimension(iim, jj_nb) :: wri_sol_ice, wri_sol_sea, wri_nsol_ice 
  real, dimension(iim, jj_nb) :: wri_nsol_sea, wri_fder_ice, wri_evap_ice
  REAL, DIMENSION(iim, jj_nb) :: wri_evap_sea, wri_rcoa, wri_rriv
  REAL, DIMENSION(iim, jj_nb) :: wri_rain, wri_snow, wri_taux, wri_tauy
  REAL, DIMENSION(iim, jj_nb) :: wri_calv
  REAL, DIMENSION(iim, jj_nb) :: wri_tauxx, wri_tauyy, wri_tauzz
  REAL, DIMENSION(iim, jj_nb) :: tmp_lon, tmp_lat
  REAL, DIMENSION(iim, jj_nb) :: wri_windsp

! 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, jj_nb, 3) :: tamp_srf
  integer, allocatable, dimension(:), save :: tamp_ind
  real, allocatable, dimension(:,:),save :: tamp_zmasq
  real, dimension(iim, jj_nb) :: 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,save :: idayref
!med  integer :: itau_w
  integer,save :: itau_w
! -- LOOP
  integer :: nb_interf_cpl
! -- LOOP

  real :: Up,Down
  integer :: ierr
  integer :: il_time_secs
  real :: tmp_field(klon)
  
#include "param_cou.h"
#include "inc_cpl.h"
#include "temps.inc"
#include "iniprint.h"

#ifdef CPP_PARA
  include 'mpif.h'
  integer :: status(MPI_STATUS_SIZE)
#endif

!
! 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
! -- LOOP
     allocate(cpl_windsp(klon,2), stat = error); sum_error = sum_error + error
! -- LOOP
    allocate(cpl_tauy(klon,2), stat = error); sum_error = sum_error + error
!ym    ALLOCATE(cpl_rriv(iim,jjm+1), stat=error); sum_error = sum_error + error
!ym    ALLOCATE(cpl_rcoa(iim,jjm+1), stat=error); sum_error = sum_error + error
!ym    ALLOCATE(cpl_rlic(iim,jjm+1), stat=error); sum_error = sum_error + error
    ALLOCATE(cpl_rriv(iim,jj_nb), stat=error); sum_error = sum_error + error
    ALLOCATE(cpl_rcoa(iim,jj_nb), stat=error); sum_error = sum_error + error
    ALLOCATE(cpl_rlic(iim,jj_nb), stat=error); sum_error = sum_error + error


!!
!ym    allocate(read_sst(iim, jjm+1), stat = error); sum_error = sum_error + error
!ym    allocate(read_sic(iim, jjm+1), stat = error); sum_error = sum_error + error
!ym    allocate(read_sit(iim, jjm+1), stat = error); sum_error = sum_error + error
!ym    allocate(read_alb_sic(iim, jjm+1), stat = error); sum_error = sum_error + error
    allocate(read_sst(iim, jj_nb), stat = error); sum_error = sum_error + error
    allocate(read_sic(iim, jj_nb), stat = error); sum_error = sum_error + error
    allocate(read_sit(iim, jj_nb), stat = error); sum_error = sum_error + error
    allocate(read_alb_sic(iim, jj_nb), stat = error); sum_error = sum_error + error
    read_sst=0.
    read_sic=0.
    read_sit=0.
    read_alb_sic=0.
    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.; cpl_rlic = 0.
! -- LOOP
     cpl_windsp = 0.
! -- LOOP

    sum_error = 0
    allocate(tamp_ind(klon), stat = error); sum_error = sum_error + error
!ym    allocate(tamp_zmasq(iim, jjm+1), stat = error); sum_error = sum_error + error    
    allocate(tamp_zmasq(iim, jj_nb), stat = error); sum_error = sum_error + error    
    tamp_zmasq=1.
    
    do ig = 1, klon
      tamp_ind(ig) = ig
    enddo
    call gath2cpl(zmasq, tamp_zmasq, klon, klon, iim, jj_nb, tamp_ind)
!
! initialisation couplage
!
    idtime = int(dtime)
#ifdef CPP_COUPLE
#ifdef CPP_PSMILE
   CALL inicma(iim, (jjm+1))
#else
   if (is_parallel) then
      abort_message='coupleur parallele uniquement avec PSMILE'
      call abort_gcm(modname,abort_message,1)
   endif
   call inicma(npas , nexca, idtime,(jjm+1)*iim)
#endif
#endif
!
! initialisation sorties netcdf
!
 !ym  IO de check deconnecte pour le moment en //
    IF (is_sequential) THEN
    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)

    ENDIF    ! monocpu
    
! pour simuler la fonte des glaciers antarctiques
! 
!ym => pour le moment, c'est en commentaire, donc je squizze

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

! Initialisations

! calcul des fluxs a passer
! -- LOOP
  nb_interf_cpl = nb_interf_cpl + 1  
  if (check) write(lunout,*)'passage dans interface_surf.F90 :  ',nb_interf_cpl 
! -- LOOP
  cpl_index = 1
  if (nisurf == is_sic) cpl_index = 2
  if (cumul) then
! -- LOOP
      if (check) write(lunout,*)'passage dans cumul '
      if (check) write(lunout,*)'valeur de cpl_index ', cpl_index 
! -- LOOP  
    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)
! -- LOOP
      IF (cpl_index .EQ. 1) THEN
      cpl_windsp(ig,cpl_index) = cpl_windsp(ig,cpl_index) &
       &                          + windsp(ig)      / FLOAT(nexca) 
      ENDIF
! -- LOOP
    enddo
    IF (cpl_index .EQ. 1) THEN 
        cpl_rriv(:,:) = cpl_rriv(:,:) + tmp_rriv(:,:) / FLOAT(nexca)
        cpl_rcoa(:,:) = cpl_rcoa(:,:) + tmp_rcoa(:,:) / FLOAT(nexca)
        cpl_rlic(:,:) = cpl_rlic(:,:) + tmp_rlic(:,:) / 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
#ifdef CPP_COUPLE
#ifdef CPP_PSMILE
      il_time_secs=(itime-1)*dtime
      CALL fromcpl(il_time_secs, iim, jj_nb,                           &
     &        read_sst, read_sic, read_sit, read_alb_sic)
#else
     if (is_parallel) then
        abort_message='coupleur parallele uniquement avec PSMILE'
        call abort_gcm(modname,abort_message,1)
     endif

      call fromcpl(itime-1,(jjm+1)*iim,                                  &
     &        read_sst, read_sic, read_sit, read_alb_sic)
#endif
#endif
!
! sorties NETCDF des champs recus
!
     if (is_sequential) THEN
       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)
     endif
! pas utile      IF (npas-itime.LT.nexca )CALL histclo(nidcs)

!ym      do j = 1, jjm + 1
       do j = 1, jj_nb
         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,jj_nb, 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,jj_nb,2), stat=error); sum_error = sum_error + error
      allocate(tmp_nsol(iim,jj_nb,2), stat=error); sum_error = sum_error + error
      allocate(tmp_rain(iim,jj_nb,2), stat=error); sum_error = sum_error + error
      allocate(tmp_snow(iim,jj_nb,2), stat=error); sum_error = sum_error + error
      allocate(tmp_evap(iim,jj_nb,2), stat=error); sum_error = sum_error + error
      allocate(tmp_tsol(iim,jj_nb,2), stat=error); sum_error = sum_error + error
      allocate(tmp_fder(iim,jj_nb,2), stat=error); sum_error = sum_error + error
      allocate(tmp_albe(iim,jj_nb,2), stat=error); sum_error = sum_error + error
      allocate(tmp_taux(iim,jj_nb,2), stat=error); sum_error = sum_error + error
      allocate(tmp_tauy(iim,jj_nb,2), stat=error); sum_error = sum_error + error
! -- LOOP
       allocate(tmp_windsp(iim,jj_nb,2), stat=error); sum_error = sum_error + error
! -- LOOP
!!$      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,jj_nb,                  knindex)
    call gath2cpl(cpl_nsol(1,cpl_index), tmp_nsol(1,1,cpl_index), klon, knon,iim,jj_nb,                  knindex)
    call gath2cpl(cpl_rain(1,cpl_index), tmp_rain(1,1,cpl_index), klon, knon,iim,jj_nb,                  knindex)
    call gath2cpl(cpl_snow(1,cpl_index), tmp_snow(1,1,cpl_index), klon, knon,iim,jj_nb,                  knindex)
    call gath2cpl(cpl_evap(1,cpl_index), tmp_evap(1,1,cpl_index), klon, knon,iim,jj_nb,                  knindex)
    call gath2cpl(cpl_tsol(1,cpl_index), tmp_tsol(1,1,cpl_index), klon, knon,iim,jj_nb,                  knindex)
    call gath2cpl(cpl_fder(1,cpl_index), tmp_fder(1,1,cpl_index), klon, knon,iim,jj_nb,                  knindex)
    call gath2cpl(cpl_albe(1,cpl_index), tmp_albe(1,1,cpl_index), klon, knon,iim,jj_nb,                  knindex)
    call gath2cpl(cpl_taux(1,cpl_index), tmp_taux(1,1,cpl_index), klon, knon,iim,jj_nb,                  knindex)
    call gath2cpl(cpl_tauy(1,cpl_index), tmp_tauy(1,1,cpl_index), klon, knon,iim,jj_nb,                  knindex)
! -- LOOP
     call gath2cpl(cpl_windsp(1,cpl_index), tmp_windsp(1,1,cpl_index), klon, knon,iim,jj_nb,            knindex)
! -- LOOP

!
! 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.
! -- LOOP
       wri_windsp = 0.
! -- LOOP      
      call gath2cpl(pctsrf(1,is_oce), tamp_srf(1,1,1), klon, klon, iim, jj_nb, tamp_ind)
      call gath2cpl(pctsrf(1,is_sic), tamp_srf(1,1,2), klon, klon, iim, jj_nb, 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)
! -- LOOP
       wri_windsp = tmp_windsp(:,:,1)
! -- LOOP

!!$PB
      wri_rriv = cpl_rriv(:,:)
      wri_rcoa = cpl_rcoa(:,:)

!ym  !! ATTENTION ICI
      
!ym      DO j = 1, jjm + 1
!ym        wri_calv(:,j) = sum(cpl_rlic(:,j)) / iim
!ym      enddo

!Essai OM+JLD : ca marche !!!! (17 mars 2006)
      tamp_srf(:,:,3)=0.
      CALL gath2cpl( pctsrf(1,is_lic), tamp_srf(1,1,3), klon, klon, iim, jj_nb, tamp_ind)

!YM pour retrouver resultat avant tamp_srf(:,3)=1.
      
      DO j = 1, jj_nb
         wri_calv(:,j) = DOT_PRODUCT (cpl_rlic(1:iim,j), tamp_srf(1:iim,j,3)) / REAL(iim)
      ENDDO

!ym      wri_calv(:,:)=0.
!ym      DO j = 1, jjphy_nb
!ym        wri_calv(:,j) = sum(cpl_rlic(:,j))/iim
!ym      enddo

      IF (is_parallel) THEN
        if (.NOT.is_north_pole) then
#ifdef CPP_PARA
          call MPI_RECV(Up,1,MPI_REAL_LMDZ,mpi_rank-1,1234,COMM_LMDZ_PHY,status,ierr)
          call MPI_SEND(wri_calv(1,1),1,MPI_REAL_LMDZ,mpi_rank-1,1234,COMM_LMDZ_PHY,ierr)
#endif
        endif
        
        if (.NOT.is_south_pole) then
#ifdef CPP_PARA
          call MPI_SEND(wri_calv(1,jj_nb),1,MPI_REAL_LMDZ,mpi_rank+1,1234,COMM_LMDZ_PHY,ierr)  
          call MPI_RECV(down,1,MPI_REAL_LMDZ,mpi_rank+1,1234,COMM_LMDZ_PHY,status,ierr)
#endif
        endif
        
        if (mpi_rank /=0 .and. ii_begin /=1) then
          Up=Up+wri_calv(iim,1)
          wri_calv(:,1)=Up
        endif
      
        if (mpi_rank /=mpi_size-1 .and. ii_end /= iim) then
          Down=Down+wri_calv(1,jj_nb)
          wri_calv(:,jj_nb)=Down         
        endif
      ENDIF
      
      where (tamp_zmasq /= 1.)
        deno =  tamp_srf(:,:,1) + tamp_srf(:,:,2)
        wri_rain = tmp_rain(:,:,1) * tamp_srf(:,:,1) / deno +    &
      &            tmp_rain(:,:,2) * tamp_srf(:,:,2) / deno
        wri_snow = tmp_snow(:,:,1) * tamp_srf(:,:,1) / deno +    &
      &            tmp_snow(:,:,2) * tamp_srf(:,:,2) / deno
        wri_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 coordonnees de la grille
!

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! ATENTION A NE PAS OUBLIER
!!! REMPLACER phy2dyn par grid1DTo2D_mpi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

       CALL Grid1DTo2D_mpi(rlon,tmp_lon)
       CALL Grid1DTo2D_mpi(rlat,tmp_lat)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     IF (is_sequential) THEN
       IF (is_north_pole) tmp_lon(:,1) = tmp_lon(:,2)
       IF (is_south_pole) tmp_lon(:,jjm + 1) = tmp_lon(:,jjm)
     ENDIF
!
! sortie netcdf des champs pour le changement de repere
!
      IF (is_sequential) THEN
        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)
      ENDIF  
!
! calcul 3 coordonnees du vent
!
      CALL atm2geo (iim , jj_nb, 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
!
    IF (is_sequential) THEN
      CALL histwrite(nidct,cl_writ(8),itau_w,wri_sol_ice,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(9),itau_w,wri_sol_sea,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(10),itau_w,wri_nsol_ice,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(11),itau_w,wri_nsol_sea,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(12),itau_w,wri_fder_ice,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(13),itau_w,wri_evap_ice,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(14),itau_w,wri_evap_sea,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(15),itau_w,wri_rain,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(16),itau_w,wri_snow,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(17),itau_w,wri_rcoa,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(18),itau_w,wri_rriv,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(19),itau_w,wri_calv,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(1),itau_w,wri_tauxx,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(2),itau_w,wri_tauyy,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(3),itau_w,wri_tauzz,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(4),itau_w,wri_tauxx,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(5),itau_w,wri_tauyy,iim*(jjm+1),ndexct)
      CALL histwrite(nidct,cl_writ(6),itau_w,wri_tauzz,iim*(jjm+1),ndexct)
! -- LOOP
      CALL histwrite(nidct,cl_writ(7),itau_w,wri_windsp,iim*(jjm+1),ndexct)
! -- LOOP
      CALL histsync(nidct)
    ENDIF
! pas utile      IF (lafin) CALL histclo(nidct)
#ifdef CPP_COUPLE
#ifdef CPP_PSMILE
      il_time_secs=(itime-1)*dtime
      
      CALL intocpl(il_time_secs, iim, jjm+1, 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,                       &
! -- LOOP
      & wri_windsp,lafin)
! -- LOOP
#else
      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,                       &
! -- LOOP
      & wri_windsp,lafin)
! -- LOOP
#endif
#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.; cpl_rlic = 0.
! -- LOOP
      cpl_windsp = 0.
! -- LOOP
!
! 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
! -- LOOP
      deallocate(tmp_windsp, stat=error); sum_error = sum_error + error
! -- LOOP
!!$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,jj_nb, knindex)
  else if (nisurf == is_sic) then
    call cpl2gath(read_sit, tsurf_new, klon, knon,iim,jj_nb, knindex)
    call cpl2gath(read_alb_sic, alb_new, klon, knon,iim,jj_nb, knindex)
  endif
  pctsrf_new(:,nisurf) = pctsrf_sav(:,nisurf)
  
  if (mod(itime, nexca) == -1) then
    tmp_field=0.
    do i = 1, knon
      ig = knindex(i)
      tmp_field(ig) = 1.
    enddo    
!ym    call WriteField_phy('knindex',tmp_field,1)
    
    tmp_field=0.
    do i = 1, knon
      ig = knindex(i)
      tmp_field(ig) = tsurf_new(i)
    enddo    
!ym    call WriteField_phy('tsurf_new',tmp_field,1)
    
    tmp_field=0.
    do i = 1, knon
      ig = knindex(i)
      tmp_field(ig) = alb_new(i)
    enddo    
!ym    call WriteField_phy('alb_new',tmp_field,1)
    
!    tmp_field=0.
!   do i = 1, knon
!      ig = knindex(i)
!      tmp_field(ig) = pctsrf_new(i,nisurf)
 !   enddo    
!ym    call WriteField_phy('pctsrf_new', pctsrf_new(:,nisurf),1)
  endif
!ym  do j=1,jjphy_nb
!ym    do i=1,iim
!ym      print *,phy_rank,'read_sst(',i,',',j,')=',read_sst(i,j)
!ym    enddo
!ym  enddo
  
!ym  do i=1,knon
!ym    print *,phy_rank,'tsurf_new(',i,')=',tsurf_new(i)
!ym  enddo
!  if (lafin) call quitcpl

  END SUBROUTINE interfoce_cpl
!
!#########################################################################
!
   SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour, &
 & radsol, fluxo, fluxg, pctsrf, &
 & tslab, seaice, pctsrf_slab)
!
! Cette routine calcule la temperature d'un slab ocean, la glace de mer 
! et les pourcentages de la maille couverte par l'ocean libre et/ou 
! la glace de mer pour un "slab" ocean de 50m
!
! Conception: Laurent Li
! Re-ecriture + adaptation LMDZ4: I. Musat
!
! input:
!   klon         nombre total de points de grille
!   debut        logical: 1er appel a la physique
!   itap         numero du pas de temps
!   dtime        pas de temps de la physique (en s)
!   ijour        jour dans l'annee en cours
!   radsol       rayonnement net au sol (LW + SW)
!   fluxo        flux turbulent (sensible + latent) sur les mailles oceaniques 
!   fluxg        flux de conduction entre la surface de la glace de mer et l'ocean
!   pctsrf       tableau des pourcentages de surface de chaque maille
! output: 
!   tslab        temperature de l'ocean libre
!   seaice       glace de mer (kg/m2)
!   pctsrf_slab  "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab
!
#include "indicesol.inc"
#include "clesphys.inc"
! Parametres d'entree
  integer, intent(IN) :: klon
  logical, intent(IN) :: debut
  INTEGER, intent(IN) :: itap
  REAL, intent(IN) :: dtime
  INTEGER, intent(IN) :: ijour
  REAL, dimension(klon), intent(IN) :: radsol
  REAL, dimension(klon), intent(IN) :: fluxo
  REAL, dimension(klon), intent(IN) :: fluxg
  real, dimension(klon, nbsrf), intent(IN) :: pctsrf
! Parametres de sortie
  real, dimension(klon), intent(INOUT) :: tslab
  real, dimension(klon), intent(INOUT)        :: seaice ! glace de mer (kg/m2)
  real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_slab
!
! Variables locales :
  REAL :: amn, amx
  INTEGER, save :: lmt_pas, julien, idayvrai
!$OMP THREADPRIVATE(lmt_pas, julien, idayvrai) 
  REAL, parameter :: unjour=86400.
  real, allocatable, dimension(:), save :: tmp_tslab, tmp_seaice
  REAL, allocatable, dimension(:), save :: slab_bils
  REAL, allocatable, dimension(:), save :: lmt_bils
!$OMP THREADPRIVATE(tmp_tslab, tmp_seaice,slab_bils,lmt_bils)
  logical,save              :: check = .false.
!$OMP THREADPRIVATE(check) 

!
  REAL, parameter :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)
  REAL, parameter :: cbing=0.334e+05        ! J/kg
  real, dimension(klon)                 :: siceh !hauteur de la glace de mer (m)
  INTEGER :: i
  integer :: sum_error, error
  REAL :: zz, za, zb
!
  character (len = 80) :: abort_message
  character (len = 20) :: modname = 'interfoce_slab'
!
  julien = MOD(ijour,360)
  sum_error = 0
  IF (debut) THEN
   allocate(slab_bils(klon), stat = error); sum_error = sum_error + error
   allocate(lmt_bils(klon), stat = error); sum_error = sum_error + error
   allocate(tmp_tslab(klon), stat = error); sum_error = sum_error + error
   allocate(tmp_seaice(klon), stat = error); sum_error = sum_error + error
   if (sum_error /= 0) then
    abort_message='Pb allocation var. slab_bils,lmt_bils,tmp_tslab,tmp_seaice'
    call abort_gcm(modname,abort_message,1)
   endif
   tmp_tslab=tslab
   tmp_seaice=seaice
   lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
!
   IF (check) THEN
    PRINT*,'interfoce_slab klon, debut, itap, dtime, ijour, &
 &          lmt_pas ', klon, debut, itap, dtime, ijour, &
 &          lmt_pas
   ENDIF !check
!
   PRINT*, '************************'
   PRINT*, 'SLAB OCEAN est actif, prenez precautions !'
   PRINT*, '************************'
!
! a mettre un slab_bils aussi en force !!!
!
    DO i = 1, klon
     slab_bils(i) = 0.0
    ENDDO
!
   ENDIF !debut
!
   IF (check ) THEN
    amn=MIN(tmp_tslab(1),1000.)
    amx=MAX(tmp_tslab(1),-1000.)
    DO i=2, klon
     amn=MIN(tmp_tslab(i),amn)
     amx=MAX(tmp_tslab(i),amx)
    ENDDO
!
    PRINT*,' debut min max tslab',amn,amx 
!
!!
    PRINT*,' itap,lmt_pas unjour',itap,lmt_pas,unjour
   ENDIF 
!!
!
   pctsrf_slab(1:klon,1:nbsrf) = pctsrf(1:klon,1:nbsrf)
!
! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee
! 
   IF (MOD(itap,lmt_pas) .EQ. 1) THEN !1er pas de temps de la journee
    idayvrai = ijour
    CALL condsurf(julien,idayvrai, lmt_bils)
   ENDIF !(MOD(itap-1,lmt_pas) .EQ. 0) THEN

   DO i = 1, klon
     IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
     &  (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
!
! fabriquer de la glace si congelation atteinte:
!
      IF (tmp_tslab(i).LT.(RTT-1.8)) THEN
       zz =  (RTT-1.8)-tmp_tslab(i)
       tmp_seaice(i) = tmp_seaice(i) + cyang/cbing * zz
       seaice(i) = tmp_seaice(i)
       tmp_tslab(i) = RTT-1.8
      ENDIF
!
! faire fondre de la glace si temperature est superieure a 0:
!
      IF ((tmp_tslab(i).GT.RTT) .AND. (tmp_seaice(i).GT.0.0)) THEN
       zz = cyang/cbing * (tmp_tslab(i)-RTT)
       zz = MIN(zz,tmp_seaice(i))
       tmp_seaice(i) = tmp_seaice(i) - zz
       seaice(i) = tmp_seaice(i)
       tmp_tslab(i) = tmp_tslab(i) - zz*cbing/cyang
      ENDIF
!
! limiter la glace de mer a 10 metres (10000 kg/m2)
!
      IF(tmp_seaice(i).GT.45.) THEN
       tmp_seaice(i) = MIN(tmp_seaice(i),10000.0)
      ELSE
       tmp_seaice(i) = 0. 
      ENDIF
      seaice(i) = tmp_seaice(i)
      siceh(i)=tmp_seaice(i)/1000. !en metres
!
! determiner la nature du sol (glace de mer ou ocean libre):
!
! on fait dependre la fraction de seaice "pctsrf(i,is_sic)" 
! de l'epaisseur de seaice :
! pctsrf(i,is_sic)=1. si l'epaisseur de la glace de mer est >= a 20cm
! et pctsrf(i,is_sic) croit lineairement avec seaice de 0. a 20cm d'epaisseur
!

     IF(.NOT.ok_slab_sicOBS) then
      pctsrf_slab(i,is_sic)=MIN(siceh(i)/0.20, &
     &                      1.-(pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)))
      pctsrf_slab(i,is_oce)=1.0 - &
     &      (pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)+pctsrf_slab(i,is_sic))
     ELSE
      IF (i.EQ.1) print*,'cas ok_slab_sicOBS TRUE : passe sur la modif.'
     ENDIF !(.NOT.ok_slab_sicOBS) then
     ENDIF !pctsrf
   ENDDO
!
! Calculer le bilan du flux de chaleur au sol :
!
   DO i = 1, klon
    za = radsol(i) + fluxo(i)
    zb = fluxg(i)
    IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
   &   (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
     slab_bils(i)=slab_bils(i)+(za*pctsrf_slab(i,is_oce) &
   &             +zb*pctsrf_slab(i,is_sic))/ FLOAT(lmt_pas)
    ENDIF
   ENDDO !klon
!
! calcul tslab 
!
   IF (MOD(itap,lmt_pas).EQ.0) THEN !fin de journee
!
! calcul tslab 
    amn=MIN(tmp_tslab(1),1000.)
    amx=MAX(tmp_tslab(1),-1000.)
    DO i = 1, klon
      IF ((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
     &    (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
       tmp_tslab(i) = tmp_tslab(i) + &
     & (slab_bils(i)-lmt_bils(i)) &
     &                         /cyang*unjour
! on remet l'accumulation a 0
       slab_bils(i) = 0.
      ENDIF !pctsrf
!
      IF (check) THEN
       IF(i.EQ.1) THEN  
        PRINT*,'i tmp_tslab MOD(itap,lmt_pas).EQ.0 sicINTER',i,itap, &
      & tmp_tslab(i), tmp_seaice(i)
       ENDIF
!
       amn=MIN(tmp_tslab(i),amn)
       amx=MAX(tmp_tslab(i),amx)
      ENDIF
    ENDDO !klon
   ENDIF !(MOD(itap,lmt_pas).EQ.0) THEN
!
   IF ( check ) THEN
    PRINT*,'fin min max tslab',amn,amx 
   ENDIF 
!
   tslab = tmp_tslab
   seaice = tmp_seaice
  END SUBROUTINE interfoce_slab
!
!#########################################################################
!
  SUBROUTINE interfoce_lim(itime, dtime, jour, &
     & klon, nisurf, knon, knindex, &
     & debut,  &
     & lmt_sst_p, pctsrf_new_p)
     
     USE mod_grid_phy_lmdz
     USE mod_phys_lmdz_para
#include "indicesol.inc"

! 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) :: klon
  integer, intent(IN) :: nisurf
  integer, intent(IN) :: knon
  integer, dimension(klon_loc), intent(in) :: knindex
  logical, intent(IN) :: debut

! Parametres de sortie
  real, intent(out), dimension(klon_loc) :: lmt_sst_p
  real, intent(out), dimension(klon_loc,nbsrf) :: pctsrf_new_p

!  real, dimension(klon) :: lmt_sst
  real, dimension(klon_glo,nbsrf) :: pctsrf_new

! Variables locales
  integer     :: ii
  INTEGER,save :: lmt_pas     ! frequence de lecture des conditions limites 
                             ! (en pas de physique)
!$OMP THREADPRIVATE(lmt_pas)
  logical,save :: deja_lu    ! pour indiquer que le jour a lire a deja
                             ! lu pour une surface precedente
!$OMP THREADPRIVATE(deja_lu)
  integer,save :: jour_lu 
!$OMP THREADPRIVATE(jour_lu)
  integer      :: ierr
  character (len = 20) :: modname = 'interfoce_lim'
  character (len = 80) :: abort_message
  character (len = 20),save :: fich ='limit.nc'
!$OMP THREADPRIVATE(fich)
  logical, save     :: newlmt = .TRUE.
!$OMP THREADPRIVATE(newlmt)
  logical, save     :: check = .FALSE.
!$OMP THREADPRIVATE(check)
! Champs lus dans le fichier de CL
  real, allocatable , save, dimension(:) :: sst_lu_p
!$OMP THREADPRIVATE(sst_lu_p)

  real, allocatable , save, dimension(:,:) :: pct_tmp_p
!$OMP THREADPRIVATE(pct_tmp_p)
  real, dimension(klon_glo,nbsrf) :: pct_tmp
  real, dimension(klon_glo) :: sst_lu
  real, dimension(klon_glo) :: nat_lu
!
! quelques variables pour netcdf
!
#include "netcdf.inc"
  integer              :: nid, nvarid
  integer, dimension(2) :: start, epais
!
! Fin declaration
!
  
  if (debut .and. .not. allocated(sst_lu_p)) then
    lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
    jour_lu = jour - 1
    allocate(sst_lu_p(klon_loc))
    allocate(pct_tmp_p(klon_loc,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
!
!$OMP MASTER
    if (is_mpi_root) then
    
    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
    epais(1) = klon_glo
    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_glo
        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)
   endif ! is_mpi_root

!$OMP END MASTER
!$OMP BARRIER
  call Scatter(sst_lu,sst_lu_p)
  call Scatter(pct_tmp(:,is_oce),pct_tmp_p(:,is_oce))
  call Scatter(pct_tmp(:,is_sic),pct_tmp_p(:,is_sic))
   deja_lu = .true.
   jour_lu = jour
  endif   

!
! Recopie des variables dans les champs de sortie
!
  
  lmt_sst_p = 999999999.
  
  do ii = 1, knon
    lmt_sst_p(ii) = sst_lu_p(knindex(ii))
  enddo

  do ii=1,klon_loc
    pctsrf_new_p(ii,is_oce)=pct_tmp_p(ii,is_oce)
    pctsrf_new_p(ii,is_sic)=pct_tmp_p(ii,is_sic)
  enddo
  

  END SUBROUTINE interfoce_lim

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

     USE mod_grid_phy_lmdz
     USE mod_phys_lmdz_para
     
! 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 rugosite 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_loc), intent(in) :: knindex
  logical, intent(IN) :: debut

! Parametres de sortie
  real, intent(out), dimension(klon_loc) :: lmt_alb_p
  real, intent(out), dimension(klon_loc) :: lmt_rug_p

!  real,  dimension(klon) :: lmt_alb
!  real,  dimension(klon) :: lmt_rug

! Variables locales
  integer     :: ii
  integer,save :: lmt_pas     ! frequence de lecture des conditions limites 
                             ! (en pas de physique)
!$OMP THREADPRIVATE(lmt_pas) 
  logical,save :: deja_lu_sur! pour indiquer que le jour a lire a deja
                             ! lu pour une surface precedente
!$OMP THREADPRIVATE(deja_lu_sur) 
  integer,save :: jour_lu_sur 
!$OMP THREADPRIVATE(jour_lu_sur) 
  integer      :: ierr
  character (len = 20) :: modname = 'interfsur_lim'
  character (len = 80) :: abort_message
  character (len = 20),save :: fich ='limit.nc'
!$OMP THREADPRIVATE(fich) 
  logical,save     :: newlmt = .false.
!$OMP THREADPRIVATE(newlmt) 
  logical,save     :: check = .false.
!$OMP THREADPRIVATE(check) 
! Champs lus dans le fichier de CL
  real, allocatable , save, dimension(:) :: alb_lu_p, rug_lu_p
!$OMP THREADPRIVATE(alb_lu_p, rug_lu_p) 
  real, dimension(klon_glo) :: alb_lu, rug_lu
!
! quelques variables pour netcdf
!
#include "netcdf.inc"
  integer ,save             :: nid, nvarid
!$OMP THREADPRIVATE(nid, nvarid) 
  integer, dimension(2),save :: start, epais
!$OMP THREADPRIVATE(start, epais) 
!
! Fin declaration
!
  
  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_p(klon_loc))
    allocate(rug_lu_p(klon_loc))
  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

!$OMP MASTER
   if (is_mpi_root) 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
    epais(1) = klon_glo
    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 rugosite !
    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)


  endif  !! is_mpi_root
!$OMP END MASTER

    call Scatter(alb_lu,alb_lu_p)
    call Scatter(rug_lu,rug_lu_p)
    
    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_p(:) = 999999.
  lmt_rug_p(:) = 999999.
  DO ii = 1, knon
    lmt_alb_p(ii) = alb_lu_p(knindex(ii))
    lmt_rug_p(ii) = rug_lu_p(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"
#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, 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 = .false.
!$OMP THREADPRIVATE(check)
  character (len = 20)  :: modname = 'calcul_fluxs'
  logical, save         :: fonte_neige = .false.
!$OMP THREADPRIVATE(fonte_neige)
  real, save            :: max_eau_sol = 150.0
!$OMP THREADPRIVATE(max_eau_sol)
  character (len = 80) :: abort_message 
  logical,save         :: first = .true.,second=.false.
!$OMP THREADPRIVATE(first,second)

  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, jmp1, knindex)
  USE mod_phys_lmdz_para
  implicit none
  
! 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, jmp1
  real, dimension(klon_loc)     :: champ_in
  integer, dimension(klon_loc)  :: knindex
! output
  real, dimension(iim,jmp1)  :: champ_out
! local
  integer                   :: i, ig, j
  real, dimension(klon_loc)     :: tamp

  tamp = 0.
  do i = 1, knon
    ig = knindex(i)
    tamp(ig) = champ_in(i)
  enddo    


  CALL Grid1Dto2D_mpi(tamp,champ_out)

  END SUBROUTINE gath2cpl
!
!#########################################################################
!
  SUBROUTINE cpl2gath(champ_in, champ_out, klon, knon, iim, jmp1, knindex)
  USE mod_phys_lmdz_para
  implicit none
! 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, jmp1
  real, dimension(iim,jmp1)     :: champ_in
  integer, dimension(klon_loc)  :: knindex
! output
  real, dimension(klon_loc)  :: champ_out
! local
  integer                   :: i, ig, j
  real, dimension(klon_loc)     :: tamp
  logical ,save                  :: check = .false.
!$OMP THREADPRIVATE(check)

  CALL Grid2Dto1D_mpi(champ_in,tamp)
  
  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
!$OMP THREADPRIVATE(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, &
     & fqcalving,fqfonte,ffonte,run_off_lic_0)

! Routine de traitement de la fonte de la neige dans le cas du traitement
! de sol simplifie !
! 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
! in/out:
!   run_off_lic_0 run off glacier du pas de temps precedent
!

#include "YOETHF.inc"
!rv#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 "perdu" par la surface et necessaire pour que limiter la
! hauteur de neige, en kg/m2/s. Et flux d'eau de fonte de la calotte.
  REAL, DIMENSION(klon), INTENT(INOUT):: fqcalving, fqfonte
  real, dimension(klon), intent(INOUT):: run_off_lic_0
! 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.
!IM cf JLD/GK
  real, parameter :: snow_max=3000.
  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, snow_evap
  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)
!IM cf JLD/ GKtest
  REAL, parameter :: chaice = 3.334E+05/(2.3867E+06*0.15)
! fin GKtest
!
  logical, save         :: check = .FALSE.
!$OMP THREADPRIVATE(check)
  character (len = 20)  :: modname = 'fonte_neige'
  logical, save         :: neige_fond = .false.
!$OMP THREADPRIVATE(neige_fond)
  real, save            :: max_eau_sol = 150.0
!$OMP THREADPRIVATE(max_eau_sol)
  character (len = 80) :: abort_message 
  logical,save         :: first = .true.,second=.false.
!$OMP THREADPRIVATE(first,second)
  real                 :: coeff_rel
#include "FCTTRE.inc"


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

! Initialisations
  coeff_rel = dtime/(tau_calv * rday)
  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)
  snow_evap = 0.
  WHERE (evap > 0. ) 
    snow_evap = MIN (snow / dtime, evap) 
    snow = snow - snow_evap * dtime
    snow = MAX(0.0, snow)
  end where
 
!  bil_eau_s = bil_eau_s + (precip_rain * dtime) - (evap - snow_evap) * dtime
  bil_eau_s = (precip_rain * dtime) - (evap - snow_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
      fqfonte(i) = fq_fonte/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     
!IM cf JLD/ GKtest fonte aussi pour la glace
      IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN
        fq_fonte = MAX((tsurf_new(i)-RTT )/chaice,0.0)
        ffonte(i) = ffonte(i) + fq_fonte * RLMLT/dtime
        IF ( ok_lic_melt ) THEN
           fqfonte(i) = fqfonte(i) + fq_fonte/dtime
           bil_eau_s(i) = bil_eau_s(i) + fq_fonte
        ENDIF
        tsurf_new(i) = RTT
      ENDIF
      d_ts(i) = tsurf_new(i) - tsurf(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 if (nisurf == is_lic) then
      run_off_lic(i) = (coeff_rel *  fqcalving(i)) + &
 &                        (1. - coeff_rel) * run_off_lic_0(i)
      run_off_lic_0(i) = run_off_lic(i)
      run_off_lic(i) = run_off_lic(i) + fqfonte(i)/dtime
    endif
  enddo

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