coupling_tebn.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     ###############################################################################
SUBROUTINE coupling_teb_n (DTCO, DST, SLT, TOP, SB, G, CHT, NT, TPN, TIR, BOP, NB, TD, GDM, GRM, &
                           HPROGRAM, HCOUPLING, PTSTEP, KYEAR, KMONTH, KDAY, PTIME, KI, KSV,&
                           KSW, PTSUN, PZENITH, PAZIM, PZREF, PUREF, PZS, PU, PV, PQA, PTA, &
                           PRHOA, PSV, PCO2, HSV, PRAIN, PSN, PLW, PDIR_SW, PSCA_SW,        &
                           PSW_BANDS, PPS, PPA, PSFTQ, PSFTH, PSFTS, PSFCO2, PSFU, PSFV,    &
                           PTRAD, PDIR_ALB, PSCA_ALB, PEMIS, PTSURF, PZ0, PZ0H, PQSURF,     &
                           PPEW_A_COEF, PPEW_B_COEF, PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, &
                           PPEQ_B_COEF, HTEST     )
!     ###############################################################################
!
!!****  *COUPLING_TEB_n * - Driver for TEB 
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!                  10/2005 (G.Pigeon) transfer of domestic heating
!!      S. Riette   06/2009 Initialisation of XT, XQ, XU and XTKE on canopy levels
!!      S. Riette   01/2010 Use of interpol_sbl to compute 10m wind diagnostic
!!      G. Pigeon   09/2012 CCH_BEM, ROUGH_WALL, ROUGH_ROOF for building conv. coef
!!      G. Pigeon   10/2012 XF_WIN_WIN as arg. of TEB_GARDEN
!!      B. Decharme 09/2012 New wind implicitation
!!      J. Escobar  09/2012 KI not allowed without-interface , replace by KI
!!      V. Masson   08/2013 adds solar panels & occupation calendar
!!      B. Decharme  04/2013 new coupling variables
!!---------------------------------------------------------------
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_dst_n, ONLY : dst_t
USE modd_slt_n, ONLY : slt_t
!
USE modd_ch_teb_n, ONLY : ch_teb_t
USE modd_canopy_n, ONLY: canopy_t
USE modd_sfx_grid_n, ONLY : grid_t
USE modd_teb_option_n, ONLY : teb_options_t
USE modd_teb_panel_n, ONLY : teb_panel_t
USE modd_teb_irrig_n, ONLY : teb_irrig_t
USE modd_teb_n, ONLY : teb_np_t
USE modd_surfex_n, ONLY : teb_diag_t
USE modd_bem_option_n, ONLY : bem_options_t
USE modd_bem_n, ONLY : bem_np_t
!
USE modd_surfex_n, ONLY : teb_garden_model_t
USE modd_surfex_n, ONLY : teb_greenroof_model_t
!
USE modd_reprod_oper, ONLY : cimplicit_wind
!
USE modd_csts,         ONLY : xrd, xcpd, xp00, xlvtt, xpi, xkarman, xg
USE modd_surf_par,     ONLY : xundef
!                            
USE modd_dst_surf
USE modd_slt_surf
!
USE mode_dslt_surf
USE mode_thermos
USE mode_sbls
!
USE modi_average_rad
USE modi_sm10
USE modi_add_forecast_to_date_surf
USE modi_diag_inline_teb_n
USE modi_cumul_diag_teb_n
USE modi_ch_aer_dep
USE modi_ch_dep_town
USE modi_dslt_dep
USE modi_teb_garden
USE modi_teb_canopy
! 
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
USE modi_canopy_evol
USE modi_canopy_grid_update
USE modi_utci_teb
USE modi_utcic_stress
USE modi_circumsolar_rad
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
TYPE(dst_t), INTENT(INOUT) :: DST
TYPE(slt_t), INTENT(INOUT) :: SLT
!
TYPE(ch_teb_t), INTENT(INOUT) :: CHT 
TYPE(canopy_t), INTENT(INOUT) :: SB
TYPE(grid_t), INTENT(INOUT) :: G
TYPE(teb_options_t), INTENT(INOUT) :: TOP
TYPE(teb_panel_t), INTENT(INOUT) :: TPN
TYPE(teb_irrig_t), INTENT(INOUT) :: TIR
TYPE(teb_np_t), INTENT(INOUT) :: NT
!
TYPE(teb_diag_t), INTENT(INOUT) :: TD
!
TYPE(bem_options_t), INTENT(INOUT) :: BOP 
TYPE(bem_np_t), INTENT(INOUT) :: NB
!
TYPE(teb_garden_model_t), INTENT(INOUT) :: GDM
TYPE(teb_greenroof_model_t), INTENT(INOUT) :: GRM
!
 CHARACTER(LEN=6),    INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
 CHARACTER(LEN=1),    INTENT(IN)  :: HCOUPLING ! type of coupling
                                              ! 'E' : explicit
                                              ! 'I' : implicit
INTEGER,             INTENT(IN)  :: KYEAR     ! current year (UTC)
INTEGER,             INTENT(IN)  :: KMONTH    ! current month (UTC)
INTEGER,             INTENT(IN)  :: KDAY      ! current day (UTC)
REAL,                INTENT(IN)  :: PTIME     ! current time since midnight (UTC, s)
INTEGER,             INTENT(IN)  :: KI        ! number of points
INTEGER,             INTENT(IN)  :: KSV       ! number of scalars
INTEGER,             INTENT(IN)  :: KSW       ! number of short-wave spectral bands
REAL, DIMENSION(KI), INTENT(IN)  :: PTSUN     ! solar time                    (s from midnight)
REAL,                INTENT(IN)  :: PTSTEP    ! atmospheric time-step                 (s)
REAL, DIMENSION(KI), INTENT(IN)  :: PZREF     ! height of T,q forcing                 (m)
REAL, DIMENSION(KI), INTENT(IN)  :: PUREF     ! height of wind forcing                (m)
!
REAL, DIMENSION(KI), INTENT(IN)  :: PTA       ! air temperature forcing               (K)
REAL, DIMENSION(KI), INTENT(IN)  :: PQA       ! air humidity forcing                  (kg/m3)
REAL, DIMENSION(KI), INTENT(IN)  :: PRHOA     ! air density                           (kg/m3)
REAL, DIMENSION(KI,KSV),INTENT(IN) :: PSV     ! scalar variables
!                                             ! chemistry:   first char. in HSV: '#'  (molecule/m3)
!                                             !
 CHARACTER(LEN=6), DIMENSION(KSV),INTENT(IN):: HSV  ! name of all scalar variables
REAL, DIMENSION(KI), INTENT(IN)  :: PU        ! zonal wind                            (m/s)
REAL, DIMENSION(KI), INTENT(IN)  :: PV        ! meridian wind                         (m/s)
REAL, DIMENSION(KI,KSW),INTENT(IN) :: PDIR_SW ! direct  solar radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KI,KSW),INTENT(IN) :: PSCA_SW ! diffuse solar radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KSW),INTENT(IN)  :: PSW_BANDS ! mean wavelength of each shortwave band (m)
REAL, DIMENSION(KI), INTENT(IN)  :: PZENITH   ! zenithal angle       (radian from the vertical)
REAL, DIMENSION(KI), INTENT(IN)  :: PAZIM     ! azimuthal angle      (radian from North, clockwise)
REAL, DIMENSION(KI), INTENT(IN)  :: PLW       ! longwave radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KI), INTENT(IN)  :: PPS       ! pressure at atmospheric model surface (Pa)
REAL, DIMENSION(KI), INTENT(IN)  :: PPA       ! pressure at forcing level             (Pa)
REAL, DIMENSION(KI), INTENT(IN)  :: PZS       ! atmospheric model orography           (m)
REAL, DIMENSION(KI), INTENT(IN)  :: PCO2      ! CO2 concentration in the air          (kg/m3)
REAL, DIMENSION(KI), INTENT(IN)  :: PSN     ! snow precipitation                    (kg/m2/s)
REAL, DIMENSION(KI), INTENT(IN)  :: PRAIN     ! liquid precipitation                  (kg/m2/s)
!
!
REAL, DIMENSION(KI), INTENT(OUT) :: PSFTH     ! flux of heat                          (W/m2)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFTQ     ! flux of water vapor                   (kg/m2/s)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFU      ! zonal momentum flux                   (Pa)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFV      ! meridian momentum flux                (Pa)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFCO2    ! flux of CO2                           (kg/m2/s)
REAL, DIMENSION(KI,KSV),INTENT(OUT):: PSFTS   ! flux of scalar var.                   (kg/m2/s)
!
REAL, DIMENSION(KI), INTENT(OUT) :: PTRAD     ! radiative temperature                 (K)
REAL, DIMENSION(KI,KSW),INTENT(OUT):: PDIR_ALB! direct albedo for each spectral band  (-)
REAL, DIMENSION(KI,KSW),INTENT(OUT):: PSCA_ALB! diffuse albedo for each spectral band (-)
REAL, DIMENSION(KI), INTENT(OUT) :: PEMIS     ! emissivity                            (-)
!
REAL, DIMENSION(KI), INTENT(OUT) :: PTSURF    ! surface effective temperature         (K)
REAL, DIMENSION(KI), INTENT(OUT) :: PZ0       ! roughness length for momentum         (m)
REAL, DIMENSION(KI), INTENT(OUT) :: PZ0H      ! roughness length for heat             (m)
REAL, DIMENSION(KI), INTENT(OUT) :: PQSURF    ! specific humidity at surface          (kg/kg)
!
REAL, DIMENSION(KI), INTENT(IN) :: PPEW_A_COEF! implicit coefficients
REAL, DIMENSION(KI), INTENT(IN) :: PPEW_B_COEF! needed if HCOUPLING='I'
REAL, DIMENSION(KI), INTENT(IN) :: PPET_A_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPEQ_A_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPET_B_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPEQ_B_COEF
 CHARACTER(LEN=2),    INTENT(IN) :: HTEST ! must be equal to 'OK'
!
!
!*      0.2    declarations of local variables
!
INTEGER                     :: JSWB        ! loop counter on shortwave spectral bands
!         
REAL, DIMENSION(KI)  :: ZQA         ! specific humidity                 (kg/kg)
REAL, DIMENSION(KI)  :: ZEXNA       ! Exner function at forcing level
REAL, DIMENSION(KI)  :: ZEXNS       ! Exner function at surface level
REAL, DIMENSION(KI)  :: ZWIND       ! wind
!
! Ouput Diagnostics:
!
REAL, DIMENSION(KI)  :: ZU_CANYON   ! wind in canyon
REAL, DIMENSION(KI)  :: ZT_CANYON   ! temperature in canyon
REAL, DIMENSION(KI)  :: ZQ_CANYON   ! specific humidity in canyon
REAL, DIMENSION(KI)  :: ZAVG_T_CANYON ! temperature in canyon for town 
REAL, DIMENSION(KI)  :: ZAVG_Q_CANYON ! specific humidity in canyon for town
REAL, DIMENSION(KI)  :: ZT_CAN      ! temperature in canyon       (evolving in TEB)
REAL, DIMENSION(KI)  :: ZQ_CAN      ! specific humidity in canyon (evolving in TEB)
!
REAL, DIMENSION(KI)  :: ZPEW_A_COEF   ! implicit coefficients
REAL, DIMENSION(KI)  :: ZPEW_B_COEF   ! needed if HCOUPLING='I'
!
REAL, DIMENSION(KI) :: ZT_LOWCAN  ! temperature at lowest canyon level (K)
REAL, DIMENSION(KI) :: ZQ_LOWCAN  ! humidity    at lowest canyon level (kg/kg)
REAL, DIMENSION(KI) :: ZU_LOWCAN  ! wind        at lowest canyon level (m/s)
REAL, DIMENSION(KI) :: ZZ_LOWCAN  ! height      of lowest canyon level (m)
!
REAL, DIMENSION(KI) :: ZPEW_A_COEF_LOWCAN   ! implicit coefficients for wind coupling
REAL, DIMENSION(KI) :: ZPEW_B_COEF_LOWCAN   ! between first canopy level and road
!
REAL, DIMENSION(KI) :: ZTA        ! temperature at canyon level just above roof (K)
REAL, DIMENSION(KI) :: ZPA        ! pressure    at canyon level just above roof (K)
REAL, DIMENSION(KI) :: ZUA        ! wind        at canyon level just above roof (m/s)
REAL, DIMENSION(KI) :: ZUREF      ! height      of canyon level just above roof (m)
REAL, DIMENSION(KI) :: ZZREF      ! height      of canyon level just above roof (m)
!
REAL, DIMENSION(KI)  :: ZDIR_SW       ! total direct SW
REAL, DIMENSION(KI)  :: ZSCA_SW       ! total diffuse SW
REAL, DIMENSION(KI) :: ZAVG_SCA_SW
REAL, DIMENSION(KI) :: ZAVG_DIR_SW 
REAL, DIMENSION(KI,SIZE(PDIR_SW,2))  :: ZDIR_SWB ! total direct SW per band
REAL, DIMENSION(KI,SIZE(PSCA_SW,2))  :: ZSCA_SWB ! total diffuse SW per band
!
!
REAL, DIMENSION(KI)  :: ZLE_WL_A  ! latent heat flux on walls
REAL, DIMENSION(KI)  :: ZLE_WL_B  ! latent heat flux on walls
REAL, DIMENSION(KI)  :: ZAVG_H_WL
!
REAL, DIMENSION(KI)  :: ZPROD_BLD        ! averaged     energy production from solar panel (W/m2 bld)
REAL, DIMENSION(KI) :: ZHU_BLD
REAL, DIMENSION(KI) :: ZAVG_TI_BLD
REAL, DIMENSION(KI) :: ZAVG_QI_BLD
!
REAL, DIMENSION(KI)  :: ZRN_GRND    ! net radiation on ground built surf
REAL, DIMENSION(KI)  :: ZH_GRND     ! sensible heat flux on ground built surf
REAL, DIMENSION(KI)  :: ZLE_GRND    ! latent heat flux on ground built surf
REAL, DIMENSION(KI)  :: ZGFLX_GRND ! storage flux in ground built surf
REAL, DIMENSION(KI)  :: ZUW_GRND      ! momentum flux for ground built surf
REAL, DIMENSION(KI)  :: ZDUWDU_GRND   !
REAL, DIMENSION(KI)  :: ZAC_GRND      ! ground built surf aerodynamical conductance
REAL, DIMENSION(KI)  :: ZAC_GRND_WAT  ! ground built surf water aerodynamical conductance
REAL, DIMENSION(KI) :: ZEMIT_LW_GRND
REAL, DIMENSION(KI) :: ZREF_SW_GRND ! total solar rad reflected from ground
REAL, DIMENSION(KI)  :: ZAVG_UW_GRND
REAL, DIMENSION(KI)  :: ZAVG_DUWDU_GRND
REAL, DIMENSION(KI)  :: ZAVG_H_GRND
REAL, DIMENSION(KI)  :: ZAVG_AC_GRND
REAL, DIMENSION(KI)  :: ZAVG_AC_GRND_WAT
REAL, DIMENSION(KI)  :: ZAVG_E_GRND
REAL, DIMENSION(KI) :: ZAVG_REF_SW_GRND
REAL, DIMENSION(KI) :: ZAVG_EMIT_LW_GRND
!
REAL, DIMENSION(KI)  :: ZLEW_RF   ! latent heat flux on snowfree roof
REAL, DIMENSION(KI)  :: ZRNSN_RF  ! net radiation over snow
REAL, DIMENSION(KI)  :: ZHSN_RF   ! sensible heat flux over snow
REAL, DIMENSION(KI)  :: ZLESN_RF  ! latent heat flux over snow
REAL, DIMENSION(KI)  :: ZGSN_RF   ! flux under the snow
REAL, DIMENSION(KI)  :: ZMELT_RF    ! snow melt
REAL, DIMENSION(KI)  :: ZUW_RF      ! momentum flux for roofs
REAL, DIMENSION(KI)  :: ZDUWDU_RF   !
REAL, DIMENSION(KI)  :: ZAVG_UW_RF
REAL, DIMENSION(KI)  :: ZAVG_DUWDU_RF
REAL, DIMENSION(KI)  :: ZAVG_H_RF
REAL, DIMENSION(KI)  :: ZAVG_E_RF
!
REAL, DIMENSION(KI)  :: ZLEW_RD   ! latent heat flux on snowfree road
REAL, DIMENSION(KI)  :: ZRNSN_RD  ! net radiation over snow
REAL, DIMENSION(KI)  :: ZHSN_RD   ! sensible heat flux over snow
REAL, DIMENSION(KI)  :: ZLESN_RD  ! latent heat flux over snow
REAL, DIMENSION(KI)  :: ZGSN_RD   ! flux under the snow
REAL, DIMENSION(KI)  :: ZMELT_RD    ! snow melt
REAL, DIMENSION(KI)  :: ZAC_RD      ! road aerodynamical conductance
REAL, DIMENSION(KI)  :: ZAC_RD_WAT  ! road water aerodynamical conductance
!
REAL, DIMENSION(KI)  :: ZAC_GD    ! green area aerodynamical conductance
REAL, DIMENSION(KI)  :: ZAC_GD_WAT! green area water aerodynamical conductance
REAL, DIMENSION(KI,1):: ZESN_GD    ! green area snow emissivity
!
REAL, DIMENSION(KI)  :: ZAC_GRF ! green roof aerodynamical conductance
REAL, DIMENSION(KI)  :: ZAC_GRF_WAT! green roof water aerodynamical conductance
!
REAL, DIMENSION(KI)  :: ZTRAD         ! radiative temperature for current patch
REAL, DIMENSION(KI)  :: ZEMIS         ! emissivity for current patch
REAL, DIMENSION(KI,TOP%NTEB_PATCH) :: ZTRAD_PATCH ! radiative temperature for each patch
REAL, DIMENSION(KI,TOP%NTEB_PATCH) :: ZEMIS_PATCH ! emissivity for each patch
!
REAL, DIMENSION(KI)  :: ZDIR_ALB      ! direct albedo of town
REAL, DIMENSION(KI)  :: ZSCA_ALB      ! diffuse albedo of town
REAL, DIMENSION(KI,KSW,TOP%NTEB_PATCH) :: ZDIR_ALB_PATCH ! direct albedo per wavelength and patch
REAL, DIMENSION(KI,KSW,TOP%NTEB_PATCH) :: ZSCA_ALB_PATCH ! diffuse albedo per wavelength and patch
REAL, DIMENSION(KI)  :: ZAVG_DIR_ALB  ! direct albedo of town
REAL, DIMENSION(KI)  :: ZAVG_SCA_ALB  ! diffuse albedo of town
!
REAL, DIMENSION(KI)  :: ZSFCO2        ! CO2 flux over town
!
REAL, DIMENSION(KI)  :: ZRI           ! Richardson number
REAL, DIMENSION(KI)  :: ZCD           ! drag coefficient
REAL, DIMENSION(KI)  :: ZCDN          ! neutral drag coefficient
REAL, DIMENSION(KI)  :: ZCH           ! heat drag
REAL, DIMENSION(KI)  :: ZRN           ! net radiation over town
REAL, DIMENSION(KI)  :: ZH            ! sensible heat flux over town
REAL, DIMENSION(KI)  :: ZLE           ! latent heat flux over town
REAL, DIMENSION(KI)  :: ZGFLX        ! flux through the ground
REAL, DIMENSION(KI)  :: ZEVAP         ! evaporation (km/m2/s)
!
REAL, DIMENSION(KI)  :: ZAVG_CD       ! aggregated drag coefficient
REAL, DIMENSION(KI)  :: ZAVG_CDN      ! aggregated neutral drag coefficient
REAL, DIMENSION(KI)  :: ZAVG_RI       ! aggregated Richardson number
REAL, DIMENSION(KI)  :: ZAVG_CH       ! aggregated Heat transfer coefficient
!
REAL, DIMENSION(KI)  :: ZUSTAR        ! friction velocity
REAL, DIMENSION(KI)  :: ZSFU          ! momentum flux for patch (U direction)
REAL, DIMENSION(KI)  :: ZSFV          ! momentum flux for patch (V direction)
!
REAL, DIMENSION(KI)  :: ZH_TRAFFIC    ! anthropogenic sensible
!                                     ! heat fluxes due to traffic
REAL, DIMENSION(KI)  :: ZLE_TRAFFIC   ! anthropogenic latent
!                                     ! heat fluxes due to traffic
!
REAL                 :: ZBEGIN_TRAFFIC_TIME ! start traffic time (solar time, s)
REAL                 :: ZEND_TRAFFIC_TIME   ! end traffic time   (solar time, s)
!
REAL, DIMENSION(KI)  :: ZRESA    ! aerodynamical resistance
!
REAL, DIMENSION(KI) :: ZEMIT_LW_FAC
REAL, DIMENSION(KI) :: ZT_RAD_IND   ! Indoor mean radiant temperature [K]
REAL, DIMENSION(KI) :: ZREF_SW_FAC  ! total solar rad reflected from facade
!
REAL, DIMENSION(KI) :: ZAVG_Z0
REAL, DIMENSION(KI) :: ZAVG_RESA
REAL, DIMENSION(KI) :: ZAVG_USTAR        ! town avegared Ustar
REAL, DIMENSION(KI) :: ZAVG_BLD          ! town averaged building fraction
REAL, DIMENSION(KI) :: ZAVG_BLD_HEIGHT   ! town averaged building height
REAL, DIMENSION(KI) :: ZAVG_WL_O_HOR     ! town averaged Wall/hor ratio
REAL, DIMENSION(KI) :: ZAVG_CAN_HW_RATIO ! town averaged road aspect ratio
REAL, DIMENSION(KI) :: ZAVG_H
REAL, DIMENSION(KI) :: ZAVG_LE
REAL, DIMENSION(KI) :: ZAVG_RN
REAL, DIMENSION(KI) :: ZAVG_GFLX
REAL, DIMENSION(KI) :: ZAVG_REF_SW_FAC
REAL, DIMENSION(KI) :: ZAVG_EMIT_LW_FAC
REAL, DIMENSION(KI) :: ZAVG_T_RAD_IND
!
! absorbed solar and infra-red radiation by road, wall and roof
!
REAL, DIMENSION(KI) :: ZU_UTCI ! wind speed for the UTCI calculation (m/s) 

REAL, DIMENSION(KI) :: ZALFAU   ! V+(1) = alfa u'w'(1) + beta
REAL, DIMENSION(KI) :: ZBETAU   ! V+(1) = alfa u'w'(1) + beta
REAL, DIMENSION(KI) :: ZALFAT   ! Th+(1) = alfa w'th'(1) + beta
REAL, DIMENSION(KI) :: ZBETAT   ! Th+(1) = alfa w'th'(1) + beta
REAL, DIMENSION(KI) :: ZALFAQ   ! Q+(1) = alfa w'q'(1) + beta
REAL, DIMENSION(KI) :: ZBETAQ   ! Q+(1) = alfa w'q'(1) + beta
!***** CANOPY  *****
REAL, DIMENSION(KI) :: ZWAKE      ! reduction of average wind speed
!                                              ! in canyon due to direction average.

!new local variables for UTCI calculation
REAL, DIMENSION(KI) :: ZF1_o_B
!
!***** CANOPY  *****
REAL, DIMENSION(KI) :: ZSFLUX_U  ! Surface flux u'w' (m2/s2)
REAL, DIMENSION(KI) :: ZSFLUX_T  ! Surface flux w'T' (mK/s)
REAL, DIMENSION(KI) :: ZSFLUX_Q  ! Surface flux w'q' (kgm2/s)
REAL, DIMENSION(KI,SB%NLVL)   :: ZFORC_U   ! tendency due to drag force for wind
REAL, DIMENSION(KI,SB%NLVL)   :: ZDFORC_UDU! formal derivative of
!                                              ! tendency due to drag force for wind
REAL, DIMENSION(KI,SB%NLVL)   :: ZFORC_E   ! tendency due to drag force for TKE
REAL, DIMENSION(KI,SB%NLVL)   :: ZDFORC_EDE! formal derivative of
!                                              ! tendency due to drag force for TKE
REAL, DIMENSION(KI,SB%NLVL)   :: ZFORC_T   ! tendency due to drag force for Temp
REAL, DIMENSION(KI,SB%NLVL)   :: ZDFORC_TDT! formal derivative of
!                                              ! tendency due to drag force for Temp
REAL, DIMENSION(KI,SB%NLVL)   :: ZFORC_Q   ! tendency due to drag force for hum
REAL, DIMENSION(KI,SB%NLVL)   :: ZDFORC_QDQ! formal derivative of
!                                           ! tendency due to drag force for hum.
REAL, DIMENSION(KI) :: ZLAMBDA_F  ! frontal density (-)
REAL, DIMENSION(KI) :: ZLMO       ! Monin-Obukhov length at canopy height (m)
REAL, DIMENSION(KI,SB%NLVL)   :: ZL         ! Mixing length generic profile at mid levels
!
REAL, DIMENSION(KI) :: ZCOEF
!
REAL                       :: ZCONVERTFACM0_SLT, ZCONVERTFACM0_DST
REAL                       :: ZCONVERTFACM3_SLT, ZCONVERTFACM3_DST
REAL                       :: ZCONVERTFACM6_SLT, ZCONVERTFACM6_DST
!
INTEGER                           :: JI
INTEGER                           :: JLAYER
INTEGER                           :: JJ
!
! number of TEB patches
!
INTEGER                    :: JP, IBEG, IEND ! loop counter
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------------
! Preliminaries:
!-------------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('COUPLING_TEB_N',0,zhook_handle)
IF (htest/='OK') THEN
  CALL abor1_sfx('COUPLING_TEBN: FATAL ERROR DURING ARGUMENT TRANSFER')
END IF

!-------------------------------------------------------------------------------------
!
! scalar fluxes
!
psfts(:,:) = 0.
!
! broadband radiative fluxes
!
zdir_sw(:) = 0.
zsca_sw(:) = 0.
DO jswb=1,ksw
  !add directionnal contrib from scattered radiation
  CALL circumsolar_rad(pdir_sw(:,jswb), psca_sw(:,jswb), pzenith, zf1_o_b)
  zdir_swb(:,jswb) = pdir_sw(:,jswb) + psca_sw(:,jswb) * zf1_o_b
  zsca_swb(:,jswb) = psca_sw(:,jswb) * (1. - zf1_o_b)
  !add directionnal contrib from scattered radiation
  DO jj=1,SIZE(pdir_sw,1)
    zdir_sw(jj) = zdir_sw(jj) + zdir_swb(jj,jswb)
    zsca_sw(jj) = zsca_sw(jj) + zsca_swb(jj,jswb)
  ENDDO
END DO
!
DO jj=1,ki
! specific humidity (conversion from kg/m3 to kg/kg)
!
  zqa(jj) = pqa(jj) / prhoa(jj)
!
! wind
!
  zwind(jj) = sqrt(pu(jj)**2+pv(jj)**2)
!
ENDDO
! method of wind coupling
!
IF (hcoupling=='I') THEN
  zpew_a_coef = ppew_a_coef
  zpew_b_coef = ppew_b_coef
ELSE
  zpew_a_coef =  0.
  zpew_b_coef =  zwind
END IF
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Time evolution
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
top%TTIME%TIME = top%TTIME%TIME + ptstep
 CALL add_forecast_to_date_surf(top%TTIME%TDATE%YEAR, top%TTIME%TDATE%MONTH,&
                                top%TTIME%TDATE%DAY, top%TTIME%TIME)
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  Anthropogenic fluxes (except building heating)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
zbegin_traffic_time = 21600.
zend_traffic_time   = 64800.
!
WHERE( ptsun>zbegin_traffic_time  .AND.  ptsun<zend_traffic_time  )
  zh_traffic(:) = nt%AL(1)%XH_TRAFFIC   (:)
  zle_traffic(:) = nt%AL(1)%XLE_TRAFFIC  (:)
ELSEWHERE
  zh_traffic(:) = 0.
  zle_traffic(:) = 0.   
END WHERE
!
!--------------------------------------------------------------------------------------
!  Canyon forcing for TEB
!--------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------
! Town averaged quantities to force canopy atmospheric layers
!-------------------------------------------------------------------------------------

DO jp=1,top%NTEB_PATCH
  CALL add_patch_contrib(jp,zavg_bld,         nt%AL(jp)%XBLD         )
  CALL add_patch_contrib(jp,zavg_bld_height,  nt%AL(jp)%XBLD_HEIGHT  )
  CALL add_patch_contrib(jp,zavg_wl_o_hor,    nt%AL(jp)%XWALL_O_HOR  )
  CALL add_patch_contrib(jp,zavg_can_hw_ratio,nt%AL(jp)%XCAN_HW_RATIO)
  CALL add_patch_contrib(jp,zavg_z0,          nt%AL(jp)%XZ0_TOWN     )
END DO
!
IF (top%LCANOPY) THEN
  !-------------------------------------------------------------------------------------
  ! Updates canopy vertical grid as a function of forcing height
  !-------------------------------------------------------------------------------------
  !
  !* determines where is the forcing level and modifies the upper levels of the canopy grid
  !
  CALL canopy_grid_update(ki, zavg_bld_height, zavg_bld_height+puref, sb)
  !
  !* Initialisations of T, Q, TKE and wind at first time step
  !
  IF(any(sb%XT(:,:) == xundef)) THEN
    DO jlayer=1,sb%NLVL
      sb%XT(:,jlayer) = pta(:)
      sb%XQ(:,jlayer) = pqa(:)
      sb%XU(:,jlayer) = 2./xpi * zwind(:)                                  &
              * log( (          2.* nt%AL(1)%XBLD_HEIGHT(:)/3.) / nt%AL(1)%XZ0_TOWN(:))   &
              / log( (puref(:)+ 2.* nt%AL(1)%XBLD_HEIGHT(:)/3.) / nt%AL(1)%XZ0_TOWN(:))
    END  DO
    sb%XTKE(:,:) = 1.
  ENDIF
  !
  !* default forcing above roof: forcing level
  zuref(:) = puref(:)
  zzref(:) = pzref(:)
  zua(:)   = sb%XU(:,sb%NLVL)
  zta(:)   = sb%XT(:,sb%NLVL)
  zqa(:)   = sb%XQ(:,sb%NLVL)/prhoa(:)
  zpa(:)   = sb%XP(:,sb%NLVL)
  !* for the time being, only one value is kept for wall in-canyon forcing, in the middle of the canyon
  zu_canyon(:) = zua(:)
  zt_canyon(:) = zta(:)
  zq_canyon(:) = zqa(:)
  DO jlayer=1,sb%NLVL-1
    DO ji=1,ki
      !* finds middle canyon layer
      IF (sb%XZ(ji,jlayer)<zavg_bld_height(ji)/2. .AND. sb%XZ(ji,jlayer+1)>=zavg_bld_height(ji)/2.) THEN
        zcoef(ji) = (zavg_bld_height(ji)/2.-sb%XZ(ji,jlayer))/(sb%XZ(ji,jlayer+1)-sb%XZ(ji,jlayer))
        zu_canyon(ji) = sb%XU(ji,jlayer) + zcoef(ji) * (sb%XU(ji,jlayer+1)-sb%XU(ji,jlayer))
        zt_canyon(ji) = sb%XT(ji,jlayer) + zcoef(ji) * (sb%XT(ji,jlayer+1)-sb%XT(ji,jlayer))
        zq_canyon(ji) =(sb%XQ(ji,jlayer) + zcoef(ji) * (sb%XQ(ji,jlayer+1)-sb%XQ(ji,jlayer)))/prhoa(ji)
      END IF
      !* finds layer just above roof (at least 1m above roof)
      IF (sb%XZ(ji,jlayer)<zavg_bld_height(ji)+1. .AND. sb%XZ(ji,jlayer+1)>=zavg_bld_height(ji)+1.) THEN
        zuref(ji) = sb%XZ(ji,jlayer+1) - zavg_bld_height(ji)
        zzref(ji) = sb%XZ(ji,jlayer+1) - zavg_bld_height(ji)
        zta(ji) = sb%XT(ji,jlayer+1)
        zqa(ji) = sb%XQ(ji,jlayer+1)/prhoa(ji)
        !ZUA  (JI) = XU(JI,JLAYER+1)
        zua(ji) = max(sb%XU(ji,jlayer+1) - 2.*sqrt(sb%XTKE(ji,jlayer+1)) , sb%XU(ji,jlayer+1)/3.)
        zpa(ji) = sb%XP(ji,jlayer+1)
        zlmo(ji) = sb%XLMO(ji,jlayer+1)
      END IF
    END DO
  END DO
  zu_canyon= max(zu_canyon,0.2)
  zu_lowcan=sb%XU(:,1)
  zt_lowcan=sb%XT(:,1)
  zq_lowcan=sb%XQ(:,1) / prhoa(:)
  zz_lowcan=sb%XZ(:,1)
  WHERE(zpa==xundef) zpa = ppa   ! security for first time step
  !
  !-------------------------------------------------------------------------------------
  ! determine the vertical profile for mixing and dissipative lengths (at full levels)
  !-------------------------------------------------------------------------------------
  !
  ! frontal density
  zlambda_f(:) = zavg_can_hw_ratio*zavg_bld / (0.5*xpi)
  !
  CALL sm10(sb%XZ, zavg_bld_height, zlambda_f, zl)
  !
  !-------------------------------------------------------------------------------------
  ! computes coefficients for implicitation
  !-------------------------------------------------------------------------------------
  !
  zavg_uw_grnd(:)    = 0.
  zavg_duwdu_grnd(:) = 0.
  zavg_uw_rf(:)      = 0.
  zavg_duwdu_rf(:)   = 0.
  zavg_h_grnd(:)     = 0.
  zavg_h_wl(:)       = 0.
  zavg_h_rf(:)       = 0.
  zavg_e_grnd(:)     = 0.
  zavg_e_rf(:)       = 0.
  zavg_ac_grnd(:)    = 0.
  zavg_ac_grnd_wat(:)= 0.
  zsflux_u(:)        = 0.
  zsflux_t(:)        = 0.
  zsflux_q(:)        = 0.
  !
  DO jlayer=1,sb%NLVL-1
    !* Monin-Obuhkov theory not used inside the urban canopy
    ! => neutral mixing  if layer is below : (roof level +1 meter)
    WHERE (sb%XZ(:,jlayer)<=zavg_bld_height(:)+1.) sb%XLMO(:,jlayer) = xundef
  ENDDO
  !
  !* computes tendencies on wind and Tke due to canopy
  CALL teb_canopy(ki, sb, zavg_bld, zavg_bld_height, zavg_wl_o_hor, ppa, prhoa, &
                  zavg_duwdu_grnd, zavg_uw_rf, zavg_duwdu_rf, zavg_h_wl,         &
                  zavg_h_rf, zavg_e_rf, zavg_ac_grnd, zavg_ac_grnd_wat, zforc_u, &
                  zdforc_udu, zforc_e, zdforc_ede, zforc_t, zdforc_tdt, zforc_q, &
                  zdforc_qdq )
  !
  !* computes coefficients for implicitation
  CALL canopy_evol(sb, ki, ptstep, 1, zl, zwind, pta, pqa, ppa, prhoa, &
                   zsflux_u, zsflux_t, zsflux_q, zforc_u, zdforc_udu,   &
                   zforc_e, zdforc_ede, zforc_t, zdforc_tdt, zforc_q,   &
                   zdforc_qdq, sb%XLM, sb%XLEPS, zavg_ustar, zalfau,  &
                   zbetau, zalfat, zbetat, zalfaq, zbetaq)
  !
  zpew_a_coef_lowcan = - zalfau / prhoa
  zpew_b_coef_lowcan = zbetau  
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
ELSE              ! no canopy case
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  DO ji=1,ki
    !* skimming flow for h/w>1 (maximum effect of direction on wind in the canyon);
    !* isolated flow for h/w<0.5 (wind is the same in large streets for all dir.)
    !* wake flow between.
    !
    zwake(ji)= 1. + (2./xpi-1.) * 2. * (zavg_can_hw_ratio(ji)-0.5)
    zwake(ji)= max(min(zwake(ji),1.),2./xpi)
    !
    !* Estimation of canyon wind speed from wind just above roof level
    !  (at 1.33h). Wind at 1.33h is estimated using the log law.
    !
    IF (zavg_bld_height(ji) .GT. 0.) THEN
      zu_canyon(ji) = zwake(ji) * exp(-zavg_can_hw_ratio(ji)/4.) * zwind(ji)     &
                  * log( (           2.* zavg_bld_height(ji)/3.) / zavg_z0(ji))   &
                  / log( (puref(ji)+ 2.* zavg_bld_height(ji)/3.) / zavg_z0(ji))
      zz_lowcan(ji) = zavg_bld_height(ji) / 2.
    ELSE
      zu_canyon(ji) = zwind(ji)
      zz_lowcan(ji) = pzref(ji)
    ENDIF
  END DO
  !
  !* Without SBL scheme, canyon air is assumed at mid height
  zu_lowcan = zu_canyon

  zt_lowcan = nt%AL(1)%XT_CANYON
  zq_lowcan = nt%AL(1)%XQ_CANYON
  zt_canyon = nt%AL(1)%XT_CANYON
  zq_canyon = nt%AL(1)%XQ_CANYON

  zuref     = puref
  zzref     = pzref
  zta       = pta
  zua       = zwind
  zpa       = ppa
  zpew_a_coef_lowcan =  0.
  zpew_b_coef_lowcan =  zu_canyon

END IF
!
! Exner functions
!
zexns(:) = (pps(:)/xp00)**(xrd/xcpd)
zexna(:) = (zpa(:)/xp00)**(xrd/xcpd)

!--------------------------------------------------------------------------------------
! Over Urban surfaces/towns:
!--------------------------------------------------------------------------------------
!
DO jp = 1,top%NTEB_PATCH
  !
  zt_can = zt_canyon
  zq_can = zq_canyon
  !
  IF (top%LCANOPY) THEN
    nt%AL(jp)%XT_CANYON(:) = zt_canyon(:)
    nt%AL(jp)%XQ_CANYON(:) = zq_canyon(:)
  END IF
  !
  zlesn_rf(:) = 0.
  zlesn_rd(:) = 0.
  td%NDMT%AL(jp)%XG_GREENROOF_ROOF(:) = 0.
  !
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! Call the physical routines of TEB (including gardens & greenroofs)
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
  CALL teb_garden(dtco, g, top, nt%AL(jp), bop, nb%AL(jp), tpn, tir, td%NDMT%AL(jp), gdm, grm, jp, &
                  cimplicit_wind, ptsun, zt_can, zq_can, zu_canyon, zt_lowcan, zq_lowcan, &
                  zu_lowcan, zz_lowcan, zpew_a_coef, zpew_b_coef, zpew_a_coef_lowcan,     &
                  zpew_b_coef_lowcan, pps, zpa, zexns, zexna, zta, zqa, prhoa, pco2, plw, &
                  zdir_swb, zsca_swb, psw_bands, ksw, pzenith, pazim, prain, psn, zzref,  &
                  zuref, zua, zh_traffic, zle_traffic, ptstep, zlew_rf, zlew_rd, zle_wl_a,&
                  zle_wl_b, zrnsn_rf, zhsn_rf, zlesn_rf, zgsn_rf, zmelt_rf, zrnsn_rd,     &
                  zhsn_rd, zlesn_rd, zgsn_rd, zmelt_rd, zrn_grnd, zh_grnd, zle_grnd,      &
                  zgflx_grnd, zrn, zh, zle, zgflx, zevap, zsfco2, zuw_grnd,               &
                  zuw_rf, zduwdu_grnd, zduwdu_rf, zustar, zcd, zcdn, zch, zri, ztrad,     &
                  zemis, zdir_alb, zsca_alb, zresa, zac_rd, zac_gd, zac_grf, zac_rd_wat,  &
                  zac_gd_wat, zac_grf_wat, kday, zemit_lw_fac, zemit_lw_grnd, zt_rad_ind, &
                  zref_sw_grnd, zref_sw_fac, zhu_bld, ptime, zprod_bld       )


  !
  IF (.NOT. top%LCANOPY) THEN

    CALL add_patch_contrib(jp,zavg_t_canyon,zt_can)
    CALL add_patch_contrib(jp,zavg_q_canyon,zq_can)
    !
    ! Momentum fluxes
    !
    zsfu = 0.
    zsfv = 0.
    DO jj=1,SIZE(pu)
      IF (zwind(jj)>0.) THEN
        zcoef(jj) = - prhoa(jj) * zustar(jj)**2 / zwind(jj)
        zsfu(jj) = zcoef(jj) * pu(jj)
        zsfv(jj) = zcoef(jj) * pv(jj)
      ENDIF
    ENDDO
    CALL add_patch_contrib(jp,psfu,zsfu)
    CALL add_patch_contrib(jp,psfv,zsfv)
    !
  ENDIF
  !
  !-------------------------------------------------------------------------------------
  ! Outputs:
  !-------------------------------------------------------------------------------------
  !
  ! Grid box average fluxes/properties: Arguments and standard diagnostics
  !
  CALL add_patch_contrib(jp,psfth,zh)
  CALL add_patch_contrib(jp,psftq,zevap)
  CALL add_patch_contrib(jp,psfco2,zsfco2)
  !
  !
  ! Albedo for each wavelength and patch
  !
  DO jswb=1,SIZE(psw_bands)
    DO jj=1,SIZE(zdir_alb)
      zdir_alb_patch(jj,jswb,jp) = zdir_alb(jj)
      zsca_alb_patch(jj,jswb,jp) = zsca_alb(jj)
    ENDDO
  END DO
  !
  ! emissivity and radiative temperature
  !
  zemis_patch(:,jp) = zemis
  ztrad_patch(:,jp) = ztrad
  !
  ! computes some aggregated diagnostics
  !
  CALL add_patch_contrib(jp,zavg_cd ,zcd )
  CALL add_patch_contrib(jp,zavg_cdn,zcdn)
  CALL add_patch_contrib(jp,zavg_ri ,zri )
  CALL add_patch_contrib(jp,zavg_ch ,zch )
  CALL add_patch_contrib(jp,zavg_rn ,zrn )
  CALL add_patch_contrib(jp,zavg_h  ,zh  )
  CALL add_patch_contrib(jp,zavg_le ,zle )
  CALL add_patch_contrib(jp,zavg_gflx ,zgflx )
  !
  !* warning: aerodynamical resistance does not yet take into account gardens
  CALL add_patch_contrib(jp,zavg_resa,1./zresa)
  IF (jp==top%NTEB_PATCH) zavg_resa = 1./zavg_resa
  !
  !-------------------------------------------------------------------------------------
  ! Diagnostics on each patch
  !-------------------------------------------------------------------------------------
  !
  IF (td%MTO%LSURF_MISC_BUDGET) THEN
    !
    ! cumulated diagnostics 
    ! ---------------------
    !
    CALL cumul_diag_teb_n(td%NDMTC%AL(jp), td%NDMT%AL(jp), gdm%VD%NDEC%AL(jp), gdm%VD%NDE%AL(jp), &
                          grm%VD%NDEC%AL(jp), grm%VD%NDE%AL(jp), top, ptstep)
    !
  END IF
  !
  !
  !-------------------------------------------------------------------------------------
  ! Computes averaged parameters necessary for UTCI
  !-------------------------------------------------------------------------------------
  !
  IF (td%O%N2M >0 .AND. td%DUT%LUTCI) THEN
    CALL add_patch_contrib(jp,zavg_ref_sw_grnd ,zref_sw_grnd )
    CALL add_patch_contrib(jp,zavg_ref_sw_fac  ,zref_sw_fac  )
    CALL add_patch_contrib(jp,zavg_sca_sw      ,zsca_sw      )
    CALL add_patch_contrib(jp,zavg_dir_sw      ,zdir_sw      )
    CALL add_patch_contrib(jp,zavg_emit_lw_fac ,zemit_lw_fac )
    CALL add_patch_contrib(jp,zavg_emit_lw_grnd,zemit_lw_grnd)
    CALL add_patch_contrib(jp,zavg_t_rad_ind   ,zt_rad_ind   )
    CALL add_patch_contrib(jp,zavg_ti_bld      ,nb%AL(jp)%XTI_BLD)
    CALL add_patch_contrib(jp,zavg_qi_bld      ,nb%AL(jp)%XQI_BLD)
  END IF
  !
  !-------------------------------------------------------------------------------------
  ! Use of the canopy version of TEB
  !-------------------------------------------------------------------------------------
  !
  IF (top%LCANOPY) THEN
    !-------------------------------------------------------------------------------------
    ! Town averaged quantities to force canopy atmospheric layers
    !-------------------------------------------------------------------------------------

    CALL add_patch_contrib(jp,zavg_duwdu_grnd, zduwdu_grnd )
    CALL add_patch_contrib(jp,zavg_uw_rf     , zuw_rf)
    CALL add_patch_contrib(jp,zavg_duwdu_rf  , zduwdu_rf)
    CALL add_patch_contrib(jp,zavg_h_wl      , 0.5*(td%NDMT%AL(jp)%XH_WALL_A+td%NDMT%AL(jp)%XH_WALL_B))
    CALL add_patch_contrib(jp,zavg_h_rf      , (td%NDMT%AL(jp)%XH_ROOF + nt%AL(jp)%XH_INDUSTRY))
    CALL add_patch_contrib(jp,zavg_e_rf      , (td%NDMT%AL(jp)%XLE_ROOF+ nt%AL(jp)%XLE_INDUSTRY)/xlvtt)
    !
    !-------------------------------------------------------------------------------------
    ! Computes the impact of canopy and surfaces on air
    !-------------------------------------------------------------------------------------
    !
    zac_grnd(:) = (nt%AL(jp)%XROAD(:)*zac_rd(:) + nt%AL(jp)%XGARDEN(:)*zac_gd(:)) / &
                      (nt%AL(jp)%XROAD(:)+nt%AL(jp)%XGARDEN(:))
    zac_grnd_wat(:) = (nt%AL(jp)%XROAD(:)*zac_rd_wat(:) + nt%AL(jp)%XGARDEN(:)*zac_gd_wat(:)) / &
                      (nt%AL(jp)%XROAD(:)+nt%AL(jp)%XGARDEN(:))
    !
    CALL add_patch_contrib(jp,zavg_ac_grnd    , zac_grnd    )
    CALL add_patch_contrib(jp,zavg_ac_grnd_wat, zac_grnd_wat)
    CALL add_patch_contrib(jp,zsflux_u        , zuw_grnd * (1.-nt%AL(jp)%XBLD))
    CALL add_patch_contrib(jp,zsflux_t        , zh_grnd  * (1.-nt%AL(jp)%XBLD)/xcpd/prhoa)
    CALL add_patch_contrib(jp,zsflux_q        , zle_grnd * (1.-nt%AL(jp)%XBLD)/xlvtt)
    !
  END IF
  !
  !-------------------------------------------------------------------------------------
  ! end of loop on TEB patches
END DO
!-------------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------------
!* Evolution of canopy air if canopy option is active
!-------------------------------------------------------------------------------------
!
IF (top%LCANOPY) THEN
  !
  !-------------------------------------------------------------------------------------
  !* Impact of TEB fluxes on the air
  !-------------------------------------------------------------------------------------
  !
  CALL teb_canopy(ki, sb, zavg_bld, zavg_bld_height, zavg_wl_o_hor, ppa, prhoa, &
                  zavg_duwdu_grnd, zavg_uw_rf, zavg_duwdu_rf, zavg_h_wl,         &
                  zavg_h_rf, zavg_e_rf, zavg_ac_grnd, zavg_ac_grnd_wat, zforc_u, &
                  zdforc_udu, zforc_e, zdforc_ede, zforc_t, zdforc_tdt, zforc_q, &
                  zdforc_qdq )
  !
  !-------------------------------------------------------------------------------------
  !* Evolution of canopy air due to these impacts
  !-------------------------------------------------------------------------------------
  !
  CALL canopy_evol(sb, ki, ptstep, 2, zl, zwind, pta, pqa, ppa, prhoa,  &
                   zsflux_u, zsflux_t, zsflux_q, zforc_u, zdforc_udu,    &
                   zforc_e, zdforc_ede, zforc_t, zdforc_tdt, zforc_q,    &
                   zdforc_qdq, sb%XLM, sb%XLEPS, zavg_ustar, zalfau,   &
                   zbetau, zalfat, zbetat, zalfaq, zbetaq      )
  !
  !-------------------------------------------------------------------------------------
  ! Momentum fluxes in the case canopy is active
  !-------------------------------------------------------------------------------------
  !
  psfu=0.
  psfv=0.
  zavg_z0(:) = min(zavg_z0(:),puref(:)*0.5)
  zavg_cdn=(xkarman/log(puref(:)/zavg_z0(:)))**2
  zavg_cd = zavg_cdn
  zavg_ri = 0.
  DO jj=1,SIZE(pu)
    IF (zwind(jj)>0.) THEN
      zcoef(jj) = - prhoa(jj) * zavg_ustar(jj)**2 / zwind(jj)
      psfu(jj) = zcoef(jj) * pu(jj)
      psfv(jj) = zcoef(jj) * pv(jj)
      zavg_cd(jj) = zavg_ustar(jj)**2 / zwind(jj)**2
      zavg_ri(jj) = -xg/pta(jj)*zsflux_t(jj)/zavg_ustar(jj)**4
    ENDIF
  ENDDO
  !
  !-------------------------------------------------------------------------------------
  ! End of specific case with canopy option
  !-------------------------------------------------------------------------------------
  !
END IF
!
!-------------------------------------------------------------------------------------
! Outputs:
!-------------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------------
!Radiative properties should be at time t+1 (see by the atmosphere) in order to close
!the energy budget between surfex and the atmosphere. It is not the case here
!for ALB and EMIS
!-------------------------------------------------------------------------------------
!
 CALL average_rad(top%XTEB_PATCH, zdir_alb_patch, zsca_alb_patch, zemis_patch, &
                  ztrad_patch, pdir_alb, psca_alb, pemis, ptrad )
!
!-------------------------------------------------------------------------------
!Physical properties see by the atmosphere in order to close the energy budget 
!between surfex and the atmosphere. All variables should be at t+1 but very 
!difficult to do. Maybe it will be done later. However, Ts can be at time t+1
!-------------------------------------------------------------------------------
!
ptsurf(:) = ptrad(:) ! Should be the surface effective temperature; not radative
pz0(:) = zavg_z0(:) ! Should account for ISBA (greenroof and garden) Z0
pz0h(:) = pz0(:) / 200.    ! Should account for ISBA (greenroof and garden) Z0
pqsurf(:) = nt%AL(1)%XQ_CANYON(:) ! Should account for ISBA (greenroof and garden) Qs
!
!-------------------------------------------------------------------------------------
! Scalar fluxes:
!-------------------------------------------------------------------------------------
!
zavg_ustar(:) = sqrt(sqrt(psfu**2+psfv**2))
!
!
IF (cht%SVT%NBEQ>0) THEN

  ibeg = cht%SVT%NSV_CHSBEG
  iend = cht%SVT%NSV_CHSEND

  IF (cht%CCH_DRY_DEP == "WES89") THEN
    CALL ch_dep_town(zavg_resa,  zavg_ustar, pta, ptrad, zavg_wl_o_hor,&
                     psv(:,ibeg:iend), cht%SVT%CSV(ibeg:iend), cht%XDEP(:,1:cht%SVT%NBEQ)  )
   
    DO ji=ibeg,iend
!cdir nodep
      DO jj=1,SIZE(psfts,1)
        psfts(jj,ji) = - psv(jj,ji) * cht%XDEP(jj,ji-ibeg+1)
      ENDDO
    ENDDO

    IF (cht%SVT%NAEREQ > 0 ) THEN

      ibeg = cht%SVT%NSV_AERBEG
      iend = cht%SVT%NSV_AEREND

      CALL ch_aer_dep(psv(:,ibeg:iend), psfts(:,ibeg:iend), &
                      zavg_ustar, zavg_resa, pta, prhoa)   
    END IF

  ELSE

    ibeg = cht%SVT%NSV_CHSBEG
    iend = cht%SVT%NSV_CHSEND

    DO ji=ibeg,iend
      psfts(:,ji) =0.
    ENDDO

    ibeg = cht%SVT%NSV_AERBEG
    iend = cht%SVT%NSV_AEREND

    IF(ibeg.LT.iend) THEN
      DO ji=ibeg,iend
        psfts(:,ji) =0.
      ENDDO
    ENDIF
  ENDIF

ENDIF

IF (cht%SVT%NDSTEQ>0) THEN
  ! Blindage à enlever lorsque que TEB aura été corrigé
  zustar(:) = min(zustar(:), 10.)
  zresa(:) = max(zresa(:), 10.)
  !
  ibeg = cht%SVT%NSV_DSTBEG
  iend = cht%SVT%NSV_DSTEND
  !
  CALL dslt_dep(psv(:,ibeg:iend), psfts(:,ibeg:iend), zustar, zresa, pta, prhoa, &
                dst%XEMISSIG_DST, dst%XEMISRADIUS_DST, jpmode_dst, xdensity_dst, &
                xmolarweight_dst, zconvertfacm0_dst, zconvertfacm6_dst,          &
                zconvertfacm3_dst, lvarsig_dst, lrgfix_dst, cvermod  )  

  CALL massflux2momentflux(         &
    psfts(:,ibeg:iend),             & !I/O ![kg/m2/sec] In: flux of only mass, out: flux of moments
    prhoa,                          & !I [kg/m3] air density
    dst%XEMISRADIUS_DST,                &!I [um] emitted radius for the modes (max 3)
    dst%XEMISSIG_DST,                   &!I [-] emitted sigma for the different modes (max 3)
    ndstmde,                        &
    zconvertfacm0_dst,              &
    zconvertfacm6_dst,              &
    zconvertfacm3_dst,              &
    lvarsig_dst, lrgfix_dst         )  
ENDIF
IF (cht%SVT%NSLTEQ>0) THEN
  !
  ibeg = cht%SVT%NSV_SLTBEG
  iend = cht%SVT%NSV_SLTEND
  !
  CALL dslt_dep(psv(:,ibeg:iend), psfts(:,ibeg:iend), zustar, zresa, pta, prhoa, &
                slt%XEMISSIG_SLT, slt%XEMISRADIUS_SLT, jpmode_slt, xdensity_slt, &
                xmolarweight_slt, zconvertfacm0_slt, zconvertfacm6_slt,          &
                zconvertfacm3_slt, lvarsig_slt, lrgfix_slt, cvermod  )  

  CALL massflux2momentflux(         &
    psfts(:,ibeg:iend),             & !I/O ![kg/m2/sec] In: flux of only mass, out: flux of moments
    prhoa,                          & !I [kg/m3] air density
    slt%XEMISRADIUS_SLT,                &!I [um] emitted radius for the modes (max 3)
    slt%XEMISSIG_SLT,                   &!I [-] emitted sigma for the different modes (max 3)
    nsltmde,                        &
    zconvertfacm0_slt,              &
    zconvertfacm6_slt,              &
    zconvertfacm3_slt,              &
    lvarsig_slt, lrgfix_slt         ) 
ENDIF
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Inline diagnostics
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
 CALL diag_inline_teb_n(td%O, td%D, sb, nt%AL(1), top%LCANOPY, &
                        pta, ptrad, zqa, ppa, pps, prhoa, pu, pv, zwind, pzref, puref, &
                        zavg_cd, zavg_cdn, zavg_ri, zavg_ch, zavg_z0, ptrad, pemis,    &
                        pdir_alb, psca_alb, plw, zdir_swb, zsca_swb,  psfth, psftq,    &
                        psfu, psfv, psfco2, zavg_rn, zavg_h, zavg_le, zavg_gflx       )
!
!-------------------------------------------------------------------------------------
! Stores Canyon air and humidity if historical option of TEB is active
!-------------------------------------------------------------------------------------
!
IF (.NOT. top%LCANOPY) THEN
  DO jp=1,top%NTEB_PATCH
    nt%AL(jp)%XT_CANYON(:) = zavg_t_canyon(:)
    nt%AL(jp)%XQ_CANYON(:) = zavg_q_canyon(:)
  END DO
END IF
!          
!-------------------------------------------------------------------------------------
! Thermal confort index
!-------------------------------------------------------------------------------------
!
IF (td%DUT%LUTCI .AND. td%O%N2M >0) THEN
  DO jj=1,ki
    IF (td%D%XZON10M(jj)/=xundef) THEN
      zu_utci(jj) = sqrt(td%D%XZON10M(jj)**2+td%D%XMER10M(jj)**2)
    ELSE
      zu_utci(jj) = zwind(jj)
    ENDIF
  ENDDO
 CALL utci_teb(nt%AL(1), td%DUT, zavg_ti_bld, zavg_qi_bld, zu_utci, pps, zavg_ref_sw_grnd, &
               zavg_ref_sw_fac, zavg_sca_sw, zavg_dir_sw, pzenith, zavg_emit_lw_fac,  &
               zavg_emit_lw_grnd, plw, zavg_t_rad_ind    )
 CALL utcic_stress(ptstep,td%DUT%XUTCI_IN      ,td%DUT%XUTCIC_IN      )
 CALL utcic_stress(ptstep,td%DUT%XUTCI_OUTSUN  ,td%DUT%XUTCIC_OUTSUN  )
 CALL utcic_stress(ptstep,td%DUT%XUTCI_OUTSHADE,td%DUT%XUTCIC_OUTSHADE)
ELSE IF (td%DUT%LUTCI) THEN
  td%DUT%XUTCI_IN         (:) = xundef
  td%DUT%XUTCI_OUTSUN     (:) = xundef
  td%DUT%XUTCI_OUTSHADE   (:) = xundef
  td%DUT%XTRAD_SUN        (:) = xundef
  td%DUT%XTRAD_SHADE      (:) = xundef
  td%DUT%XUTCIC_IN      (:,:) = xundef
  td%DUT%XUTCIC_OUTSUN  (:,:) = xundef
  td%DUT%XUTCIC_OUTSHADE(:,:) = xundef
ENDIF

!
IF (lhook) CALL dr_hook('COUPLING_TEB_N',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
CONTAINS
SUBROUTINE add_patch_contrib(JP,PAVG,PFIELD)
INTEGER, INTENT(IN) :: JP
REAL, DIMENSION(:), INTENT(INOUT) :: PAVG
REAL, DIMENSION(:), INTENT(IN)    :: PFIELD
!
IF (jp==1) pavg = 0.
pavg = pavg + top%XTEB_PATCH(:,jp) * pfield(:)
!
END SUBROUTINE add_patch_contrib
!-------------------------------------------------------------------------------------
!
END SUBROUTINE coupling_teb_n