!!
!! This module computes surface conditions
!! - albedo
!! - roughness
!! - emissivity
!!
!! @author Marie-Alice Foujols and Jan Polcher
!! @Version : $Revision$, $Date$
!! 
!< $HeadURL$
!< $Date$
!< $Author$
!< $Revision$
!! IPSL (2006)
!!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!! 
MODULE condveg
  !
  USE ioipsl
  !
  ! modules used :
  USE constantes
  USE pft_parameters
  USE interpol_help
  USE parallel

  IMPLICIT NONE

  ! public routines :
  ! condveg_main only
  PRIVATE
  PUBLIC :: condveg_main,condveg_clear 

  !
  ! variables used inside condveg module : declaration and initialisation
  !
  LOGICAL, SAVE                     :: l_first_condveg=.TRUE.           !! To keep first call's trace
  !
  REAL(r_std), ALLOCATABLE, SAVE     :: soilalb_dry(:,:)                 !! albedo for the dry bare soil
  REAL(r_std), ALLOCATABLE, SAVE     :: soilalb_wet(:,:)                 !! albedo for the wet bare soil
  ! Ajout Nathalie pour autre calcul soilalbedo
  REAL(r_std), ALLOCATABLE, SAVE     :: soilalb_moy(:,:)                 !! mean soil albedo.
  ! Ajouts Nathalie - Septembre 2008
  REAL(r_std), ALLOCATABLE, SAVE :: alb_bare(:,:)                         !! Soil Albedo, vis(1) and nir(2)
  REAL(r_std), ALLOCATABLE, SAVE :: alb_veget(:,:)                        !! Vegetation Albedo, vis(1) and nir(2)
  !
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)  :: albedo_snow      !! Snow albedo
  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)  :: albedo_glob      !! Mean albedo
  !

CONTAINS

  !! 
  !! Main routine for *condveg* module
  !! - called only one time for initialisation
  !! - called every time step
  !! - called one more time at last time step for writing _restart_ file
  !!
  !! Algorithm:
  !! - call condveg_init for initialisation
  !! - call condveg_var_init for initialisation done every time step
  !! - call condveg_snow for computing the modification of the albedo induced by snow cover
  !!
  !!  @call condveg_init
  !!  @call condveg_var_init
  !!  @call condveg_snow
  !!
  !!
  SUBROUTINE condveg_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index,&
       & lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
       & zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
       & drysoil_frac, height, deadleaf_cover, emis, albedo, z0, roughheight, soiltype, rest_id, hist_id, hist2_id)

    ! interface description:
    ! input scalar 
    INTEGER(i_std), INTENT(in)                       :: kjit             !! Time step number 
    INTEGER(i_std), INTENT(in)                       :: kjpindex         !! Domain size
    INTEGER(i_std),INTENT (in)                       :: rest_id          !! _Restart_ file identifier
    INTEGER(i_std),INTENT (in)                       :: hist_id          !! _History_ file identifier
    INTEGER(i_std), OPTIONAL, INTENT (in)            :: hist2_id          !! _History_ file 2 identifier
    REAL(r_std), INTENT (in)                         :: dtradia          !! Time step in seconds
    LOGICAL, INTENT(in)                             :: ldrestart_read   !! Logical for restart file to read
    LOGICAL, INTENT(in)                             :: ldrestart_write  !! Logical for restart file to write
    ! input fields
    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index            !! Indeces of the points on the map
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)  :: lalo             !! Geographical coordinates
    INTEGER(i_std),DIMENSION (kjpindex,8), INTENT(in):: neighbours       !! neighoring grid points if land
    REAL(r_std), DIMENSION (kjpindex,2), INTENT(in)  :: resolution       !! size in x an y of the grid (m)
    REAL(r_std), DIMENSION (kjpindex), INTENT(in)    :: contfrac         ! Fraction of land in each grid box.
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget            !! Fraction of vegetation types
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max        !! Fraction of vegetation type
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio    !! Fraction of continental ice, lakes, ...
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: totfrac_nobio    !! total fraction of continental ice+lakes+...
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)    :: zlev             !! Height of first layer
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: snow             !! Snow mass [Kg/m^2]
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: snow_age         !! Snow age
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio    !! Snow mass [Kg/m^2] on ice, lakes, ...
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: snow_nobio_age   !! Snow age on ice, lakes, ...
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: drysoil_frac     !! Fraction of visibly Dry soil(between 0 and 1)
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height           !! Vegetation Height (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: deadleaf_cover   !! Fraction of soil covered by dead leaves
    ! output scalar
    ! output fields
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)   :: emis             !! Emissivity
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo           !! Albedo, vis(1) and nir(2)
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)   :: z0               !! Roughness
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)   :: roughheight      !! Effective height for roughness
    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype       !! fraction of soil types
    ! local
    CHARACTER(LEN=80)                               :: var_name         !! To store variables names for I/O
    !
    IF (l_first_condveg) THEN

        IF (long_print) WRITE (numout,*) ' l_first_condveg : call condveg_init '

        CALL condveg_init (kjit, ldrestart_read, kjpindex, index, veget, &
                           lalo, neighbours, resolution, contfrac, rest_id)
        CALL condveg_var_init (ldrestart_read, kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, &
             & drysoil_frac, zlev, height, deadleaf_cover, emis, albedo, z0, roughheight)
        ! Nathalie - Fevrier 2007 - c'est veget_max qu'il faut passer
        ! Sauf en cas de DGVM .... pour le moment car la somme des maxvegetfrac depasse 1!
        ! Si DGVM alors on passe veget.
        IF (control%ok_dgvm) THEN
           CALL condveg_snow (ldrestart_read, kjpindex, veget, frac_nobio, totfrac_nobio, &
                snow, snow_age, snow_nobio, snow_nobio_age, albedo, albedo_snow, albedo_glob)
        ELSE
           CALL condveg_snow (ldrestart_read, kjpindex, veget_max, frac_nobio, totfrac_nobio, &
                snow, snow_age, snow_nobio, snow_nobio_age, albedo, albedo_snow, albedo_glob)
        ENDIF

        RETURN

    ENDIF

    CALL condveg_var_update (ldrestart_read, kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, &
         & drysoil_frac, zlev, height, deadleaf_cover, emis, albedo, z0, roughheight)
    ! Nathalie - Fevrier 2007 - c'est veget_max qu'il faut passer
    ! Sauf en cas de DGVM .... pour le moment car la somme des maxvegetfrac depasse 1!
    ! Si DGVM alors on passe veget.
    IF (control%ok_dgvm) THEN
       CALL condveg_snow (ldrestart_read, kjpindex, veget, frac_nobio, totfrac_nobio, &
            snow, snow_age, snow_nobio, snow_nobio_age, albedo, albedo_snow, albedo_glob)
    ELSE 
       CALL condveg_snow (ldrestart_read, kjpindex, veget_max, frac_nobio, totfrac_nobio, &
            snow, snow_age, snow_nobio, snow_nobio_age, albedo, albedo_snow, albedo_glob)
    ENDIF

    IF (ldrestart_write) THEN
       !
       var_name = 'soilalbedo_dry'
       CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_dry, 'scatter',  nbp_glo, index_g)
       !
       var_name = 'soilalbedo_wet'
       CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_wet, 'scatter',  nbp_glo, index_g)
       !
       var_name = 'soilalbedo_moy'
       CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_moy, 'scatter',  nbp_glo, index_g)
       !
       RETURN
       !
    ENDIF

    IF ( almaoutput ) THEN
       CALL histwrite(hist_id, 'Albedo', kjit, albedo_glob, kjpindex, index)
       CALL histwrite(hist_id, 'SAlbedo', kjit, albedo_snow, kjpindex, index)
       IF ( hist2_id > 0 ) THEN
          CALL histwrite(hist2_id, 'Albedo', kjit, albedo_glob, kjpindex, index)
          CALL histwrite(hist2_id, 'SAlbedo', kjit, albedo_snow, kjpindex, index)
       ENDIF
    ELSE
       CALL histwrite(hist_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
       CALL histwrite(hist_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
       CALL histwrite(hist_id, 'vegalb_vis',  kjit, alb_veget(:,1), kjpindex, index)
       CALL histwrite(hist_id, 'vegalb_nir',  kjit, alb_veget(:,2), kjpindex, index)
       IF ( hist2_id > 0 ) THEN
          CALL histwrite(hist2_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
          CALL histwrite(hist2_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
          CALL histwrite(hist2_id, 'vegalb_vis',  kjit, alb_veget(:,1), kjpindex, index)
          CALL histwrite(hist2_id, 'vegalb_nir',  kjit, alb_veget(:,2), kjpindex, index)
       ENDIF
    ENDIF

    IF (long_print) WRITE (numout,*)' condveg_main done '

  END SUBROUTINE condveg_main

  !! Algorithm:
  !! - dynamic allocation for local array 
  !!
  SUBROUTINE condveg_init  (kjit, ldrestart_read, kjpindex, index, veget, &
                            lalo, neighbours, resolution, contfrac, rest_id)

    ! interface description
    ! input scalar 
    INTEGER(i_std), INTENT(in)                       :: kjit             !! Time step number 
    LOGICAL, INTENT(in)                             :: ldrestart_read   !! Logical for restart file to read
    INTEGER(i_std), INTENT(in)                       :: kjpindex         !! Domain size
    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index            !! Indeces of the points on the map
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in):: veget            !! Vegetation distribution
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (in)  :: lalo             !! Geographical coordinates
    INTEGER(i_std),DIMENSION (kjpindex,8), INTENT(in):: neighbours       !! neighoring grid points if land
    REAL(r_std), DIMENSION (kjpindex,2), INTENT(in)  :: resolution       !! size in x an y of the grid (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: contfrac         ! Fraction of land in each grid box.
    INTEGER(i_std), INTENT(in)                       :: rest_id          !! Restart file identifier
    ! input fields
    ! output scalar
    ! output fields

    ! local declaration
    INTEGER(i_std)                                  :: ji
    INTEGER(i_std)                                  :: ier
    CHARACTER(LEN=80)                               :: var_name         !! To store variables names for I/O
    ! initialisation
    IF (l_first_condveg) THEN 
       !       
       l_first_condveg=.FALSE.
       !
       !  Allocate variables which have to
       ! 
       ALLOCATE (soilalb_dry(kjpindex,2),stat=ier)
       IF (ier.NE.0) THEN
          WRITE (numout,*) ' error in soilalb_dry allocation. We stop.We need kjpindex*2 words = ',kjpindex*2
          STOP 'condveg_init'
       END IF
       soilalb_dry(:,:) = val_exp

       ALLOCATE (soilalb_wet(kjpindex,2),stat=ier)
       IF (ier.NE.0) THEN
          WRITE (numout,*) ' error in soilalb_wet allocation. We stop.We need kjpindex*2 words = ',kjpindex*2
          STOP 'condveg_init'
       END IF
       soilalb_wet(:,:) = val_exp

        ! Ajout Nathalie
       ALLOCATE (soilalb_moy(kjpindex,2),stat=ier)
        IF (ier.NE.0) THEN
          WRITE (numout,*) ' error in soilalb_moy allocation. We stop.We need kjpindex*2 words = ',kjpindex*2
          STOP 'condveg_init'
        END IF
       soilalb_moy(:,:) = val_exp

       ALLOCATE (albedo_snow(kjpindex),stat=ier)
        IF (ier.NE.0) THEN
          WRITE (numout,*) ' error in albedo_snow allocation. We stop.We need kjpindex words = ',kjpindex
          STOP 'condveg_init'
        END IF

       ALLOCATE (albedo_glob(kjpindex),stat=ier)
        IF (ier.NE.0) THEN
          WRITE (numout,*) ' error in albedo_glob allocation. We stop.We need kjpindex words = ',kjpindex
          STOP 'condveg_init'
        END IF

       ALLOCATE (alb_bare(kjpindex,2),stat=ier)
       IF (ier.NE.0) THEN
          WRITE (numout,*) ' error in alb_bare allocation. We stop.We need kjpindex*2 words = ',kjpindex*2
          STOP 'condveg_init'
       END IF
       alb_bare(:,:) = val_exp

       ALLOCATE (alb_veget(kjpindex,2),stat=ier)
       IF (ier.NE.0) THEN
          WRITE (numout,*) ' error in alb_veget allocation. We stop.We need kjpindex*2 words = ',kjpindex*2
          STOP 'condveg_init'
       END IF
       alb_veget(:,:) = val_exp
       !
       !  Get the bare soil albedo
       !
       !
       var_name= 'soilalbedo_dry'
       CALL ioconf_setatt('UNITS', '-')
       CALL ioconf_setatt('LONG_NAME','Dry bare soil albedo')
       CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_dry, "gather", nbp_glo, index_g)

       var_name= 'soilalbedo_wet'
       CALL ioconf_setatt('UNITS', '-')
       CALL ioconf_setatt('LONG_NAME','Wet bare soil albedo')
       CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_wet, "gather", nbp_glo, index_g)

       var_name= 'soilalbedo_moy'
       CALL ioconf_setatt('UNITS', '-')
       CALL ioconf_setatt('LONG_NAME','Mean bare soil albedo')
       CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_moy, "gather", nbp_glo, index_g)
       !
       IF ( MINVAL(soilalb_wet) .EQ. MAXVAL(soilalb_wet) .AND. MAXVAL(soilalb_wet) .EQ. val_exp .OR.&
            & MINVAL(soilalb_dry) .EQ. MAXVAL(soilalb_dry) .AND. MAXVAL(soilalb_dry) .EQ. val_exp .OR.&
            & MINVAL(soilalb_moy) .EQ. MAXVAL(soilalb_moy) .AND. MAXVAL(soilalb_moy) .EQ. val_exp) THEN
          CALL condveg_soilalb(kjpindex, lalo, neighbours, resolution, contfrac, soilalb_dry,soilalb_wet)
          WRITE(numout,*) '---> val_exp ', val_exp
          WRITE(numout,*) '---> ALBEDO_wet VIS:', MINVAL(soilalb_wet(:,ivis)), MAXVAL(soilalb_wet(:,ivis))
          WRITE(numout,*) '---> ALBEDO_wet NIR:', MINVAL(soilalb_wet(:,inir)), MAXVAL(soilalb_wet(:,inir))
          WRITE(numout,*) '---> ALBEDO_dry VIS:', MINVAL(soilalb_dry(:,ivis)), MAXVAL(soilalb_dry(:,ivis))
          WRITE(numout,*) '---> ALBEDO_dry NIR:', MINVAL(soilalb_dry(:,inir)), MAXVAL(soilalb_dry(:,inir))
          WRITE(numout,*) '---> ALBEDO_moy VIS:', MINVAL(soilalb_moy(:,ivis)), MAXVAL(soilalb_moy(:,ivis))
          WRITE(numout,*) '---> ALBEDO_moy NIR:', MINVAL(soilalb_moy(:,inir)), MAXVAL(soilalb_moy(:,inir))
       ENDIF
       !
    ELSE 
       WRITE (numout,*) ' l_first_condveg false . we stop '
       STOP 'condveg_init'
    ENDIF
    !! test de commentaires

    IF (long_print) WRITE (numout,*) ' condveg_init done '

  END SUBROUTINE condveg_init
  !!
  !!
  SUBROUTINE condveg_clear  ()

      l_first_condveg=.TRUE.
       
       IF (ALLOCATED (soilalb_dry)) DEALLOCATE (soilalb_dry)
       IF (ALLOCATED(soilalb_wet))  DEALLOCATE (soilalb_wet)
       ! Ajout Nathalie
       IF (ALLOCATED(soilalb_moy))  DEALLOCATE (soilalb_moy)
       IF (ALLOCATED(albedo_snow))  DEALLOCATE (albedo_snow)
       IF (ALLOCATED(albedo_glob))  DEALLOCATE (albedo_glob)
       IF (ALLOCATED(alb_bare))  DEALLOCATE (alb_bare)
       IF (ALLOCATED(alb_veget))  DEALLOCATE (alb_veget)
       !
  END SUBROUTINE condveg_clear

  !! Algorithm:
  !! - initialisation of local arry
  !! - reads map for emissivity, albedo or roughness
  !!
  SUBROUTINE condveg_var_init  (ldrestart_read, kjpindex, veget, veget_max, frac_nobio, totfrac_nobio,&
       & drysoil_frac, zlev, height, deadleaf_cover, emis, albedo, z0, roughheight)

    ! interface description
    ! input scalar 
    INTEGER(i_std), INTENT(in)                       :: kjpindex         !! Domain size
    LOGICAL, INTENT(in)                             :: ldrestart_read   !! Logical for _restart_ file to read
    ! input fields
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget            !! Fraction of vegetation types
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max        !! Fraction of vegetation type
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio    !! Fraction of continental ice, lakes, ...
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: totfrac_nobio    !! Total fraction of continental ice+lakes+...
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: drysoil_frac     !! Fraction of visibly Dry soil(between 0 and 1)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)    :: zlev             !! Height of first layer
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height           !! Vegetation Height (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: deadleaf_cover   !! Fraction of soil covered by dead leaves
    ! output scalar
    ! output fields
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)   :: emis             !! Emissivity
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo           !! Albedo, vis(1) and nir(2)
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)   :: z0               !! Roughness
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)   :: roughheight      !! Effective height for roughness
    !
    ! local declaration
    INTEGER(i_std)                                 :: ier               !! Error code
    INTEGER(i_std)                                 :: jv
    ! initialisation of variables

    !!
    !! calculs de emis
    !!

    IF ( impaze ) THEN
       !
       emis(:) = emis_scal
       !
    ELSE
       !  Some day it will be moisture dependent
       emis_scal = un
       
       emis(:) = emis_scal
    ENDIF

    !!
    !! calculs de albedo
    !!
    !
    IF ( impaze ) THEN
       !
       albedo(:,ivis) = albedo_scal(ivis)
       albedo(:,inir) = albedo_scal(inir)
       !
    ELSE
       !
       CALL condveg_albcalc (kjpindex,deadleaf_cover,veget,drysoil_frac,albedo)
       !
    ENDIF
    !!
    !! calculs de z0
    !!
    !!
    IF ( impaze ) THEN
      !
      z0(:) = z0_scal
      roughheight(:) = roughheight_scal
      !
    ELSE
      !
       IF ( z0cdrag_ave ) THEN
          CALL condveg_z0cdrag(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
               &               height, z0, roughheight)
       ELSE
          CALL condveg_z0logz(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
               &              z0, roughheight)
       ENDIF
      !
    ENDIF
    !!
    !!
    IF (long_print) WRITE (numout,*) ' condveg_var_init done '

  END SUBROUTINE condveg_var_init
  !!
  !! Algorithm:
  !! - Simply update the emissivity, albedo and roughness fields
  !!
  SUBROUTINE condveg_var_update  (ldrestart_read, kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, &
       & drysoil_frac, zlev, height, deadleaf_cover, emis, albedo, z0, roughheight)

    ! interface description
    ! input scalar 
    INTEGER(i_std), INTENT(in)                       :: kjpindex         !! Domain size
    LOGICAL, INTENT(in)                             :: ldrestart_read   !! Logical for _restart_ file to read
    ! input fields
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget            !! Fraction of vegetation types
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max        !! Fraction of vegetation type
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio    !! Fraction of continental ice, lakes, ...
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: totfrac_nobio    !! Total fraction of continental ice+lakes+...
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: drysoil_frac     !! Fraction of visibly Dry soil(between 0 and 1)
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)    :: zlev             !! Height of first layer
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height           !! Vegetation Height (m)
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)     :: deadleaf_cover   !! Fraction of soil covered by dead leaves
    ! output scalar
    ! output fields
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)   :: emis             !! Emissivity
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo           !! Albedo, vis(1) and nir(2)
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)   :: z0               !! Roughness
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)   :: roughheight      !! Effective height for roughness
    !
    ! local
    INTEGER(i_std)                                  :: ji, jv
    !!
    !! calculs de emis
    !!

    emis(:) = emis_scal
    
    !!
    !! calculs de albedo
    !!
    !
    IF ( impaze ) THEN
       !
       albedo(:,ivis) = albedo_scal(ivis)
       albedo(:,inir) = albedo_scal(inir)
       !
    ELSE
       !
       CALL condveg_albcalc (kjpindex,deadleaf_cover,veget,drysoil_frac,albedo)
       !
    ENDIF
    !
    !!
    !! Calculs de la rugosite
    !!
    
    IF ( impaze ) THEN

       DO ji = 1, kjpindex
         z0(ji) = z0_scal
         roughheight(ji) = roughheight_scal
      ENDDO

    ELSE
       !
       IF ( z0cdrag_ave ) THEN
          CALL condveg_z0cdrag (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, height, &
               &                z0, roughheight)
       ELSE
          CALL condveg_z0logz (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
               &               z0, roughheight)
       ENDIF
       !
    ENDIF

    IF (long_print) WRITE (numout,*) ' condveg_var_update done '

  END SUBROUTINE condveg_var_update

  !! Algorithm:
  !! - The snow albedo is updated by the snow within the mesh
  !! This is done as a function of snow mass, snow age and vegetation type
  !! The model is described in Chalita 1992
  !!
  SUBROUTINE condveg_snow  (ldrestart_read, kjpindex, veget, frac_nobio, totfrac_nobio, &
                            snow, snow_age, snow_nobio, snow_nobio_age, albedo, albedo_snow, albedo_glob)

    ! interface description
    ! input scalar 
    INTEGER(i_std), INTENT(in)                          :: kjpindex         !! Domain size
    LOGICAL, INTENT(in)                                :: ldrestart_read   !! Logical for _restart_ file to read
    ! input fields
    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)   :: veget            !! Fraction of vegetation types
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio       !! Fraction of continental ice, lakes, ...
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)        :: totfrac_nobio    !! Total fraction of continental ice+lakes+ ...
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)        :: snow             !! Snow mass [Kg/m^2] in vegetation
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio       !! Snow mass [Kg/m^2]on ice, lakes, ...
    REAL(r_std),DIMENSION (kjpindex), INTENT(in)        :: snow_age         !! Snow age
    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age   !! Snow age on ice, lakes, ...
    ! output scalar
    ! output fields
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (inout)  :: albedo           !! Albedo
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)      :: albedo_snow      !! Albedo de la neige
    REAL(r_std),DIMENSION (kjpindex), INTENT (out)      :: albedo_glob      !! Mean albedo
    !
    ! local declaration
    INTEGER(i_std)                                  :: ji, jv, jb       !! Loop index
    REAL(r_std), DIMENSION(kjpindex)                 :: frac_snow_veg    !! Fraction of snow on vegetation
    REAL(r_std), DIMENSION(kjpindex,nnobio)          :: frac_snow_nobio  !! Fraction of snow on ice, lakes, ...
    REAL(r_std), DIMENSION(kjpindex)                 :: snowa_veg        !! Total albedo of snow covered area on vegetation
    REAL(r_std), DIMENSION(kjpindex,nnobio)          :: snowa_nobio      !! albedo of snow covered area on ice, lakes, ...
    REAL(r_std), DIMENSION(kjpindex)                 :: fraction_veg     !! total vegetation fraction
    REAL(r_std), DIMENSION(kjpindex)                 :: agefunc_veg      !! age dependency of snow albedo on vegetation
    REAL(r_std), DIMENSION(kjpindex,nnobio)          :: agefunc_nobio    !! age dependency of snow albedo on ice, lakes, ..;
    REAL(r_std)                                      :: alb_nobio
    !

    DO ji = 1, kjpindex
     albedo_snow(ji) = zero
     albedo_glob(ji) = zero
    ENDDO

    IF (ABS(fixed_snow_albedo - undef_sechiba) .GT. EPSILON(undef_sechiba)) THEN
      snowa_veg(:) = fixed_snow_albedo
      snowa_nobio(:,:) = fixed_snow_albedo
    ELSE
      !
      ! calculate first age dependence
      !
      DO ji = 1, kjpindex
        agefunc_veg(ji) = EXP(-snow_age(ji)/tcst_snowa)
      ENDDO
      !
      !
      DO jv = 1, nnobio
        DO ji = 1, kjpindex
          agefunc_nobio(ji,jv) = EXP(-snow_nobio_age(ji,jv)/tcst_snowa)
        ENDDO
      ENDDO
      !
      ! snow albedo on vegetated surfaces
      !
      fraction_veg(:) = un - totfrac_nobio(:)
      snowa_veg(:) = zero
      DO jv = 1, nvm
        DO ji = 1, kjpindex
          IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
            snowa_veg(ji) = snowa_veg(ji) + &
              veget(ji,jv)/fraction_veg(ji) * ( snowa_ini(jv)+snowa_dec(jv)*agefunc_veg(ji) )
          ENDIF
        ENDDO
      ENDDO
      !
      ! snow albedo on other surfaces
      !
      DO jv = 1, nnobio
        DO ji = 1, kjpindex
          snowa_nobio(ji,jv) = ( snowa_ini(1) + snowa_dec(1) * agefunc_nobio(ji,jv) ) 
        ENDDO
      ENDDO
    ENDIF

    frac_snow_veg(:) = MIN(MAX(snow(:),zero)/(MAX(snow(:),zero)+snowcri_alb),un)
    DO jv = 1, nnobio
      frac_snow_nobio(:,jv) = MIN(MAX(snow_nobio(:,jv),zero)/(MAX(snow_nobio(:,jv),zero)+snowcri_alb),un)
    ENDDO


    DO jb = 1, 2
      !
      albedo(:,jb) = ( fraction_veg(:) ) * &
                         ( (un-frac_snow_veg(:)) * albedo(:,jb) + &
                           ( frac_snow_veg(:)  ) * snowa_veg(:)    )
      albedo_snow(:) =  albedo_snow(:) + (fraction_veg(:)) * (frac_snow_veg(:)) * snowa_veg(:)
      !
      DO jv = 1, nnobio
        ! 
        IF ( jv .EQ. iice ) THEN
          alb_nobio = alb_ice(jb)
        ELSE
          WRITE(numout,*) 'jv=',jv
          STOP 'DO NOT KNOW ALBEDO OF THIS SURFACE TYPE'
        ENDIF
        !
        albedo(:,jb) = albedo(:,jb) + &
                       ( frac_nobio(:,jv) ) * &
                         ( (un-frac_snow_nobio(:,jv)) * alb_nobio + &
                           ( frac_snow_nobio(:,jv)  ) * snowa_nobio(:,jv)   )
        albedo_snow(:) = albedo_snow(:) + &
                         ( frac_nobio(:,jv) ) * ( frac_snow_nobio(:,jv) ) * &
                           snowa_nobio(:,jv)
        albedo_glob(:) = albedo_glob(:) + albedo(:,jb)
        !
      ENDDO
      !
    END DO
    !
    DO ji = 1, kjpindex
      albedo_snow(ji) = (albedo_snow(ji))/2.
      albedo_glob(ji) = (albedo_glob(ji))/2.
    ENDDO
    
    IF (long_print) WRITE (numout,*) ' condveg_snow done '

  END SUBROUTINE condveg_snow
  !
  !
  !
  SUBROUTINE condveg_soilalb(nbpt, lalo, neighbours, resolution, contfrac, soilalb_dry,soilalb_wet)
    !
    !
    !   This subroutine should read the soil color maps from the Henderson-Sellers & Wilson database. This data
    !   is then interpolated to the models resolution and transformed into dry and wet albedos. We have chosen to
    !  do both these operations at the same time as one can average the albedo but not the color types.
    !
    !  We make the assumption in this code that the grid of the data is regular and that it covers the globe. 
    !  For the model grid we base the calculation of the borders of the grid on the resolution.
    !
    !
    !  0.1 INPUT
    !
  INTEGER(i_std), INTENT(in)    :: nbpt                ! Number of points for which the data needs to be interpolated
  REAL(r_std), INTENT(in)       :: lalo(nbpt,2)        ! Vector of latitude and longitudes (beware of the order !)
  INTEGER(i_std), INTENT(in)    :: neighbours(nbpt,8)  ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
  REAL(r_std), INTENT(in)       :: resolution(nbpt,2)  ! The size in km of each grid-box in X and Y
  REAL(r_std), INTENT(in)       :: contfrac(nbpt)     ! Fraction of land in each grid box.
  REAL(r_std), INTENT(inout)    :: soilalb_dry(nbpt,2) ! albedo for the dry bare soil
  REAL(r_std), INTENT(inout)    :: soilalb_wet(nbpt,2) ! albedo for the wet bare soil

  !
  !
  !  0.3 LOCAL
  !
  INTEGER(i_std)    :: nbvmax
  !
  CHARACTER(LEN=80) :: filename
  INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, fopt, ilf, lastjp, nbexp
  REAL(r_std) :: lev(1), date, dt, coslat, sgn
  INTEGER(i_std) :: itau(1)
  REAL(r_std), ALLOCATABLE, DIMENSION(:,:)        :: lat_rel, lon_rel, soilcol
  REAL(r_std), ALLOCATABLE, DIMENSION(:,:)        :: sub_area
  INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:)  :: sub_index
  INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:)    :: mask
  CHARACTER(LEN=30) :: callsign
  !   
  LOGICAL                  :: ok_interpol = .FALSE. ! optionnal return of aggregate_2d
  !   
  INTEGER                  :: ALLOC_ERR
  !
  !
  !  Needs to be a configurable variable
  !
  !
  !Config Key   = SOILALB_FILE
  !Config Desc  = Name of file from which the bare soil albedo
  !Config Def   = soils_param.nc
  !Config If    = NOT(IMPOSE_AZE)
  !Config Help  = The name of the file to be opened to read the soil types from 
  !Config         which we derive then the bare soil albedos. This file is 1x1 
  !Config         deg and based on the soil colors defined by Wilson and Henderson-Seller.
  !Config Units = [FILE]
  !
  filename = 'soils_param.nc'
  CALL getin_p('SOILALB_FILE',filename)
  !
  IF (is_root_prc) CALL flininfo(filename,iml, jml, lml, tml, fid)
  CALL bcast(iml)
  CALL bcast(jml)
  CALL bcast(lml)
  CALL bcast(tml)
  !
  ! soils_param.nc file is 1° soit texture file.
  !
  ALLOC_ERR=-1
  ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
     WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
     STOP 
  ENDIF
  ALLOC_ERR=-1
  ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
     WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
     STOP 
  ENDIF
  ALLOC_ERR=-1
  ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
     WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
     STOP 
  ENDIF
  ALLOC_ERR=-1
  ALLOCATE(soilcol(iml,jml), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
     WRITE(numout,*) "ERROR IN ALLOCATION of soiltext : ",ALLOC_ERR
     STOP 
  ENDIF
  !
  IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid)
  CALL bcast(lon_rel)
  CALL bcast(lat_rel)
  CALL bcast(lev)
  CALL bcast(itau)
  CALL bcast(date)
  CALL bcast(dt)
  !
  IF (is_root_prc) CALL flinget(fid, 'soilcolor', iml, jml, lml, tml, 1, 1, soilcol)
  CALL bcast(soilcol)
  !
  IF (is_root_prc) CALL flinclo(fid)
  !
  ! Mask of permitted variables.
  !
  mask(:,:) = zero
  DO ip=1,iml
     DO jp=1,jml
        IF (soilcol(ip,jp) > min_sechiba) THEN
           mask(ip,jp) = un
        ENDIF
     ENDDO
  ENDDO
  !
  nbvmax = 200
  !
  callsign = 'Soil color map'
  !
  DO WHILE ( .NOT. ok_interpol )
     WRITE(numout,*) "Projection arrays for ",callsign," : "
     WRITE(numout,*) "nbvmax = ",nbvmax
     !
     ALLOC_ERR=-1
     ALLOCATE(sub_area(nbpt,nbvmax), STAT=ALLOC_ERR)
     IF (ALLOC_ERR/=0) THEN
         WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
        STOP 
     ENDIF
     sub_area(:,:)=zero
     ALLOC_ERR=-1
     ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
     IF (ALLOC_ERR/=0) THEN
        WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
        STOP 
     ENDIF
     sub_index(:,:,:)=0
     !
     CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
          &                iml, jml, lon_rel, lat_rel, mask, callsign, &
          &                nbvmax, sub_index, sub_area, ok_interpol)
     !
     IF ( .NOT. ok_interpol ) THEN
        DEALLOCATE(sub_area)
        DEALLOCATE(sub_index)
     ENDIF
     !
     nbvmax = nbvmax * 2
  ENDDO
  !
  nbexp = 0
  !
  !    Check that we found some points
  !
  soilalb_dry(:,:) = zero
  soilalb_wet(:,:) = zero
  ! Ajout Nathalie
  soilalb_moy(:,:) = zero
  !
  DO ib=1,nbpt
     !
     ! GO through the point we have found
     !
     !
     fopt =  COUNT(sub_area(ib,:) > zero)
     !
     IF ( fopt .EQ. 0) THEN
        nbexp = nbexp + 1
        soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
        soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
        soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
        soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
        ! Ajout Nathalie
        soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
        soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
     ELSE
        sgn = zero
        !
        !   Compute the average bare soil albedo parameters
        !
        DO ilf = 1,fopt
           !
           ip = sub_index(ib,ilf,1)
           jp = sub_index(ib,ilf,2)
           !
           IF ( NINT(soilcol(ip,jp)) .LE. classnb) THEN
              soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis) + vis_dry(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
              soilalb_dry(ib,inir) = soilalb_dry(ib,inir) + nir_dry(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
              soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis) + vis_wet(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
              soilalb_wet(ib,inir) = soilalb_wet(ib,inir) + nir_wet(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
              ! Ajout Nathalie
              soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis) + albsoil_vis(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
              soilalb_moy(ib,inir) = soilalb_moy(ib,inir) + albsoil_nir(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
              sgn = sgn + sub_area(ib,ilf)
           ELSE
              WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
              STOP
           ENDIF
           !
        ENDDO
        !
        ! Normalize the surface
        !
        IF ( sgn .LT. min_sechiba) THEN
           nbexp = nbexp + 1
           soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
           soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
           soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
           soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
           ! Ajout Nathalie
           soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
           soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
        ELSE
           soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis)/sgn
           soilalb_dry(ib,inir) = soilalb_dry(ib,inir)/sgn
           soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis)/sgn
           soilalb_wet(ib,inir) = soilalb_wet(ib,inir)/sgn
           ! Ajout Nathalie
           soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis)/sgn
           soilalb_moy(ib,inir) = soilalb_moy(ib,inir)/sgn           
        ENDIF
        !
     ENDIF
     !
  ENDDO
  !
  IF ( nbexp .GT. 0 ) THEN
     WRITE(numout,*) 'CONDVEG_soilalb : The interpolation of the bare soil albedo had ', nbexp
     WRITE(numout,*) 'CONDVEG_soilalb : points without data. This are either coastal points or'
     WRITE(numout,*) 'CONDVEG_soilalb : ice covered land.'
     WRITE(numout,*) 'CONDVEG_soilalb : The problem was solved by using the average of all soils'
     WRITE(numout,*) 'CONDVEG_soilalb : in dry and wet conditions'
  ENDIF
  !
  DEALLOCATE (lat_rel)
  DEALLOCATE (lon_rel)
  DEALLOCATE (mask)
  DEALLOCATE (sub_index)
  DEALLOCATE (sub_area)
  DEALLOCATE (soilcol)
  !
  !
  RETURN
  !
  END SUBROUTINE condveg_soilalb
  !
  !
  !
  SUBROUTINE condveg_z0logz (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
       &                     z0, roughheight)

    !
    ! 0. Declarations
    !

    ! 0.1 Input
    INTEGER(i_std), INTENT(in)                       :: kjpindex
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max
    REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio
    REAL(r_std), DIMENSION(kjpindex), INTENT(in)     :: totfrac_nobio
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height
    
    ! 0.2 Output
    REAL(r_std), DIMENSION(kjpindex), INTENT(out)    :: z0
    REAL(r_std), DIMENSION(kjpindex), INTENT(out)    :: roughheight
    
    ! 0.3 Local
    INTEGER(i_std)                                 :: jv
    REAL(r_std), DIMENSION(kjpindex)                 :: sumveg, ave_height
    REAL(r_std), DIMENSION(kjpindex)                 :: d_veg, zhdispl
    REAL(r_std)                                      :: z0_nobio
    
    z0(:) = veget(:,1) * LOG(z0_bare)
    sumveg(:) = veget(:,1)
    ave_height(:) = zero
    !
    DO jv = 2, nvm
       !
       IF ( is_tree(jv) ) THEN
          ! tree trunks influence the atmosphere even when there are no leaves
          d_veg(:) = veget_max(:,jv)
       ELSE
          ! grasses only have an influence if they are really there!
          d_veg(:) = veget(:,jv)
       ENDIF
       !
       z0(:) = z0(:) + d_veg(:) * &
            LOG( MAX(height(:,jv)*z0_over_height,z0_bare) )
       sumveg(:) = sumveg(:) + d_veg(:)
       !
       ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
       !
    ENDDO
    !
    WHERE ( sumveg(:) > zero ) z0(:) = z0(:) / sumveg(:)
    !
    z0(:) = (un - totfrac_nobio(:)) * z0(:)
    !
    DO jv = 1, nnobio
       !
       IF ( jv .EQ. iice ) THEN
          z0_nobio = z0_ice
       ELSE
          WRITE(numout,*) 'jv=',jv
          STOP 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
       ENDIF
       !
       z0(:) = z0(:) + frac_nobio(:,jv) * LOG(z0_nobio)
       !
    ENDDO
    !
    z0(:) = EXP( z0(:) )
    !
    ! Temporarily we compute the zero plane displacement height
    !
    zhdispl(:) = ave_height(:) * height_displacement
    !
    ! In order to get a variable independent of the height of the
    ! vegetation we compute what we call the effective roughness height.
    ! This is the height over which the roughness acts. It combines the
    ! zero plane displacement height and the vegetation height.
    !
    roughheight(:) = ave_height(:) - zhdispl(:)
    !
  END SUBROUTINE condveg_z0logz
  !
  !
  !
  SUBROUTINE condveg_z0cdrag (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, &
       &                      z0, roughheight)
    
    !
    ! 0. Declarations
    !
    
    ! 0.1 Input
    INTEGER(i_std), INTENT(in)                       :: kjpindex
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max
    REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio
    REAL(r_std), DIMENSION(kjpindex), INTENT(in)     :: totfrac_nobio
    REAL(r_std),DIMENSION (kjpindex), INTENT (in)    :: zlev             !! Height of first layer
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height    
    
    ! 0.2 Output
    REAL(r_std), DIMENSION(kjpindex), INTENT(out)    :: z0
    REAL(r_std), DIMENSION(kjpindex), INTENT(out)    :: roughheight
    
    ! 0.3 Local
    INTEGER(i_std)                                 :: jv
    REAL(r_std), DIMENSION(kjpindex)                 :: sumveg, ztmp, ave_height
    REAL(r_std), DIMENSION(kjpindex)                 :: d_veg, zhdispl
    REAL(r_std)                                      :: z0_nobio
    !
    ! The grid average z0 is computed by averaging the neutral drag coefficients.
    ! This is pretty straight forward except for the reference level which needs
    ! to be chosen.
    !
    ! We need a reference lever high enough above the canopy else we get into
    ! singularities of the LOG.
    !
    ztmp(:) = MAX(10., zlev(:))
    !
    z0(:) = veget(:,1) * (ct_karman/LOG(ztmp(:)/z0_bare))**2
    sumveg(:) = veget(:,1)
    ave_height(:) = zero
    !
    DO jv = 2, nvm
       !
       IF ( is_tree(jv) ) THEN
          ! tree trunks influence the atmosphere even when there are no leaves
          d_veg(:) = veget_max(:,jv)
       ELSE
          ! grasses only have an influence if they are really there!
          d_veg(:) = veget(:,jv)
       ENDIF
       !
       z0(:) = z0(:) + d_veg(:) * (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)*z0_over_height,z0_bare)))**2
       sumveg(:) = sumveg(:) + d_veg(:)
       !
       ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
       !
    ENDDO
    !
    WHERE ( sumveg(:) .GT. zero ) z0(:) = z0(:) / sumveg(:)
    !
    z0(:) = (un - totfrac_nobio(:)) * z0(:)
    !
    DO jv = 1, nnobio
       !
       IF ( jv .EQ. iice ) THEN
          z0_nobio = z0_ice
       ELSE
          WRITE(numout,*) 'jv=',jv
          STOP 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
       ENDIF
       !
       z0(:) = z0(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
       !
    ENDDO
    !
    z0(:) = ztmp(:) / EXP(ct_karman/SQRT(z0(:)))
    !
    ! Temporarily we compute the zero plane displacement height
    !
    zhdispl(:) = ave_height(:) * height_displacement
    !
    ! In order to get a variable independent of the height of the
    ! vegetation we compute what we call the effective roughness height.
    ! This is the height over which the roughness acts. It combines the
    ! zero plane displacement height and the vegetation height.
    !
    roughheight(:) = ave_height(:) - zhdispl(:)
    !
  END SUBROUTINE condveg_z0cdrag
  !
  !
  !
  SUBROUTINE condveg_albcalc (kjpindex,deadleaf_cover,veget,drysoil_frac,albedo)

    !
    ! 0. Declarations
    !

    ! 0.1 Input
    INTEGER(i_std), INTENT(in)                       :: kjpindex        !! Domain size
    REAL(r_std), DIMENSION(kjpindex), INTENT(in)     :: deadleaf_cover  !! Fraction of soil covered by dead leaves
    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget           !! Fraction of vegetation types
    REAL(r_std), DIMENSION(kjpindex), INTENT(in)     :: drysoil_frac    !! Fraction of visibly Dry soil(between 0 and 1)
    
    ! 0.2 Output
    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo  !! Albedo, vis(1) and nir(2)
    ! 0.3 Local
!    REAL(r_std),DIMENSION (kjpindex)                 :: alb_bare
    REAL(r_std),DIMENSION (nvm,2)                    :: alb_leaf_tmp
    INTEGER(i_std)                                   :: ks, jv
    !
    alb_leaf_tmp(:,ivis) = alb_leaf_vis(:)
    alb_leaf_tmp(:,inir) = alb_leaf_nir(:)
    !
    !
    DO ks = 1, 2
       !
       IF ( alb_bare_model ) THEN
          alb_bare(:,ks) = soilalb_wet(:,ks) + drysoil_frac(:) * (soilalb_dry(:,ks) -  soilalb_wet(:,ks))
       ELSE
          ! Nouvelle formulation Nathalie, sans dependance en drysoil_frac.
          alb_bare(:,ks) = soilalb_moy(:,ks)
       ENDIF
       !
       ! Correction Nathalie le 12 Avril 2006 - suppression de la dependance en deadleaf_cover
       !albedo(:,ks) = veget(:,1) * ( (un-deadleaf_cover(:))*alb_bare(:) + &
       !                              deadleaf_cover(:)*alb_deadleaf(ks)    )
       albedo(:,ks) = veget(:,1) * alb_bare(:,ks)

       ! vegetation
       alb_veget(:,ks) = zero
       DO jv = 2, nvm
          albedo(:,ks) = albedo(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
          alb_veget(:,ks) = alb_veget(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
       ENDDO
       !
    ENDDO
    
  END SUBROUTINE condveg_albcalc

END MODULE condveg