!!
!! This subroutine is the interface between Fortran 77 main program
!! and SECHIBA.
!! - Input fields are gathered to keep just continental points
!! - call sechiba_main That's SECHIBA process.
!! - Output fields are scattered to complete global fields
!!
!! @call sechiba_main
!! @Version : $Revision: 1.9 $, $Date: 2000/04/05 14:48:08 $
!!
!! @author Marie-Alice Foujols and Jan Polcher
!! 
!! $Header: /home/ssipsl/CVSREP/sechiba/sechiba/intersurf.f90,v 1.9 2000/04/05 14:48:08 ssipsl Exp $
!!
!f90doc MODULE intersurf
MODULE intersurf

  USE IOIPSL

  USE reqdprec
  USE sechiba
  USE constantes
  USE constantes_veg

  IMPLICIT NONE

  PRIVATE
  PUBLIC :: intersurf_main

  INTERFACE intersurf_main
    MODULE PROCEDURE intersurf_main_2d, intersurf_main_1d
  END INTERFACE
  !
  !  Global variables
  !
  INTEGER(long),PARAMETER                           :: max_hist_level = 10
  !
  LOGICAL, SAVE                                     :: l_first_intersurf=.TRUE. !! Initialisation has to be done one time
  !
  INTEGER(long), SAVE                               :: hist_id, rest_id        !! IDs for history and restart files
  INTEGER(long), SAVE                               :: hist_id_stom, rest_id_stom !! Dito for STOMATE
  !
  INTEGER(long), SAVE                               :: itau_offset  !! This offset is used to phase the 
                                                                    !! calendar of the GCM or the driver.
  !
  TYPE(control_type), SAVE                          :: control_flags !! Flags that (de)activate parts of the model
  !
  REAL(stnd), ALLOCATABLE, DIMENSION (:,:), SAVE    :: lalo          !! Geographical coordinates
  INTEGER(long), ALLOCATABLE, DIMENSION (:,:), SAVE :: neighbours    !! neighoring grid points if land
  REAL(stnd), ALLOCATABLE, DIMENSION (:,:), SAVE    :: resolution    !! size in x an y of the grid (m)
  !
CONTAINS
  !
  !f90doc CONTAINS
  ! 
  SUBROUTINE intersurf_main_2d (kjit, iim, jjm, kjpindex, kindex, xrdt, &
     & lrestart_read, lrestart_write, coupling, lon, lat, date0, &
! First level conditions
     & zlev, zlflu,  u, v, qair, temp_air, epot_air, ccanopy, &
! Variables for the implicit coupling
     & cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
! Rain, snow, radiation and surface pressure
     & precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
! Output : Fluxes
     & vevapp, fluxsens, fluxlat, run_off_tot, &
! Surface temperatures and surface properties
     & tsol_rad, temp_sol_new, albedo, emis, z0) 
    
    ! routines called : sechiba_main
    !
    !    
    ! interface description for dummy arguments
    ! input scalar 
    INTEGER(long),INTENT (in)                            :: kjit          !! Time step number
    INTEGER(long),INTENT (in)                            :: iim, jjm      !! Dimension of input fields
    INTEGER(long),INTENT (in)                            :: kjpindex      !! Number of continental points
    REAL(stnd),INTENT (in)                               :: xrdt          !! Time step in seconds
    LOGICAL, INTENT (in)                                 :: lrestart_read !! Logical for _restart_ file to read
    LOGICAL, INTENT (in)                                 :: lrestart_write!! Logical for _restart_ file to write'
    CHARACTER*10, INTENT (in)                            :: coupling      !! Describes the type of coupling  
    REAL(stnd), INTENT (in)                              :: date0         !! Date at which kjit = 0
    ! input fields
    INTEGER(long),DIMENSION (kjpindex), INTENT (in)      :: kindex        !! Index for continental points
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: u             !! Lowest level wind speed
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: v             !! Lowest level wind speed 
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: zlev          !! Height of first layer
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: zlflu         !! Height of between the two first fluxes
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: qair          !! Lowest level specific humidity
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: precip_rain   !! Rain precipitation
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: precip_snow   !! Snow precipitation
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: lwdown        !! Down-welling long-wave flux 
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: swnet         !! Net surface short-wave flux
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: swdown        !! Downwelling surface short-wave flux
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: temp_air      !! Air temperature in Kelvin
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: epot_air      !! Air humidity
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: ccanopy       !! CO2 concentration in the canopy
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: petAcoef      !! Coeficients A from the PBL resolution
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: peqAcoef      !! One for T and another for q
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: petBcoef      !! Coeficients B from the PBL resolution
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: peqBcoef      !! One for T and another for q
    REAL(stnd),DIMENSION (iim,jjm), INTENT(inout)          :: cdrag         !! Cdrag
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: pb            !! Lowest level pressure
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in)             :: lon, lat      !! Geographical coordinates
    ! output fields
    REAL(stnd),DIMENSION (iim,jjm), INTENT(out)            :: z0            !! Surface roughness
    REAL(stnd),DIMENSION (iim,jjm), INTENT(out)            :: run_off_tot   !! Complete runoff
    REAL(stnd),DIMENSION (iim,jjm), INTENT(out)            :: tsol_rad      !!
    REAL(stnd),DIMENSION (iim,jjm), INTENT(out)            :: vevapp        !! Total of evaporation
    REAL(stnd),DIMENSION (iim,jjm), INTENT(out)            :: temp_sol_new  !! New soil temperature
    REAL(stnd),DIMENSION (iim,jjm,2), INTENT(out)          :: albedo        !! Albedo
    REAL(stnd),DIMENSION (iim,jjm), INTENT(out)            :: fluxsens      !! Sensible chaleur flux
    REAL(stnd),DIMENSION (iim,jjm), INTENT(out)            :: fluxlat       !! Latent chaleur flux
    REAL(stnd),DIMENSION (iim,jjm), INTENT(out)            :: emis          !! Emissivity
    ! LOCAL declaration
    ! work arrays to scatter and/or gather information just before/after sechiba_main call's
    ! and to keep output value for next call
    REAL(stnd),DIMENSION (kjpindex)                      :: zu            !! Work array to keep u
    REAL(stnd),DIMENSION (kjpindex)                      :: zv            !! Work array to keep v
    REAL(stnd),DIMENSION (kjpindex)                      :: zzlev         !! Work array to keep zlev
    REAL(stnd),DIMENSION (kjpindex)                      :: zzlflu        !! Work array to keep zlflu
    REAL(stnd),DIMENSION (kjpindex)                      :: zqair         !! Work array to keep qair
    REAL(stnd),DIMENSION (kjpindex)                      :: zprecip_rain  !! Work array to keep precip_rain
    REAL(stnd),DIMENSION (kjpindex)                      :: zprecip_snow  !! Work array to keep precip_snow
    REAL(stnd),DIMENSION (kjpindex)                      :: zlwdown       !! Work array to keep lwdown
    REAL(stnd),DIMENSION (kjpindex)                      :: zswnet        !! Work array to keep swnet
    REAL(stnd),DIMENSION (kjpindex)                      :: zswdown       !! Work array to keep swdown
    REAL(stnd),DIMENSION (kjpindex)                      :: ztemp_air     !! Work array to keep temp_air
    REAL(stnd),DIMENSION (kjpindex)                      :: zepot_air     !! Work array to keep epot_air
    REAL(stnd),DIMENSION (kjpindex)                      :: zccanopy      !! Work array to keep ccanopy
    REAL(stnd),DIMENSION (kjpindex)                      :: zpetAcoef     !! Work array to keep petAcoef
    REAL(stnd),DIMENSION (kjpindex)                      :: zpeqAcoef     !! Work array to keep peqAcoef
    REAL(stnd),DIMENSION (kjpindex)                      :: zpetBcoef     !! Work array to keep petBcoef
    REAL(stnd),DIMENSION (kjpindex)                      :: zpeqBcoef     !! Work array to keep peqVcoef
    REAL(stnd),DIMENSION (kjpindex)                      :: zcdrag        !! Work array to keep cdrag
    REAL(stnd),DIMENSION (kjpindex)                      :: zpb           !! Work array to keep pb
    REAL(stnd),DIMENSION (kjpindex)                      :: zz0           !! Work array to keep z0
    REAL(stnd),DIMENSION (kjpindex)                      :: zrun_off_tot  !! Work array to keep run_off_tot
    REAL(stnd),DIMENSION (kjpindex)                      :: ztsol_rad     !! Work array to keep tsol_rad
    REAL(stnd),DIMENSION (kjpindex)                      :: zvevapp       !! Work array to keep vevapp
    REAL(stnd),DIMENSION (kjpindex)                      :: ztemp_sol_new !! Work array to keep temp_sol_new
    REAL(stnd),DIMENSION (kjpindex,2)                    :: zalbedo       !! Work array to keep albedo
    REAL(stnd),DIMENSION (kjpindex)                      :: zfluxsens     !! Work array to keep fluxsens
    REAL(stnd),DIMENSION (kjpindex)                      :: zfluxlat      !! Work array to keep fluxlat
    REAL(stnd),DIMENSION (kjpindex)                      :: zemis         !! Work array to keep emis

    INTEGER(long)                                        :: i, j, ik
    INTEGER(long)                                        :: itau_sechiba
    LOGICAL :: check = .FALSE.

    IF (l_first_intersurf) THEN
       !
       IF ( check ) WRITE(*,*) 'Initialisation of intersurf'
       !       
       !  Configuration of SSL specific parameters
       !
       CALL intsurf_config(control_flags)
       !
       !  Create the internal coordinate table
       !
       IF ( (.NOT.ALLOCATED(lalo))) THEN
          ALLOCATE(lalo(kjpindex,2))
       ENDIF
       !
       DO ik=1,kjpindex
          
          j = ((kindex(ik)-1)/iim) + 1
          i = (kindex(ik) - (j-1)*iim)

          lalo(ik,1) = lat(i,j)
          lalo(ik,2) = lon(i,j)

       ENDDO
       !-
       !- Compute variable to help describe the grid
       !- once the points are gathered.
       !-
       IF ( (.NOT.ALLOCATED(neighbours))) THEN
          ALLOCATE(neighbours(kjpindex,4))
       ENDIF
       IF ( (.NOT.ALLOCATED(resolution))) THEN
          ALLOCATE(resolution(kjpindex,2))
       ENDIF
       !
       CALL grid_stuff (kjpindex, iim, jjm, lon, lat, kindex, neighbours, resolution)
       !
       CALL intsurf_restart(kjit, iim, jjm, lon, lat, date0, xrdt, control_flags, rest_id, rest_id_stom, itau_offset)
       !
       CALL intsurf_history(iim, jjm, lon, lat, kjit+itau_offset, date0, xrdt, control_flags, hist_id, hist_id_stom)
       !
       IF ( check ) WRITE(*,*) 'End of Initialisation of intersurf'
       !
    ENDIF

    !
    ! 1. gather input fields from kindex array
    !    Warning : I'm not sure this interface with one dimension array is the good one

    DO ik=1, kjpindex
      
       j = ((kindex(ik)-1)/iim) + 1
       i = (kindex(ik) - (j-1)*iim)
       
       zu(ik)           = u(i,j)
       zv(ik)           = v(i,j)
       zzlev(ik)        = zlev(i,j)
       zzlflu(ik)       = zlflu(i,j)
       zqair(ik)        = qair(i,j)
       zprecip_rain(ik) = precip_rain(i,j)
       zprecip_snow(ik) = precip_snow(i,j)
       zlwdown(ik)      = lwdown(i,j)
       zswnet(ik)       = swnet(i,j)
       zswdown(ik)      = swdown(i,j)
       ztemp_air(ik)    = temp_air(i,j)
       zepot_air(ik)    = epot_air(i,j)
       zccanopy(ik)     = ccanopy(i,j)
       zpetAcoef(ik)    = petAcoef(i,j)
       zpeqAcoef(ik)    = peqAcoef(i,j)
       zpetBcoef(ik)    = petBcoef(i,j)
       zpeqBcoef(ik)    = peqBcoef(i,j)
       zcdrag(ik)       = cdrag(i,j)
       zpb(ik)          = pb(i,j)
       !
       ! keep output variables from previous sechiba calls (necessary for stack storage)
       !
       zemis(ik)        = emis(i,j)
       zalbedo(ik,1)    = albedo(i,j,1)
       zalbedo(ik,2)    = albedo(i,j,2)
       zz0(ik)          = z0(i,j)
       
    ENDDO
    
    !
    ! 3. call sechiba for continental points only
    !
    IF ( check ) WRITE(*,*) 'Calling sechiba'
    !
    !  Shift the time step to phase the two models
    !
    itau_sechiba = kjit + itau_offset
    !
    CALL sechiba_main (itau_sechiba, iim*jjm, kjpindex, kindex, xrdt, &
       & lrestart_read, lrestart_write, coupling, control_flags, &
       & lalo, neighbours, resolution, &
! First level conditions
       & zzlev, zzlflu, zu, zv, zqair, ztemp_air, zepot_air, zccanopy, &
! Variables for the implicit coupling
       & zcdrag, zpetAcoef, zpeqAcoef, zpetBcoef, zpeqBcoef, &
! Rain, snow, radiation and surface pressure
       & zprecip_rain ,zprecip_snow,  zlwdown, zswnet, zswdown, zpb, &
! Output : Fluxes
       & zvevapp, zfluxsens, zfluxlat, zrun_off_tot, &
! Surface temperatures and surface properties
       & ztsol_rad, ztemp_sol_new, zalbedo, zemis, zz0, &
! File ids
       & rest_id, hist_id, rest_id_stom, hist_id_stom ) 
    
    !
    IF ( check ) WRITE(*,*) 'out of SECHIBA'
    !
    ! 4. scatter output fields
    !

    DO ik=1, kjpindex
      
    
       j = ((kindex(ik)-1)/iim) + 1
       i = (kindex(ik) - (j-1)*iim)

       z0(i,j)           = zz0(ik)
       run_off_tot(i,j)  = zrun_off_tot(ik)
       tsol_rad(i,j)     = ztsol_rad(ik)
       vevapp(i,j)       = zvevapp(ik)
       temp_sol_new(i,j) = ztemp_sol_new(ik)
       albedo(i,j,1)     = zalbedo(ik,1)
       albedo(i,j,2)     = zalbedo(ik,2)
       fluxsens(i,j)     = zfluxsens(ik)
       fluxlat(i,j)      = zfluxlat(ik)
       emis(i,j)         = zemis(ik)
       cdrag(i,j)        = zcdrag(ik)
       
    ENDDO
    !
    IF ( .NOT. l_first_intersurf) THEN
       !
       CALL histwrite (hist_id, 'evap',     itau_sechiba, vevapp, iim*jjm, kindex)
       CALL histwrite (hist_id, 'temp_sol', itau_sechiba, temp_sol_NEW, iim*jjm, kindex)
       CALL histwrite (hist_id, 'tsol_max', itau_sechiba, temp_sol_NEW, iim*jjm, kindex)
       CALL histwrite (hist_id, 'tsol_min', itau_sechiba, temp_sol_NEW, iim*jjm, kindex)
       CALL histwrite (hist_id, 'fluxsens', itau_sechiba, fluxsens, iim*jjm, kindex)
       CALL histwrite (hist_id, 'swnet',    itau_sechiba, swnet, iim*jjm, kindex)
       CALL histwrite (hist_id, 'alb_vis',  itau_sechiba, albedo(:,:,1), iim*jjm, kindex)
       CALL histwrite (hist_id, 'alb_nir',  itau_sechiba, albedo(:,:,2), iim*jjm, kindex)
       CALL histwrite (hist_id, 'tair',     itau_sechiba, temp_air, iim*jjm, kindex)
       CALL histwrite (hist_id, 'qair',     itau_sechiba, qair, iim*jjm, kindex)
       !
    ENDIF
    !
    l_first_intersurf = .FALSE.
    !
    IF (long_print) WRITE (numout,*) ' intersurf_main done '
    
  END SUBROUTINE intersurf_main_2d
!
  SUBROUTINE intersurf_main_1d (kjit, iim, jjm, kjpindex, kindex, xrdt, &
     & lrestart_read, lrestart_write, coupling, lon, lat, date0, &
! First level conditions
     & zlev, zlflu,  u, v, qair, temp_air, epot_air, ccanopy, &
! Variables for the implicit coupling
     & cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
! Rain, snow, radiation and surface pressure
     & precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
! Output : Fluxes
     & vevapp, fluxsens, fluxlat, run_off_tot, &
! Surface temperatures and surface properties
     & tsol_rad, temp_sol_new, albedo, emis, z0) 
    
    ! routines called : sechiba_main
    !
    !    
    ! interface description for dummy arguments
    ! input scalar 
    INTEGER(long),INTENT (in)                            :: kjit          !! Time step number
    INTEGER(long),INTENT (in)                            :: iim, jjm      !! Dimension of input fields
    INTEGER(long),INTENT (in)                            :: kjpindex      !! Number of continental points
    REAL(stnd),INTENT (in)                               :: xrdt          !! Time step in seconds
    LOGICAL, INTENT (in)                                 :: lrestart_read !! Logical for _restart_ file to read
    LOGICAL, INTENT (in)                                 :: lrestart_write!! Logical for _restart_ file to write'
    CHARACTER*10, INTENT (in)                            :: coupling      !! Describes the type of coupling  
    REAL(stnd), INTENT (in)                              :: date0         !! Date at which kjit = 0
    ! input fields
    INTEGER(long),DIMENSION (kjpindex), INTENT (in)      :: kindex        !! Index for continental points
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: u             !! Lowest level wind speed
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: v             !! Lowest level wind speed 
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: zlev          !! Height of first layer
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: zlflu         !! Height of between the two first fluxes
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: qair          !! Lowest level specific humidity
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: precip_rain   !! Rain precipitation
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: precip_snow   !! Snow precipitation
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: lwdown        !! Down-welling long-wave flux 
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: swnet         !! Net surface short-wave flux
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: swdown        !! Downwelling surface short-wave flux
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: temp_air      !! Air temperature in Kelvin
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: epot_air      !! Air humidity
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: ccanopy       !! CO2 concentration in the canopy
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: petAcoef      !! Coeficients A from the PBL resolution
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: peqAcoef      !! One for T and another for q
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: petBcoef      !! Coeficients B from the PBL resolution
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: peqBcoef      !! One for T and another for q
    REAL(stnd),DIMENSION (iim*jjm), INTENT(inout)          :: cdrag         !! Cdrag
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: pb            !! Lowest level pressure
    REAL(stnd),DIMENSION (iim*jjm), INTENT(in)             :: lon, lat      !! Geographical coordinates
    ! output fields
    REAL(stnd),DIMENSION (iim*jjm), INTENT(out)            :: z0            !! Surface roughness
    REAL(stnd),DIMENSION (iim*jjm), INTENT(out)            :: run_off_tot   !! Complete runoff
    REAL(stnd),DIMENSION (iim*jjm), INTENT(out)            :: tsol_rad      !!
    REAL(stnd),DIMENSION (iim*jjm), INTENT(out)            :: vevapp        !! Total of evaporation
    REAL(stnd),DIMENSION (iim*jjm), INTENT(out)            :: temp_sol_new  !! New soil temperature
    REAL(stnd),DIMENSION (iim*jjm,2), INTENT(out)          :: albedo        !! Albedo
    REAL(stnd),DIMENSION (iim*jjm), INTENT(out)            :: fluxsens      !! Sensible chaleur flux
    REAL(stnd),DIMENSION (iim*jjm), INTENT(out)            :: fluxlat       !! Latent chaleur flux
    REAL(stnd),DIMENSION (iim*jjm), INTENT(out)            :: emis          !! Emissivity
    ! LOCAL declaration
    ! work arrays to scatter and/or gather information just before/after sechiba_main call's
    ! and to keep output value for next call
    REAL(stnd),DIMENSION (kjpindex)                      :: zu            !! Work array to keep u
    REAL(stnd),DIMENSION (kjpindex)                      :: zv            !! Work array to keep v
    REAL(stnd),DIMENSION (kjpindex)                      :: zzlev         !! Work array to keep zlev
    REAL(stnd),DIMENSION (kjpindex)                      :: zzlflu        !! Work array to keep zlflu
    REAL(stnd),DIMENSION (kjpindex)                      :: zqair         !! Work array to keep qair
    REAL(stnd),DIMENSION (kjpindex)                      :: zprecip_rain  !! Work array to keep precip_rain
    REAL(stnd),DIMENSION (kjpindex)                      :: zprecip_snow  !! Work array to keep precip_snow
    REAL(stnd),DIMENSION (kjpindex)                      :: zlwdown       !! Work array to keep lwdown
    REAL(stnd),DIMENSION (kjpindex)                      :: zswnet        !! Work array to keep swnet
    REAL(stnd),DIMENSION (kjpindex)                      :: zswdown       !! Work array to keep swdown
    REAL(stnd),DIMENSION (kjpindex)                      :: ztemp_air     !! Work array to keep temp_air
    REAL(stnd),DIMENSION (kjpindex)                      :: zepot_air     !! Work array to keep epot_air
    REAL(stnd),DIMENSION (kjpindex)                      :: zccanopy      !! Work array to keep ccanopy
    REAL(stnd),DIMENSION (kjpindex)                      :: zpetAcoef     !! Work array to keep petAcoef
    REAL(stnd),DIMENSION (kjpindex)                      :: zpeqAcoef     !! Work array to keep peqAcoef
    REAL(stnd),DIMENSION (kjpindex)                      :: zpetBcoef     !! Work array to keep petBcoef
    REAL(stnd),DIMENSION (kjpindex)                      :: zpeqBcoef     !! Work array to keep peqVcoef
    REAL(stnd),DIMENSION (kjpindex)                      :: zcdrag        !! Work array to keep cdrag
    REAL(stnd),DIMENSION (kjpindex)                      :: zpb           !! Work array to keep pb
    REAL(stnd),DIMENSION (kjpindex)                      :: zz0           !! Work array to keep z0
    REAL(stnd),DIMENSION (kjpindex)                      :: zrun_off_tot  !! Work array to keep run_off_tot
    REAL(stnd),DIMENSION (kjpindex)                      :: ztsol_rad     !! Work array to keep tsol_rad
    REAL(stnd),DIMENSION (kjpindex)                      :: zvevapp       !! Work array to keep vevapp
    REAL(stnd),DIMENSION (kjpindex)                      :: ztemp_sol_new !! Work array to keep temp_sol_new
    REAL(stnd),DIMENSION (kjpindex,2)                    :: zalbedo       !! Work array to keep albedo
    REAL(stnd),DIMENSION (kjpindex)                      :: zfluxsens     !! Work array to keep fluxsens
    REAL(stnd),DIMENSION (kjpindex)                      :: zfluxlat      !! Work array to keep fluxlat
    REAL(stnd),DIMENSION (kjpindex)                      :: zemis         !! Work array to keep emis

    INTEGER(long)                                        :: i, j, ik
    INTEGER(long)                                        :: itau_sechiba
    REAL(stnd), ALLOCATABLE, DIMENSION(:,:)              :: tmp_lon, tmp_lat
    LOGICAL :: check = .FALSE.

    IF (l_first_intersurf) THEN
       !
       IF ( check ) WRITE(*,*) 'Initialisation of intersurf'
       !       
       !  Configuration of SSL specific parameters
       !
       CALL intsurf_config(control_flags)
       !
       !  Create the internal coordinate table
       !
       IF ( (.NOT.ALLOCATED(tmp_lon))) THEN
          ALLOCATE(tmp_lon(iim,jjm))
       ENDIF
       IF ( (.NOT.ALLOCATED(tmp_lat))) THEN
          ALLOCATE(tmp_lat(iim,jjm))
       ENDIF
       !
       IF ( (.NOT.ALLOCATED(lalo))) THEN
          ALLOCATE(lalo(kjpindex,2))
       ENDIF
       !
       DO i=1,iim
          DO j=1,jjm
             ik = (j-1)*iim + i
             tmp_lon(i,j) = lon(ik)
             tmp_lat(i,j) = lat(ik)
          ENDDO
       ENDDO
       !
       lalo(:,1) = lat(:)
       lalo(:,2) = lon(:)
       !-
       !- Compute variable to help describe the grid
       !- once the points are gathered.
       !-
       IF ( (.NOT.ALLOCATED(neighbours))) THEN
          ALLOCATE(neighbours(kjpindex,4))
       ENDIF
       IF ( (.NOT.ALLOCATED(resolution))) THEN
          ALLOCATE(resolution(kjpindex,2))
       ENDIF
       !
       CALL grid_stuff (kjpindex, iim, jjm, tmp_lon, tmp_lat, kindex, neighbours, resolution)
       !
       CALL intsurf_restart(kjit, iim, jjm, tmp_lon, tmp_lat, date0, xrdt, control_flags, rest_id, rest_id_stom, itau_offset)
       !
       CALL intsurf_history(iim, jjm, tmp_lon, tmp_lat, kjit+itau_offset, date0, xrdt, control_flags, hist_id, hist_id_stom)
       !
       IF ( check ) WRITE(*,*) 'End of Initialisation of intersurf'
       !
    ENDIF

    !
    ! 1. gather input fields from kindex array
    !    Warning : I'm not sure this interface with one dimension array is the good one

    DO ik=1, kjpindex
      
       zu(ik)           = u(kindex(ik))
       zv(ik)           = v(kindex(ik))
       zzlev(ik)        = zlev(kindex(ik))
       zzlflu(ik)       = zlflu(kindex(ik))
       zqair(ik)        = qair(kindex(ik))
       zprecip_rain(ik) = precip_rain(kindex(ik))
       zprecip_snow(ik) = precip_snow(kindex(ik))
       zlwdown(ik)      = lwdown(kindex(ik))
       zswnet(ik)       = swnet(kindex(ik))
       zswdown(ik)      = swdown(kindex(ik))
       ztemp_air(ik)    = temp_air(kindex(ik))
       zepot_air(ik)    = epot_air(kindex(ik))
       zccanopy(ik)     = ccanopy(kindex(ik))
       zpetAcoef(ik)    = petAcoef(kindex(ik))
       zpeqAcoef(ik)    = peqAcoef(kindex(ik))
       zpetBcoef(ik)    = petBcoef(kindex(ik))
       zpeqBcoef(ik)    = peqBcoef(kindex(ik))
       zcdrag(ik)       = cdrag(kindex(ik))
       zpb(ik)          = pb(kindex(ik))
       !
       ! keep output variables from previous sechiba calls (necessary for stack storage)
       !
       zemis(ik)        = emis(kindex(ik))
       zalbedo(ik,1)    = albedo(kindex(ik),1)
       zalbedo(ik,2)    = albedo(kindex(ik),2)
       zz0(ik)          = z0(kindex(ik))
       
    ENDDO
    
    !
    ! 3. call sechiba for continental points only
    !
    IF ( check ) WRITE(*,*) 'Calling sechiba'
    !
    !  Shift the time step to phase the two models
    !
    itau_sechiba = kjit + itau_offset
    !
    CALL sechiba_main (itau_sechiba, iim*jjm, kjpindex, kindex, xrdt, &
       & lrestart_read, lrestart_write, coupling, control_flags, &
       & lalo, neighbours, resolution, &
! First level conditions
       & zzlev, zzlflu, zu, zv, zqair, ztemp_air, zepot_air, zccanopy, &
! Variables for the implicit coupling
       & zcdrag, zpetAcoef, zpeqAcoef, zpetBcoef, zpeqBcoef, &
! Rain, snow, radiation and surface pressure
       & zprecip_rain ,zprecip_snow,  zlwdown, zswnet, zswdown, zpb, &
! Output : Fluxes
       & zvevapp, zfluxsens, zfluxlat, zrun_off_tot, &
! Surface temperatures and surface properties
       & ztsol_rad, ztemp_sol_new, zalbedo, zemis, zz0, &
! File ids
       & rest_id, hist_id, rest_id_stom, hist_id_stom ) 
    
    !
    IF ( check ) WRITE(*,*) 'out of SECHIBA'
    !
    ! 4. scatter output fields
    !

    DO ik=1, kjpindex
      
       z0(kindex(ik))           = zz0(ik)
       run_off_tot(kindex(ik))  = zrun_off_tot(ik)
       tsol_rad(kindex(ik))     = ztsol_rad(ik)
       vevapp(kindex(ik))       = zvevapp(ik)
       temp_sol_new(kindex(ik)) = ztemp_sol_new(ik)
       albedo(kindex(ik),1)     = zalbedo(ik,1)
       albedo(kindex(ik),2)     = zalbedo(ik,2)
       fluxsens(kindex(ik))     = zfluxsens(ik)
       fluxlat(kindex(ik))      = zfluxlat(ik)
       emis(kindex(ik))         = zemis(ik)
       cdrag(kindex(ik))        = zcdrag(ik)
       
    ENDDO
    !
    IF ( .NOT. l_first_intersurf) THEN
       !
       CALL histwrite (hist_id, 'evap',     itau_sechiba, vevapp, iim*jjm, kindex)
       CALL histwrite (hist_id, 'temp_sol', itau_sechiba, temp_sol_NEW, iim*jjm, kindex)
       CALL histwrite (hist_id, 'tsol_max', itau_sechiba, temp_sol_NEW, iim*jjm, kindex)
       CALL histwrite (hist_id, 'tsol_min', itau_sechiba, temp_sol_NEW, iim*jjm, kindex)
       CALL histwrite (hist_id, 'fluxsens', itau_sechiba, fluxsens, iim*jjm, kindex)
       CALL histwrite (hist_id, 'swnet',    itau_sechiba, swnet, iim*jjm, kindex)
       CALL histwrite (hist_id, 'alb_vis',  itau_sechiba, albedo(:,1), iim*jjm, kindex)
       CALL histwrite (hist_id, 'alb_nir',  itau_sechiba, albedo(:,2), iim*jjm, kindex)
       CALL histwrite (hist_id, 'tair',     itau_sechiba, temp_air, iim*jjm, kindex)
       CALL histwrite (hist_id, 'qair',     itau_sechiba, qair, iim*jjm, kindex)
       !
    ENDIF
    !
    l_first_intersurf = .FALSE.
    !
    IF (long_print) WRITE (numout,*) ' intersurf_main done '
    
  END SUBROUTINE intersurf_main_1d
!
!
  SUBROUTINE intsurf_config(control_flags)
    !
    !  This subroutine reads all the configuration flags which control the behaviour of the model
    !
    TYPE(control_type), INTENT(out)            :: control_flags !! Flags that (de)activate parts of the model
    !
    !
    !Config Key  = STOMATE_OK_CO2
    !Config Desc = Activate CO2?
    !Config Def  = n
    !Config Help = set to TRUE if photosynthesis is to be activated
    !
    control_flags%ok_co2 = .FALSE.
    CALL getin('STOMATE_OK_CO2', control_flags%ok_co2)
    WRITE(*,*) 'photosynthesis: ', control_flags%ok_co2
    !
    !Config Key  = STOMATE_OK_STOMATE
    !Config Desc = Activate STOMATE?
    !Config Def  = n
    !Config Help = set to TRUE if STOMATE is to be activated
    !
    control_flags%ok_stomate = .FALSE.
    CALL getin('STOMATE_OK_STOMATE',control_flags%ok_stomate)
    WRITE(*,*) 'STOMATE is activated: ',control_flags%ok_stomate
    !
    !Config Key  = STOMATE_OK_DGVM
    !Config Desc = Activate DGVM?
    !Config Def  = n
    !Config Help = set to TRUE if DGVM is to be activated
    !
    control_flags%ok_dgvm = .FALSE.
    CALL getin('STOMATE_OK_DGVM',control_flags%ok_dgvm)
    WRITE(*,*) 'LPJ is activated: ',control_flags%ok_dgvm
    !
    !Config Key  = STOMATE_WATCHOUT
    !Config Desc = STOMATE does minimum service
    !Config Def  = n
    !Config Help = set to TRUE if you want STOMATE to read
    !Config        and write its start files and keep track
    !Config        of longer-term biometeorological variables.
    !Config        This is useful if OK_STOMATE is not set,
    !Config        but if you intend to activate STOMATE later.
    !Config        In that case, this run can serve as a 
    !Config        spinup for longer-term biometeorological
    !Config        variables.
    !
    control_flags%stomate_watchout = .FALSE.
    CALL getin('STOMATE_WATCHOUT',control_flags%stomate_watchout)
    WRITE(*,*) 'STOMATE keeps an eye open: ',control_flags%stomate_watchout
    !
    !
    RETURN
    !
  END SUBROUTINE intsurf_config
  !
  !
  !
  SUBROUTINE intsurf_restart(istp, iim, jjm, lon, lat, date0, dt, control_flags, rest_id, rest_id_stom, itau_offset)
    !
    !  This subroutine initialized the restart file for the land-surface scheme
    !
    INTEGER(long), INTENT(in)                  :: istp      !! Time step of the restart file
    INTEGER(long), INTENT(in)                  :: iim, jjm  !! Size in x and y of the data to be handeled
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: lon, lat  !! Logitude and latitude of the data points
    REAL(stnd)                                 :: date0     !! The date at which itau = 0
    REAL(stnd)                                 :: dt        !! Time step
    INTEGER(long), INTENT(out)                 :: rest_id, rest_id_stom   !! ID of the restart file
    INTEGER(long), INTENT(out)                 :: itau_offset
    !
    TYPE(control_type), INTENT(in)             :: control_flags !! Flags that (de)activate parts of the model
    !
    !  LOCAL
    !
    CHARACTER(LEN=80)          :: restname_in, restname_out, stom_restname_in, stom_restname_out
    REAL(stnd)                 :: dt_rest, date0_rest
    INTEGER(long)              :: itau_dep
    INTEGER(long),PARAMETER    :: llm=1
    REAL(stnd), DIMENSION(llm) :: lev
    LOGICAL                    :: overwrite_time
    REAL(stnd)                 :: in_julian, rest_julian
    REAl(stnd)                 :: un_an, un_jour
    !
    !Config  Key  = SECHIBA_restart_in
    !Config  Desc = Name of restart to READ for initial conditions
    !Config  Def  = NONE
    !Config  Help = This is the name of the file which will be opened
    !Config         to extract the initial values of all prognostic
    !Config         values of the model. This has to be a netCDF file.
    !Config         Not truly COADS compliant. NONE will mean that
    !Config         no restart file is to be expected.
!-
    restname_in = 'NONE'
    CALL getin('SECHIBA_restart_in',restname_in)
    WRITE(*,*) 'INPUT RESTART_FILE', restname_in
    !-
    !Config Key  = SECHIBA_rest_out
    !Config Desc = Name of restart files to be created by SECHIBA
    !Config Def  = sechiba_rest_out.nc
    !Config Help = This variable give the name for
    !Config        the restart files. The restart software within
    !Config        IOIPSL will add .nc if needed.
    !
    restname_out = 'restart_out.nc'
    CALL getin('SECHIBA_rest_out', restname_out)
    !
    !Config Key  = SECHIBA_reset_time
    !Config Desc = Option to overrides the time of the restart
    !Config Def  = n
    !Config Help = This option allows the model to override the time
    !Config        found in the restart file of SECHIBA with the time
    !Config        of the first call. That is the restart time of the GCM.
    !
    overwrite_time = .FALSE.
    CALL getin('SECHIBA_reset_time', overwrite_time)
    !
    CALL ioget_calendar(un_an, un_jour)
    !
    in_julian = itau2date(istp, date0, dt)
    date0_rest = date0
    dt_rest = dt
    !
    CALL restini( restname_in, iim, jjm, lon, lat, llm, lev, &
         &  restname_out, itau_dep, date0_rest, dt_rest, rest_id, overwrite_time)
    !
    !  itau_dep of SECHIBA is phased with the GCM if needed
    !
    rest_julian = itau2date(itau_dep, date0_rest, dt_rest)
    !
    IF ( ABS( in_julian - rest_julian) .GT. dt/un_jour ) THEN
       IF ( overwrite_time ) THEN
          WRITE(*,*) 'The SECHIBA restart is not for the same timestep as the GCM,'
          WRITE(*,*) 'the two are synchronized. The land-surface conditions can not impose'
          WRITE(*,*) 'the chronology of the simulation.'
          WRITE(*,*) 'Time step of the GCM :', istp, 'Julian day : ', in_julian
          WRITE(*,*) 'Time step of SECHIBA :', itau_dep, 'Julian day : ', rest_julian
          itau_offset = itau_dep - istp
       ELSE
          WRITE(*,*) 'IN -> OUT :', istp, '->', itau_dep
          WRITE(*,*) 'IN -> OUT :', in_julian, '->', rest_julian
          WRITE(*,*) 'SECHIBA''s restart file is not consistent with the one of the GCM'
          WRITE(*,*) 'Correct the time axis of the restart file or force the code to change it.'
          STOP
       ENDIF
    ELSE
       itau_offset = 0
    ENDIF
    !
    !=====================================================================
    !- 1.5 Restart file for STOMATE
    !=====================================================================
    IF ( control_flags%ok_stomate .OR. control_flags%stomate_watchout ) THEN 
       !-
       ! STOMATE IS ACTIVATED
       !-
       !Config  Key  = STOMATE_RESTART_FILEIN
       !Config  Desc = Name of restart to READ for initial conditions
       !Config         of STOMATE
       !Config  If   = STOMATE_OK_STOMATE || STOMATE_WATCHOUT
       !Config  Def  = NONE
       !Config  Help = This is the name of the file which will be opened
       !Config         to extract the initial values of all prognostic
       !Config         values of STOMATE.
       !-
       stom_restname_in = 'NONE'
       CALL getin('STOMATE_RESTART_FILEIN',stom_restname_in)
       WRITE(*,*) 'STOMATE INPUT RESTART_FILE', stom_restname_in
       !-
       !Config Key  = STOMATE_RESTART_FILEOUT
       !Config Desc = Name of restart files to be created by STOMATE
       !Config  If   = STOMATE_OK_STOMATE || STOMATE_WATCHOUT
       !Config Def  = stomate_restart.nc
       !Config Help = This is the name of the file which will be opened
       !Config        to write the final values of all prognostic values
       !Config        of STOMATE.
       !-
       stom_restname_out = 'stomate_restart.nc'
       CALL getin('STOMATE_RESTART_FILEOUT', stom_restname_out)
       WRITE(*,*) 'STOMATE OUTPUT RESTART_FILE', stom_restname_out
       !-
       CALL restini (stom_restname_in, iim, jjm, lon, lat, llm, lev, &
            &  stom_restname_out, itau_dep, date0_rest, dt_rest, rest_id_stom, overwrite_time)
       !-
    ENDIF
    !
  END SUBROUTINE intsurf_restart
  
  SUBROUTINE intsurf_history(iim, jjm, lon, lat, istp_old, date0, dt, control_flags, hist_id, hist_id_stom)
    !
    !   
    !  This subroutine initialized the history files for the land-surface scheme
    !
    !
    INTEGER(long), INTENT(in)                  :: iim, jjm  !! Size in x and y of the data to be handeled
    REAL(stnd),DIMENSION (iim,jjm), INTENT(in) :: lon, lat  !! Logitude and latitude of the data points
    INTEGER(long), INTENT(in)                  :: istp_old  !! Time step counter
    REAL(stnd), INTENT(in)                     :: date0     !! Julian day at which istp=0
    REAL(stnd), INTENT(in)                     :: dt        !! Time step of the counter in seconds
    !
    TYPE(control_type), INTENT(in)             :: control_flags !! Flags that (de)activate parts of the model
    !
    INTEGER(long), INTENT(out)                 :: hist_id, hist_id_stom !! History file identification for SECHIBA and STOMATE
    !
    !  LOCAL
    !
    CHARACTER(LEN=80) :: histname,stom_histname       !! Name of history file
    REAL(stnd)        :: dw                           !! frequency of history write (sec.)
    CHARACTER(LEN=30)   :: flux_op                    !! Operations to be performed on fluxes
    CHARACTER(LEN=30)   :: flux_sc                    !! Operations which do not include a scatter
    INTEGER(long)     :: hist_level                   !! history output level (default is 10 => maximum output)
    CHARACTER(LEN=40),DIMENSION(max_hist_level) :: &
         & ave, avecels, avescatter, tmax, fluxop, fluxop_sc, &
         & tmincels, tmaxcels, once                   !! The various operation to be performed
    INTEGER(long)     :: i
    INTEGER(long)     :: hori_id                      !! ID of the default horizontal longitude and latitude map.
    INTEGER(long)     :: vegax_id, solax_id, soltax_id !! ID's for two vertical coordinates
    INTEGER(long)     :: hist_PFTaxis_id
    !
    REAL(stnd),DIMENSION(nvm)   :: veg
    REAL(stnd),DIMENSION(ngrnd) :: sol
    REAL(stnd),DIMENSION(nbsoil):: soltyp
    !
    REAL(stnd)                  :: hist_days_stom     !!- GK time step in days for this history file
    REAL(stnd)                  :: hist_dt_stom       !!- GK time step in seconds for this history file
    INTEGER(long),PARAMETER     :: npft = nvm-1       !!- GK Number of PFTs
    REAL(stnd),DIMENSION(npft)  :: hist_PFTaxis       !!- GK An axis we need for the history files
    !
!=====================================================================
!- 3.0 Setting up the history files
!=====================================================================
!- 3.1 SECHIBA
!=====================================================================
   !Config  Key  = OUTPUT_FILE
   !Config  Desc = Name of file in which the output is going
   !Config         to be written
   !Config  Def  = cabauw_out.nc
   !Config  Help = This file is going to be created by the model
   !Config         and will contain the output from the model.
   !Config         This file is a truly COADS compliant netCDF file.
   !Config         It will be generated by the hist software from
   !Config         the IOIPSL package.
!-
   histname='cabauw_out.nc'
   CALL getin('OUTPUT_FILE', histname)
   WRITE(*,*) 'OUTPUT_FILE', histname
!-
   !Config  Key  = WRITE_STEP
   !Config  Desc = Frequency in seconds at which to WRITE output
   !Config  Def  = 86400.0
   !Config  Help = This variables gives the frequency the output of
   !Config         the model should be written into the netCDF file.
   !Config         It does not affect the frequency at which the
   !Config         operations such as averaging are done.
   !Config         That is IF the coding of the calls to histdef
   !Config         are correct !
!-
   dw = un_jour
   CALL getin('WRITE_STEP', dw)
!-
!
   veg(1:nvm)   = (/ (FLOAT(i),i=1,nvm) /)
   sol(1:ngrnd) = (/ (FLOAT(i),i=1,ngrnd) /)   
   soltyp(1:nbsoil) = (/ (FLOAT(i),i=1,nbsoil) /)
!
!- We need to flux averaging operation as when the data is written
!- from within SECHIBA a scatter is needed. In the driver on the other
!- hand the data is 2D and can be written is it is.
!-
   WRITE(flux_op,'("ave(scatter(X*",F8.1,"))")') un_jour/dt
   WRITE(flux_sc,'("ave(X*",F8.1,")")') un_jour/dt
   WRITE(*,*) flux_op, un_jour/dt, dt, dw
!-
   !Config  Key  = SECHIBA_HISTLEVEL
   !Config  Desc = SECHIBA history output level (0..10)
   !Config  Def  = 5
   !Config  Help = Chooses the list of variables in the history file. 
   !Config         Values between 0: nothing is written; 10: everything is 
   !Config         written are available More details can be found on the web under documentation.
   !Config         web under documentation.
!-
   hist_level = 5
   CALL getin('SECHIBA_HISTLEVEL', hist_level)
!-
   WRITE(*,*) 'SECHIBA history level: ',hist_level
   IF ( (hist_level > max_hist_level).OR.(hist_level < 0) ) THEN
     STOP 'This history level is not allowed'
   ENDIF
!-
!- define operations as a function of history level.
!- Above hist_level, operation='never'
!-
   ave(1:10) = 'ave(X)'
   IF (hist_level < 10) THEN
     ave(hist_level+1:10) = 'never'
   ENDIF
   avecels(1:10) = 'ave(cels(X))'
   IF (hist_level < 10) THEN
     avecels(hist_level+1:10) = 'never'
   ENDIF
   avescatter(1:10) = 'ave(scatter(X))'
   IF (hist_level < 10) THEN
     avescatter(hist_level+1:10) = 'never'
   ENDIF
   tmincels(1:10) = 't_min(cels(X))'
   IF (hist_level < 10) THEN
     tmincels(hist_level+1:10) = 'never'
   ENDIF
   tmaxcels(1:10) = 't_max(cels(X))'
   IF (hist_level < 10) THEN
     tmaxcels(hist_level+1:10) = 'never'
   ENDIF
   tmax(1:10) = 't_max(X)'
   IF (hist_level < 10) THEN
     tmax(hist_level+1:10) = 'never'
   ENDIF
   fluxop(1:10) = flux_op
   IF (hist_level < 10) THEN
     fluxop(hist_level+1:10) = 'never'
   ENDIF
   fluxop_sc(1:10) = flux_sc
   IF (hist_level < 10) THEN
     fluxop_sc(hist_level+1:10) = 'never'
   ENDIF
   once(1:10) = 'once(scatter(X))'
   IF (hist_level < 10) THEN
     once(hist_level+1:10) = 'never'
   ENDIF
!-
   CALL histbeg(histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
   &     istp_old, date0, dt, hori_id, hist_id)
!-
   CALL histvert(hist_id, 'veget', 'Vegetation types', '-', &
   &    nvm,   veg, vegax_id)
   CALL histvert(hist_id, 'solth', 'Soil levels',      'm', &
   &    ngrnd, sol, solax_id)
   CALL histvert(hist_id, 'soiltyp', 'Soil types',      '-', &
   &    nbsoil, soltyp, soltax_id)
!-
!- SECHIBA_HISTLEVEL = 1
!-
   CALL histdef(hist_id, 'evap', 'Evaporation', 'mm/d', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop_sc(1), dt,dw)
   CALL histdef(hist_id, 'temp_sol', 'Surface Temperature', 'C', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avecels(1), dt,dw)
   CALL histdef(hist_id, 'rain', 'Rainfall', 'mm/d',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(1), dt,dw)
   CALL histdef(hist_id, 'snowf', 'Snowfall', 'mm/d',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(1), dt,dw)
   CALL histdef(hist_id, 'netrad', 'Net radiation', 'W/m^2',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(1), dt,dw)
   CALL histdef(hist_id, 'lai', 'Leaf Area Index', '-', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(1), dt,dw)
!-
!- SECHIBA_HISTLEVEL = 2
!-
   CALL histdef(hist_id, 'subli', 'Sublimation', 'mm/d', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(2), dt,dw)
   CALL histdef(hist_id, 'evapnu', 'Bare soil evaporation', 'mm/d', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(2), dt,dw)
   CALL histdef(hist_id, 'runoff', 'Surface runoff', 'mm/d', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, fluxop(2), dt,dw)
   CALL histdef(hist_id, 'tair', 'Air Temperature', 'C',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, ave(2), dt,dw)
   CALL histdef(hist_id, 'qair', 'Air humidity', 'g/g',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, ave(2), dt,dw)
   CALL histdef(hist_id, 'alb_vis', 'Albedo visible', '-', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, ave(2), dt,dw)
   CALL histdef(hist_id, 'alb_nir', 'Albedo near infrared', '-', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, ave(2), dt,dw)
   CALL histdef(hist_id, 'z0', 'Surface roughness', '-',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(2), dt,dw)
   CALL histdef(hist_id, 'transpir', 'Transpiration', 'mm/d', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, fluxop(2), dt,dw)
   CALL histdef(hist_id, 'inter', 'Interception loss', 'mm/d', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, fluxop(2), dt,dw)
!-
!- SECHIBA_HISTLEVEL = 3
!-
   CALL histdef(hist_id, 'tsol_max', 'Maximum Surface Temperature',&
   & 'C', iim,jjm, hori_id, 1,1,1, -99, 32, tmaxcels(3), dt,dw)
   CALL histdef(hist_id, 'tsol_min', 'Minimum Surface Temperature',&
   & 'C', iim,jjm, hori_id, 1,1,1, -99, 32, tmincels(3), dt,dw)
   CALL histdef(hist_id, 'fluxsens', 'Sensible Heat Flux', 'W/m^2',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, ave(3), dt,dw)
   CALL histdef(hist_id, 'snow', 'Snow mass', 'kg/m^2', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(3), dt,dw)
   CALL histdef(hist_id, 'snowage', 'Snow age', '?', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(3), dt,dw)
   CALL histdef(hist_id, 'snowice', 'Snow on ice', 'kg/m^2', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(3), dt,dw)
   CALL histdef(hist_id, 'snowiceage', 'Snow age on ice', 'd', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(3), dt,dw)
   CALL histdef(hist_id, 'vegetfrac', 'Fraction of vegetation', '-', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(3), dt,dw)
!-
!- SECHIBA_HISTLEVEL = 4
!-
   CALL histdef(hist_id, 'gqsb', 'Upper Soil Moisture', '-',  &
   & iim,jjm, hori_id, 1, 1, 1, -99, 32, avescatter(4), dt,dw)
   CALL histdef(hist_id, 'bqsb', 'Lower Soil Moisture', '-',  &
   & iim,jjm, hori_id, 1, 1, 1, -99, 32, avescatter(4), dt,dw)
   CALL histdef(hist_id, 'qsintveg', 'Water on canopy', 'Kg/m^2', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(4), dt,dw)
   CALL histdef(hist_id, 'rstruct', 'Structural resistance', 'mm/d', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(4), dt,dw)
   CALL histdef(hist_id, 'gpp', 'CO2 Assimilation', 'gC/m2/d', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, fluxop(4), dt,dw)
   CALL histdef(hist_id, 'precisol', 'Throughfall', 'mm/d',  &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, fluxop(4), dt,dw)
   CALL histdef(hist_id, 'hdry', 'Dry soil height', 'm',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(4), dt,dw)
!-
!- SECHIBA_HISTLEVEL = 5
!-
   CALL histdef(hist_id, 'swnet', 'Net solar radiation', 'W/m^2',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, ave(5), dt,dw)
   CALL histdef(hist_id, 'temp_pheno', 'Temperature for Pheno', 'K',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(5), dt,dw)
!-
!- SECHIBA_HISTLEVEL = 6
!-
   CALL histdef(hist_id, 'ptn', 'Deep ground temperature', 'K', &
   & iim,jjm, hori_id, ngrnd, 1, ngrnd, solax_id, 32, ave(6),  dt,dw)
!-
!- SECHIBA_HISTLEVEL = 7
!-
!-
!- SECHIBA_HISTLEVEL = 8
!-
   CALL histdef(hist_id, 'beta', 'Beta Function', '-',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw)
   CALL histdef(hist_id, 'qcdrag', 'Diffusion coeff.', '-',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw)
   CALL histdef(hist_id, 'vbeta1', 'Beta for sublimation', '-',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw)
   CALL histdef(hist_id, 'vbeta4', 'Beta for bare soil', '-',  &
   & iim,jjm, hori_id, 1,1,1, -99, 32, avescatter(8), dt,dw)
   CALL histdef(hist_id, 'vbetaco2', 'beta for CO2', 'mm/d', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(8), dt,dw)
   CALL histdef(hist_id, 'dtdt', 'Max value of dT/dt', 'K/s', &
   & iim,jjm, hori_id, 1,1,1, -99, 32, tmax(8), dt,dw)
   CALL histdef(hist_id, 'soiltype', 'Fraction of soil textures', '%', &
   & iim,jjm, hori_id, nbsoil, 1, nbsoil, soltax_id, 32, once(8),  dt,dw)
!-
!- SECHIBA_HISTLEVEL = 10
!-
   CALL histdef(hist_id, 'cimean', 'Stomatal CO2 concentation', 'mmole/m2/s', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(10), dt,dw)
   CALL histdef(hist_id, 'vbeta3', 'Beta for Transpiration', 'mm/d', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(10), dt,dw)
   CALL histdef(hist_id, 'rveget', 'Surface resistance', 'mm/d', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(10), dt,dw)
   CALL histdef(hist_id,'vbeta2','Beta for Interception loss','mm/d', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(10), dt,dw)
   CALL histdef(hist_id, 'qsintmax', 'Maximum Interception capacity', 'Kg/m^2', &
   & iim,jjm, hori_id, nvm, 1, nvm, vegax_id, 32, avescatter(10), dt,dw)




!-
   CALL histend(hist_id)
!=====================================================================
!- 3.2 STOMATE's history file
!=====================================================================
   IF ( control_flags%ok_co2 ) THEN
!-
     ! STOMATE IS ACTIVATED
!-
     !Config  Key  = STOMATE_OUTPUT_FILE
     !Config  Desc = Name of file in which STOMATE's output is going
     !Config         to be written
     !Config  Def  = stomate_history.nc
     !Config  Help = This file is going to be created by the model
     !Config         and will contain the output from the model.
     !Config         This file is a truly COADS compliant netCDF file.
     !Config         It will be generated by the hist software from
     !Config         the IOIPSL package.
!-
     stom_histname='stomate_history.nc'
     CALL getin('STOMATE_OUTPUT_FILE', stom_histname)
     WRITE(*,*) 'STOMATE_OUTPUT_FILE', stom_histname
!-
     !Config  Key  = STOMATE_HIST_DT
     !Config  Desc = STOMATE history time step (d)
     !Config  Def  = 10.
     !Config  Help = Time step of the STOMATE history file
!-
     hist_days_stom = 10.
     CALL getin('STOMATE_HIST_DT', hist_days_stom)
     hist_dt_stom = NINT( hist_days_stom ) * un_jour
     WRITE(*,*) 'output frequency for STOMATE history file (d): ', &
                hist_dt_stom/un_jour
!- initialize
     CALL histbeg(stom_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
     &     istp_old, date0, dt, hori_id, hist_id_stom)
!- define PFT axis
     hist_PFTaxis = (/ ( FLOAT(i), i=1,npft ) /)
!- declare this axis
     CALL histvert (hist_id_stom, 'PFT', 'Plant functional type', &
     & '-', npft, hist_PFTaxis, hist_PFTaxis_id)
!- define STOMATE history file
     CALL stom_define_history (hist_id_stom, npft, iim, jjm, &
     & dt, hist_dt_stom, hori_id, hist_PFTaxis_id)
!- end definition
     CALL histend(hist_id_stom)
!-
   ENDIF

   RETURN

  END SUBROUTINE intsurf_history
  
SUBROUTINE stom_define_history &
& (hist_id_stom, npft, iim, jjm, dt, &
&  hist_dt, hist_hori_id, hist_PFTaxis_id)
!---------------------------------------------------------------------
!- Tell ioipsl which variables are to be written
!- and on which grid they are defined
!---------------------------------------------------------------------

!-
!- Input
!-
!- File id
   INTEGER,INTENT(in) :: hist_id_stom
!- number of PFTs
   INTEGER,INTENT(in) :: npft
!- Domain size
   INTEGER,INTENT(in) :: iim, jjm
!- Time step of STOMATE (seconds)
   REAL,INTENT(in)    :: dt
!- Time step of history file (s)
   REAL,INTENT(in)    :: hist_dt
!- id horizontal grid
   INTEGER,INTENT(in) :: hist_hori_id
!- id of PFT axis
   INTEGER,INTENT(in) :: hist_PFTaxis_id
!-
!- 1 local
!-
!- maximum history level
   INTEGER, PARAMETER  :: max_hist_level = 10
!- output level (between 0 and 10)
!-  ( 0:nothing is written, 10:everything is written)
   INTEGER             :: hist_level
!- Character strings to define operations for histdef
   CHARACTER(LEN=40),DIMENSION(max_hist_level) :: ave
!-
!- Define variables to histdef
!-
   INTEGER :: i,k1,k2,k3,k4,k5
!- Number of variables
   INTEGER,PARAMETER  :: n_v = 68
!- Names of variables
   CHARACTER(LEN=20),DIMENSION(n_v) :: v_n = (/ &
 'SPACE_NAT           ','LITTER_STR_AB_NAT   ','LITTER_STR_AB_AGRI  ',&
 'LITTER_MET_AB_NAT   ','LITTER_MET_AB_AGRI  ','LITTER_STR_BE_NAT   ',&
 'LITTER_STR_BE_AGRI  ','LITTER_MET_BE_NAT   ','LITTER_MET_BE_AGRI  ',&
 'DEADLEAF_COVER      ','CARBON_ACTIVE_NAT   ','CARBON_ACTIVE_AGRI  ',&
 'CARBON_SLOW_NAT     ','CARBON_SLOW_AGRI    ','CARBON_PASSIVE_NAT  ',&
 'CARBON_PASSIVE_AGRI ','T2M_MONTH           ','T2M_WEEK            ',&
 'HET_RESP_NAT        ','HET_RESP_AGRI       ','BLACK_CARBON        ',&
 'FIREFRAC_NAT        ','FIREFRAC_AGRI       ','FIREINDEX_NAT       ',&
 'LITTERHUM           ','LAI                 ','VEGET               ',&
 'VEGET_MAX           ','NPP                 ','GPP                 ',&
 'IND                 ','LEAF_M              ','SAP_M_AB            ',&
 'SAP_M_BE            ','HEART_M_AB          ','HEART_M_BE          ',&
 'ROOT_M              ','FRUIT_M             ','RESERVE_M           ',&
 'LEAF_TURN           ','SAP_AB_TURN         ','ROOT_TURN           ',&
 'FRUIT_TURN          ','LEAF_BM_LITTER      ','SAP_AB_BM_LITTER    ',&
 'SAP_BE_BM_LITTER    ','HEART_AB_BM_LITTER  ','HEART_BE_BM_LITTER  ',&
 'ROOT_BM_LITTER      ','FRUIT_BM_LITTER     ','RESERVE_BM_LITTER   ',&
 'MAINT_RESP          ','GROWTH_RESP         ','AGE                 ',&
 'HEIGHT              ','MOISTRESS           ','VCMAX               ',&
 'LEAF_AGE            ','MORTALITY           ','FIREDEATH           ',&
 'IND_ESTAB           ','LIGHT_DEATH         ','BM_ALLOC_LEAF       ',&
 'BM_ALLOC_SAP_AB     ','BM_ALLOC_SAP_BE     ','BM_ALLOC_ROOT       ',&
 'BM_ALLOC_FRUIT      ','BM_ALLOC_RES        ' /)
!---------------------------------
!- Descriptions of variables
   CHARACTER(LEN=50),DIMENSION(n_v) :: v_d = (/ &
 'fraction of total space that is natural           ',&
 'structural litter above nat. ground               ',&
 'structural litter above agric. ground             ',&
 'metabolic litter above nat. ground                ',&
 'metabolic litter above agric. ground              ',&
 'structural litter below nat. ground               ',&
 'structural litter below agric. ground             ',&
 'metabolic litter below nat. ground                ',&
 'metabolic litter below agric. ground              ',&
 'fraction of soil covered by dead leaves           ',&
 'active soil carbon in nat. ground                 ',&
 'active soil carbon in agric. ground               ',&
 'slow soil carbon in nat. ground                   ',&
 'slow soil carbon in agric. ground                 ',&
 'passive soil carbon in nat. ground                ',&
 'passive soil carbon in agric. ground              ',&
 'Monthly 2 m temperature                           ',&
 'Weekly 2 m temperature                            ',&
 'heterotr. resp. from nat. ground                  ',&
 'heterotr. resp. from agric. ground                ',&
 'black carbon on average total ground              ',&
 'Fire fraction on natural ground                   ',&
 'Fire fraction on agricultural ground              ',&
 'Fire index on natural ground                      ',&
 'litter humidity                                   ',&
 'Leaf Area Index                                   ',&
 'Vegetation fraction                               ',&
 'Maximum vegetation fraction (LAI -> infinity)     ',&
 'Net primary productivity                          ',&
 'Gross primary productivity                        ',&
 'Density of individuals                            ',&
 'Leaf mass                                         ',&
 'Sap mass above ground                             ',&
 'Sap mass below ground                             ',&
 'Heartwood mass above ground                       ',&
 'Heartwood mass below ground                       ',&
 'Root mass                                         ',&
 'Fruit mass                                        ',&
 'Carbohydrate reserve mass                         ',&
 'Leaf turnover                                     ',&
 'Sap turnover above                                ',&
 'Root turnover                                     ',&
 'Fruit turnover                                    ',&
 'Leaf death                                        ',&
 'Sap death above ground                            ',&
 'Sap death below ground                            ',&
 'Heartwood death above ground                      ',&
 'Heartwood death below ground                      ',&
 'Root death                                        ',&
 'Fruit death                                       ',&
 'Carbohydrate reserve death                        ',&
 'Maintenance respiration                           ',&
 'Growth respiration                                ',&
 'age                                               ',&
 'height                                            ',&
 'weekly moisture stress                            ',&
 'Maximum rate of carboxylation                     ',&
 'leaf age                                          ',&
 'Fraction of trees that dies (gap)                 ',&
 'Fraction of plants killed by fire                 ',&
 'Density of newly established saplings             ',&
 'Fraction of plants that dies (light competition)  ',&
 'biomass allocated to leaves                       ',&
 'biomass allocated to sapwood above ground         ',&
 'biomass allocated to sapwood below ground         ',&
 'biomass allocated to roots                        ',&
 'biomass allocated to fruits                       ',&
 'biomass allocated to carbohydrate reserve         ' /)
!- Units of variables
   CHARACTER(LEN=20),DIMENSION(n_v) :: v_u = (/ &
 '-                   ','gC/m^2 (n/a)        ','gC/m^2 (n/a)        ',&
 'gC/m^2 (n/a)        ','gC/m^2 (n/a)        ','gC/m^2 (n/a)        ',&
 'gC/m^2 (n/a)        ','gC/m^2 (n/a)        ','gC/m^2 (n/a)        ',&
 '-                   ','gC/m^2 (n/a)        ','gC/m^2 (n/a)        ',&
 'gC/m^2 (n/a)        ','gC/m^2 (n/a)        ','gC/m^2 (n/a)        ',&
 'gC/m^2 (n/a)        ','K                   ','K                   ',&
 'gC/m^2 tot/day      ','gC/m^2 tot/day      ','gC/m^2 tot          ',&
 '1/day               ','1/day               ','-                   ',&
 '-                   ','-                   ','-                   ',&
 '-                   ','gC/day/(m^2 (n/a))  ','gC/day/(m^2 (n/a))  ',&
 '1/(m^2 (n/a))       ','gC/m^2 (n/a)        ','gC/m^2 (n/a)        ',&
 'gC/m^2 (n/a)        ','gC/m^2 (n/a)        ','gC/m^2 (n/a)        ',&
 'gC/m^2 (n/a)        ','gC/m^2 (n/a)        ','gC/m^2 (n/a)        ',&
 'gC/m^2 (n/a)/day    ','gC/m^2 (n/a)/day    ','gC/m^2 (n/a)/day    ',&
 'gC/m^2 (n/a)/day    ','gC/m^2 (n/a)/day    ','gC/m^2 (n/a)/day    ',&
 'gC/m^2 (n/a)/day    ','gC/m^2 (n/a)/day    ','gC/m^2 (n/a)/day    ',&
 'gC/m^2 (n/a)/day    ','gC/m^2 (n/a)/day    ','gC/m^2 (n/a)/day    ',&
 'gC/m^2 tot/day      ','gC/m^2 tot/day      ','years               ',&
 'm                   ','-                   ','-                   ',&
 'days                ','1/day               ','1/day               ',&
 '1/day               ','1/day               ','gC/m**2 nat/agri/dt ',&
 'gC/m**2 nat/agri/dt ','gC/m**2 nat/agri/dt ','gC/m**2 nat/agri/dt ',&
 'gC/m**2 nat/agri/dt ','gC/m**2 nat/agri/dt ' /)
!- Dimensions of variables
   INTEGER,DIMENSION(n_v) :: v_t = (/ (2,i=1,25), (3,i=26,n_v) /)
!- Level of operations for variables
   INTEGER,DIMENSION(n_v) :: v_l = (/ &
    1,  5,  5,  5,  5,  5,  5,  5,  5,  5, &
    5,  5,  5,  5,  5,  5,  1,  1,  3,  3, &
   10,  1,  1, 10,  5,  1,  1,  1,  1,  1, &
    3,  2,  2,  2,  2,  2,  2,  2,  2,  4, &
    4,  4,  4,  4,  4,  4,  4,  4,  4,  4, &
    4,  2,  2,  7,  7,  3,  6,  6,  6,  6, &
    6,  6,  5,  5,  5,  5,  5,  5 /)
!---------------------------------------------------------------------
!=====================================================================
!- 1 history level
!=====================================================================
!- 1.1 define history level
!=====================================================================
   !Config  Key  = STOMATE_HISTLEVEL
   !Config  Desc = STOMATE history output level (0..10)
   !Config  Def  = 10
   !Config  Help = 0: nothing is written; 10: everything is written
!-
   hist_level = 10
   CALL getin('STOMATE_HISTLEVEL', hist_level)
!-
   WRITE(*,*) 'STOMATE history level: ',hist_level
   IF ( (hist_level > max_hist_level).OR.(hist_level < 0) ) THEN
     STOP 'This history level is not allowed'
   ENDIF
!=====================================================================
!- 1.2 define operations according to output level
!=====================================================================
   ave(1:10) = (/ ('ave(scatter(X))',i=1,hist_level), &
    &             ('never          ',i=hist_level+1,10) /)
!=====================================================================
!- 2 surface fields (2d)
!- 3 PFT: 3rd dimension
!=====================================================================
   k2 = 1; k5 = 32;
   DO i=1,n_v
     IF      (v_t(i) == 2) THEN
       k1 = 1; k3 = 1; k4 = -99;
     ELSE IF (v_t(i) == 3) THEN
       k1 = npft; k3 = npft; k4 = hist_PFTaxis_id;
     ELSE
       CYCLE
     ENDIF
     CALL histdef (hist_id_stom, &
     & TRIM(v_n(i)),TRIM(v_d(i)),TRIM(v_u(i)), iim,jjm, &
     & hist_hori_id, k1,k2,k3,k4,k5, ave(v_l(i)), dt, hist_dt)
   ENDDO
!---------------------------------
END SUBROUTINE stom_define_history

SUBROUTINE grid_stuff (npts, iim, jjm, grid_lon, grid_lat, kindex, &
                         neighbours, resolution)

  !
  ! 0 interface
  !
  
  !
  ! 0.1 input
  !
  
  ! Domain size
  INTEGER, INTENT(in)                                  :: npts
  ! Size of cartesian grid
  INTEGER, INTENT(in)                                  :: iim, jjm
  ! Longitudes on cartesian grid
  REAL, DIMENSION(iim,jjm), INTENT(in)                 :: grid_lon
  ! Latitudes on cartesian grid
  REAL, DIMENSION(iim,jjm), INTENT(in)                 :: grid_lat
  ! Indices of STOMATE points on cartesian grid
  INTEGER, DIMENSION(npts), INTENT(in)                 :: kindex

  !
  ! 0.2 output
  !
  
  ! indices of the 4 neighbours of each grid point (1=N, 2=E, 3=S, 4=W)
  !   a zero or negative index means that this neighbour is not a land point
  INTEGER, DIMENSION(npts,4), INTENT(out)              :: neighbours
  ! resolution at each grid point in m (1=E-W, 2=N-S)
  REAL, DIMENSION(npts,2), INTENT(out)                 :: resolution
  
  !
  ! 0.3 local
  !
  
  ! default resolution (m)
  REAL                                                 :: default_resolution
  ! which STOMATE point corresponds to the given point on the cartesian grid
  INTEGER, DIMENSION(iim,jjm)                          :: correspondance
  ! cosine of the latitude
  REAL                                                 :: coslat
  ! number of points where default resolution is used
  INTEGER                                              :: ndefault_lon, ndefault_lat
  ! Indices
  INTEGER                                              :: i,ip,jp
  
  REAL, PARAMETER                          :: R_Earth = 6378000.
  INTEGER, SAVE                            :: bavard=1
  REAL, PARAMETER                          :: min_stomate = 1.E-8
  REAL, SAVE                               :: pi
  
  ! =========================================================================
  
  pi = 4. * ATAN(1.)
  
  IF ( bavard .GE. 4 ) WRITE(*,*) 'Entering grid_stuff'
  
  ! default resolution
  default_resolution = 250000.
  IF ( bavard .GT. 1 ) WRITE(*,*) 'grid stuff: default resolution (m): ',default_resolution
  
  ndefault_lon = 0
  ndefault_lat = 0
  
  ! initialize output
  neighbours(:,:) = -1
  resolution(:,:) = 0.
  
  correspondance(:,:) = -1
  
  DO i = 1, npts
     
     !
     ! 1 find numbers of the latitude and longitude of each point
     !
     
     ! index of latitude
     jp = INT( (kindex(i)-1) /iim ) + 1
     
     ! index of longitude
     ip = kindex(i) - ( jp-1 ) * iim
     
     correspondance(ip,jp) = kindex(i)
     
  ENDDO
  
  DO i = 1, npts
     
     ! index of latitude
     jp = INT( (kindex(i)-1) /iim ) + 1
     
     ! index of longitude
     ip = kindex(i) - ( jp-1 ) * iim
     
     !
     ! 2 resolution
     !
     
     !
     ! 2.1 longitude
     !
     
     ! prevent infinite resolution at the pole
     coslat = AMAX1( COS( grid_lat(ip,jp) * pi/180. ), 0.001 )
     
     IF ( iim .GT. 1 ) THEN
        
        IF ( ip .EQ. 1 ) THEN
           resolution(i,1) = ABS( grid_lon(ip+1,jp) - grid_lon(ip,jp) ) * &
                pi/180. * R_Earth * coslat
        ELSEIF ( ip .EQ. iim ) THEN
           resolution(i,1) = ABS( grid_lon(ip,jp) - grid_lon(ip-1,jp) ) * &
                pi/180. * R_Earth * coslat
        ELSE
           resolution(i,1) = ABS( grid_lon(ip+1,jp) - grid_lon(ip-1,jp) ) / 2. * &
                pi/180. * R_Earth * coslat
        ENDIF
        
     ELSE
        
        resolution(i,1) = default_resolution
        
        ndefault_lon = ndefault_lon + 1
        
     ENDIF
     
     !
     ! 2.2 latitude
     !
     
     IF ( jjm .GT. 1 ) THEN
        
        IF ( jp .EQ. 1 ) THEN
           resolution(i,2) = ABS( grid_lat(ip,jp) - grid_lat(ip,jp+1) ) * pi/180. * R_Earth
        ELSEIF ( jp .EQ. jjm ) THEN
           resolution(i,2) = ABS( grid_lat(ip,jp-1) - grid_lat(ip,jp) ) * pi/180. * R_Earth
        ELSE
           resolution(i,2) = ABS( grid_lat(ip,jp-1) - grid_lat(ip,jp+1) ) / 2. * pi/180. * R_Earth
        ENDIF
        
     ELSE
        
        resolution(i,2) = default_resolution
        
        ndefault_lat = ndefault_lat + 1
        
     ENDIF
     
     !
     ! 3 find neighbours 
     !
     
     ! North
     IF ( jp .GT. 1 ) THEN
        
        neighbours(i,1) = correspondance(ip,jp-1)
        
     ELSE
        
        ! special treatment for the pole if we are really in a 2d grid
        
        IF ( ( iim .GT. 1 ) .AND. ( jjm .GT. 1 ) ) THEN
           IF ( ABS( ( grid_lat(ip,jp) ) - 90. ) .LT. min_stomate ) THEN
              ! the grid point sits exactly on the pole. The neighbour is situated
              ! at a lower latitude.
              neighbours(i,1) = correspondance( MOD(ip+iim/2-1,iim)+1, jp+1 )
           ELSE
              ! look across the North Pole
              neighbours(i,1) = correspondance( MOD(ip+iim/2-1,iim)+1, jp )
           ENDIF
        ENDIF
        
     ENDIF
     
     ! South
     IF ( jp .LT. jjm ) THEN
        
        neighbours(i,3) = correspondance(ip,jp+1)
        
     ELSE
        
        ! special treatment for the pole if we are really in a 2d grid
        
        IF ( ( iim .GT. 1 ) .AND. ( jjm .GT. 1 ) ) THEN
           IF ( ( ABS( grid_lat(ip,jp) ) - 90. ) .LT. min_stomate ) THEN
              ! the grid point sits exactly on the pole. The neighbour is situated
              ! at a lower latitude.
              neighbours(i,3) = correspondance( MOD(ip+iim/2-1,iim)+1, jp-1 )
           ELSE
              ! look across the South Pole
              neighbours(i,3) = correspondance( MOD(ip+iim/2-1,iim)+1, jp )
           ENDIF
        ENDIF
        
     ENDIF
     
     ! East
     IF ( ip .LT. iim ) THEN
        neighbours(i,2) = correspondance(ip+1,jp)
     ELSEIF ( iim .GT. 1 ) THEN
        neighbours(i,2) = correspondance(1,jp)
     ENDIF
     
     ! West 
     IF ( ip .GT. 1 ) THEN
        neighbours(i,4) = correspondance(ip-1,jp)
     ELSEIF ( iim .GT. 1 ) THEN
        neighbours(i,4) = correspondance(iim,jp)
     ENDIF
     
    ENDDO
    
    IF ( bavard .GT. 1 ) THEN
       WRITE(*,*) '  > Total number of points: ',npts
       WRITE(*,*) '  > Using default zonal resolution at',ndefault_lon,' points.'
       WRITE(*,*) '  > Using default meridional resolution at',ndefault_lat,' points.'
    ENDIF
    
    IF ( bavard .GT. 1 ) WRITE(*,*) 'Leaving grid_stuff'
    
  END SUBROUTINE grid_stuff

END MODULE intersurf
