teb.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 teb  (TOP, T, BOP, B, TIR, DMT, HIMPLICIT_WIND, PTSUN,                   &
                     PT_CANYON, PQ_CANYON, PU_CANYON, PT_LOWCAN, PQ_LOWCAN, PU_LOWCAN,  &
                     PZ_LOWCAN, PPEW_A_COEF, PPEW_B_COEF, PPEW_A_COEF_LOWCAN,           &
                     PPEW_B_COEF_LOWCAN, PPS, PPA, PEXNS, PEXNA, PTA, PQA, PRHOA,       &
                     PLW_RAD, PRR, PSR, PZREF, PUREF, PVMOD, PH_TRAFFIC, PLE_TRAFFIC,   &
                     PTSTEP, PDF_RF, PDN_RF, PDF_RD, PDN_RD, PQSAT_RF, PQSAT_RD,        &
                     PDELT_RF, PDELT_RD, PTS_GARDEN, PLEW_RF, PUW_GR, PLEW_RD, PLE_WL_A,&
                     PLE_WL_B, PRNSN_RF, PHSN_RF, PLESN_RF, PGSN_RF, PMELT_RF, PRN_GR, &
                     PH_GR, PLE_GR, PGFLUX_GR, PDRAIN_GR, PRUNOFF_GR, PRNSN_RD,    &
                     PHSN_RD, PLESN_RD, PGSN_RD, PMELT_RD, PUW_RD, PUW_RF, PDUWDU_RD,   &
                     PDUWDU_RF, PUSTAR_TWN, PCD, PCDN, PCH_TWN, PRI_TWN, PRESA_TWN,     &
                     PAC_RF, PAC_RD, PAC_WL, PAC_TOP, PAC_GARDEN, PAC_RF_WAT,           &
                     PAC_RD_WAT, PLW_WA_TO_WB, PLW_WA_TO_R, PLW_WB_TO_R, PLW_WA_TO_NR,  &
                     PLW_WB_TO_NR, PLW_R_TO_WA, PLW_R_TO_WB, PLW_G_TO_WA, PLW_G_TO_WB,  &
                     PLW_S_TO_WA, PLW_S_TO_WB, PLW_S_TO_R, PLW_S_TO_NR, PLW_NR_TO_WA,   &
                     PLW_NR_TO_WB, PLW_NR_TO_WIN, PLW_WA_TO_WIN, PLW_WB_TO_WIN,         &
                     PLW_G_TO_WIN, PLW_R_TO_WIN, PLW_S_TO_WIN, PLW_WIN_TO_WA,           &
                     PLW_WIN_TO_WB, PLW_WIN_TO_R, PLW_WIN_TO_NR, KDAY, PEMIT_LW_FAC,    &
                     PEMIT_LW_RD, PT_RAD_IND, PHU_BLD, PTIME, PE_SHADING )
!   ##########################################################################
!
!!****  *TEB*  
!!
!!    PURPOSE
!!    -------
!
!     Computes the evoultion of prognostic variables and the fluxes
!     over artificial surfaces as towns, taking into account the canyon like
!     geometry of urbanized areas.
!         
!     
!!**  METHOD
!     ------
!
!     The prognostic variables are:
!       - the surface temperature for roofs, roads, and walls
!       - the water reservoir, whose maximum value is 10mm
!
!
!    1 : Warning about snow
!        ******************
!
!     Except for snow mantel evolution, all other computation with snow
!   variables must be performed with these variables at previous time-step,
!   and NOT new time-step. This insure coherence between snow fractions
!   (computed at the begining) and other snow characteristics (albedo, Ts).
!
!
!    2 : computation of input solar radiation on each surface
!        ****************************************************
!
!      Those are now done in subroutine urban_solar_abs.F90
!
!    3 : drag coefficient for momentum 
!        *****************************
!
!
!    4 : aerodynamical resistance for heat transfers
!        *******************************************
!
!
!    5 : equation for evolution of Ts_roof
!        *********************************
!
!
!       Rn = (dir_Rg + sca_Rg) (1-a) + emis * ( Rat - sigma Ts**4 (t+dt) )
!
!       H  = rho Cp CH V ( Ts (t+dt) - Tas )
!
!       LE = rho Lv CH V ( qs (t+dt) - qas )
!
!      where the as subscript denotes atmospheric values at ground level
!      (and not at first half level)
!
!
!    6 : equations for evolution of Ts_road and Ts_wall simultaneously
!        *************************************************************
!
!
!
!   Rn_w = abs_Rg_w 
!  - sigma * emis_w                                                   * Ts_w**4 (t+dt)
!  +         emis_w                       *      SVF_w                * Rat
!  + sigma * emis_w * emis_r              *      SVF_w                * Ts_r**4 (t+dt)
!  + sigma * emis_w * emis_w              * (1-2*SVF_w)               * Ts_w**4 (t+dt)
!  + sigma * emis_w * emis_w * (1-emis_r) *      SVF_w  * (1-  SVF_r) * Ts_w**4 (t+dt)
!  + sigma * emis_w * emis_w * (1-emis_w) * (1-2*SVF_w) * (1-2*SVF_w) * Ts_w**4 (t+dt)
!  + sigma * emis_w * emis_r * (1-emis_w) *      SVF_w  * (1-2*SVF_w) * Ts_r**4 (t+dt)
!
!   Rn_r = abs_Rg_r
!  - sigma * emis_r                                                   * Ts_r**4 (t+dt)
!  +         emis_r                       *    SVF_r                  * Rat
!  + sigma * emis_r * emis_w              * (1-SVF_r)                 * Ts_w**4 (t+dt)
!  + sigma * emis_r * emis_w * (1-emis_w) * (1-SVF_r)   * (1-2*SVF_w) * Ts_w**4 (t+dt)
!  + sigma * emis_r * emis_r * (1-emis_w) * (1-SVF_r)   *      SVF_w  * Ts_r**4 (t+dt)
!
!  H_w  = rho Cp CH V ( Ts_w (t+dt) - Ta_canyon )
!
!  LE_w = rho Lv CH V ( qs_w (t+dt) - qa_canyon )
!
!  H_r  = rho Cp CH V ( Ts_r (t+dt) - Ta_canyon )
!
!  LE_r = rho Lv CH V ( qs_r (t+dt) - qa_canyon )
!
! with again
!                AC_can * Swall/Sroad * Twall + AC_can * Troad + AC_top * Ta + H_traffic/Cp/rho/Sroad
!   Ta_canyon = -------------------------------------------------------------------------------------
!                AC_can * Swall/Sroad         + AC_can         + AC_top
!
!
!                 AC_can * delt_road * Hu_road * qsat(Troad) + AC_top * qa + LE_traffic/Lv/rho/Sroad
!   qa_canyon = ------------------------------------------------------------------------------------
!                 AC_can * delt_road                        + AC_top
!
!
!
!
!    7 : computation of fluxes for each surface type
!        *******************************************
!
!
!    8 : averaging of the fluxes
!        ***********************
!
!   This is done on the total exchange surface (roof + wall + road),
!  which is bigger than the horizontal surface (roof+road), leading
!  to bigger fluxes.
!
!   The fluxes due to industrial activity are directly added into the 
!  atmosphere
!
!
!    9 : road reservoir evolution
!        ************************
!
!   The roof reservoir runoff goes directly into the road reservoir.
!
!   Runoff occurs for road reservoir (too much water), as well as drainage
!   (evacuation system, typical time scale: 1 day)
!
!
!------------------------
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    MODD_CST
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      V. Masson           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    23/01/98 
!!     21 / 10 / 2003   P. Tulet    output aerodynamical resistance
!!     01 / 07 / 2005   P.Le Moigne Exner functions as arguments to urban_fluxes
!!     17 / 10 / 2005   (G. Pigeon) computation of anthropogenic heat from domestic heating
!!          01 / 2012   V. Masson   Separates the 2 walls 
!!     25 / 09 / 2012   B. Decharme new wind implicitation
!!          07 / 2013   V. Masson   Adds road watering
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!     ------------
!
USE modd_teb_option_n, ONLY : teb_options_t
USE modd_teb_n, ONLY : teb_t
USE modd_bem_option_n, ONLY : bem_options_t
USE modd_bem_n, ONLY : bem_t
USE modd_teb_irrig_n, ONLY : teb_irrig_t
USE modd_diag_misc_teb_n, ONLY : diag_misc_teb_t
!
USE modd_type_date_surf,ONLY: date_time
USE modd_csts,         ONLY : xtt, xstefan, xcpd, xlvtt
USE modd_surf_par,     ONLY : xundef
USE modd_snow_par,     ONLY : xemissn, xansmax_roof, &
                          xansmax_road,xwcrn_roof,xwcrn_road
!
USE mode_thermos
USE mode_surf_snow_frac
!
USE modi_snow_cover_1layer
USE modi_urban_drag
USE modi_urban_snow_evol
USE modi_roof_layer_e_budget
USE modi_road_layer_e_budget
USE modi_facade_e_budget
USE modi_urban_fluxes
USE modi_urban_hydro
USE modi_bld_e_budget
USE modi_wind_threshold
USE modi_bem
USE modi_teb_irrig
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    Declarations of arguments
!
TYPE(teb_options_t), INTENT(INOUT) :: TOP
TYPE(teb_t), INTENT(INOUT) :: T
TYPE(bem_options_t), INTENT(INOUT) :: BOP
TYPE(bem_t), INTENT(INOUT) :: B
TYPE(teb_irrig_t), INTENT(INOUT) :: TIR
TYPE(diag_misc_teb_t), INTENT(INOUT) :: DMT
!
 CHARACTER(LEN=*),     INTENT(IN)  :: HIMPLICIT_WIND   ! wind implicitation option
!                                                     ! 'OLD' = direct
!                                                     ! 'NEW' = Taylor serie, order 1
REAL, DIMENSION(:),   INTENT(IN)    :: PTSUN              ! solar time   (s from midnight)
REAL, DIMENSION(:), INTENT(INOUT) :: PT_CANYON     ! canyon air temperature
REAL, DIMENSION(:), INTENT(INOUT) :: PQ_CANYON     ! canyon air specific humidity
REAL, DIMENSION(:), INTENT(IN)    :: PU_CANYON     ! canyon hor. wind
REAL, DIMENSION(:), INTENT(IN)    :: PU_LOWCAN     ! wind near the road
REAL, DIMENSION(:), INTENT(IN)    :: PT_LOWCAN     ! temp. near the road
REAL, DIMENSION(:), INTENT(IN)    :: PQ_LOWCAN     ! hum. near the road
REAL, DIMENSION(:), INTENT(IN)    :: PZ_LOWCAN     ! height of atm. var. near the road
REAL, DIMENSION(:), INTENT(IN)    :: PPEW_A_COEF   ! implicit coefficients
REAL, DIMENSION(:), INTENT(IN)    :: PPEW_B_COEF   ! for wind coupling
REAL, DIMENSION(:), INTENT(IN)    :: PPEW_A_COEF_LOWCAN ! implicit coefficients for wind coupling
REAL, DIMENSION(:), INTENT(IN)    :: PPEW_B_COEF_LOWCAN ! between low canyon wind and road
REAL, DIMENSION(:), INTENT(IN)    :: PPS           ! pressure at the surface
REAL, DIMENSION(:), INTENT(IN)    :: PPA           ! pressure at the first atmospheric level
REAL, DIMENSION(:), INTENT(IN)    :: PEXNS         ! surface exner function
REAL, DIMENSION(:), INTENT(IN)    :: PTA           ! temperature at the lowest level
REAL, DIMENSION(:), INTENT(IN)    :: PQA           ! specific humidity
                                                   ! at the lowest level
REAL, DIMENSION(:), INTENT(IN)    :: PVMOD         ! module of the horizontal wind
REAL, DIMENSION(:), INTENT(IN)    :: PH_TRAFFIC    ! anthropogenic sensible
!                                                  ! heat fluxes due to traffic
REAL, DIMENSION(:), INTENT(IN)    :: PLE_TRAFFIC   ! anthropogenic latent
!                                                  ! heat fluxes due to traffic
REAL, DIMENSION(:), INTENT(IN)    :: PEXNA         ! exner function
                                                   ! at the lowest level
REAL, DIMENSION(:), INTENT(IN)    :: PRHOA         ! air density
                                                   ! at the lowest level
REAL, DIMENSION(:), INTENT(IN)    :: PLW_RAD       ! atmospheric infrared radiation
REAL, DIMENSION(:), INTENT(IN)    :: PRR           ! rain rate
REAL, DIMENSION(:), INTENT(IN)    :: PSR           ! snow rate
REAL, DIMENSION(:), INTENT(IN)    :: PZREF         ! reference height of the first
                                                   ! atmospheric level (temperature)
REAL, DIMENSION(:), INTENT(IN)    :: PUREF         ! reference height of the first
                                                   ! atmospheric level (wind)
REAL,               INTENT(IN)    :: PTSTEP        ! time step
!
REAL, DIMENSION(:), INTENT(INOUT) :: PDF_RF      ! snow-free    fraction on roofs
REAL, DIMENSION(:), INTENT(INOUT) :: PDN_RF      ! snow-covered fraction on roofs
REAL, DIMENSION(:), INTENT(INOUT) :: PDF_RD      ! snow-free    fraction on roads
REAL, DIMENSION(:), INTENT(INOUT) :: PDN_RD      ! snow-covered fraction on roads
REAL, DIMENSION(:), INTENT(OUT)   :: PQSAT_RF    ! hum at saturation over roof
REAL, DIMENSION(:), INTENT(OUT)   :: PQSAT_RD    ! hum at saturation over road
REAL, DIMENSION(:), INTENT(OUT)   :: PDELT_RF    ! water fraction on roof
REAL, DIMENSION(:), INTENT(OUT)   :: PDELT_RD    ! water fraction on road
!
REAL, DIMENSION(:), INTENT(IN)    :: PTS_GARDEN    ! GARDEN area surf temp.
! greenroof
REAL, DIMENSION(:), INTENT(OUT)   :: PLEW_RF    ! latent heat flux over roof (snow)
REAL, DIMENSION(:), INTENT(OUT)   :: PLEW_RD    ! latent heat flux over road (snow)
REAL, DIMENSION(:), INTENT(OUT)   :: PLE_WL_A   ! latent heat flux over wall
REAL, DIMENSION(:), INTENT(OUT)   :: PLE_WL_B   ! latent heat flux over wall
!
REAL, DIMENSION(:), INTENT(IN)    :: PUW_GR     ! Momentum flux for greenroofs
!
REAL, DIMENSION(:), INTENT(OUT)   :: PRNSN_RF ! net radiation over snow
REAL, DIMENSION(:), INTENT(OUT)   :: PHSN_RF  ! sensible heat flux over snow
REAL, DIMENSION(:), INTENT(OUT)   :: PLESN_RF ! latent heat flux over snow
REAL, DIMENSION(:), INTENT(OUT)   :: PGSN_RF  ! flux under the snow
REAL, DIMENSION(:), INTENT(OUT)   :: PMELT_RF   ! snow melt
REAL, DIMENSION(:), INTENT(OUT)   :: PRNSN_RD ! net radiation over snow
REAL, DIMENSION(:), INTENT(OUT)   :: PHSN_RD  ! sensible heat flux over snow
REAL, DIMENSION(:), INTENT(OUT)   :: PLESN_RD ! latent heat flux over snow
REAL, DIMENSION(:), INTENT(OUT)   :: PGSN_RD  ! flux under the snow
REAL, DIMENSION(:), INTENT(OUT)   :: PMELT_RD   ! snow melt
!
REAL, DIMENSION(:), INTENT(OUT)   :: PUW_RD     ! Momentum flux for roads
REAL, DIMENSION(:), INTENT(OUT)   :: PUW_RF     ! Momentum flux for roofs
REAL, DIMENSION(:), INTENT(OUT)   :: PDUWDU_RD  !
REAL, DIMENSION(:), INTENT(OUT)   :: PDUWDU_RF  !
REAL, DIMENSION(:), INTENT(OUT)   :: PUSTAR_TWN ! friciton velocity over town
!
REAL, DIMENSION(:), INTENT(IN)    :: PRN_GR     ! net radiation over greenroof
REAL, DIMENSION(:), INTENT(IN)    :: PH_GR      ! sensible heat flux over greenroof
REAL, DIMENSION(:), INTENT(IN)    :: PLE_GR     ! latent heat flux over greenroof
REAL, DIMENSION(:), INTENT(IN)    :: PGFLUX_GR  ! flux through the greenroof
REAL, DIMENSION(:), INTENT(IN)    :: PRUNOFF_GR ! runoff over green roofs
REAL, DIMENSION(:), INTENT(IN)    :: PDRAIN_GR  ! outlet drainage at base of green roofs
!
REAL, DIMENSION(:), INTENT(OUT)   :: PCD          ! town averaged drag coefficient
REAL, DIMENSION(:), INTENT(OUT)   :: PCDN         ! town averaged neutral drag coefficient
REAL, DIMENSION(:), INTENT(OUT)   :: PCH_TWN     ! town averaged heat transfer
!                                                 ! coefficient
REAL, DIMENSION(:), INTENT(OUT)   :: PRI_TWN      ! town averaged Richardson number
REAL, DIMENSION(:), INTENT(OUT)   :: PRESA_TWN    ! town aerodynamical resistance
REAL, DIMENSION(:), INTENT(OUT)   :: PAC_RF      ! roof conductance
REAL, DIMENSION(:), INTENT(INOUT) :: PAC_RD      ! road conductance
REAL, DIMENSION(:), INTENT(OUT)   :: PAC_WL      ! wall conductance
REAL, DIMENSION(:), INTENT(OUT)   :: PAC_TOP       ! top conductance
REAL, DIMENSION(:), INTENT(IN)    :: PAC_GARDEN    ! garden conductance
REAL, DIMENSION(:), INTENT(OUT)   :: PAC_RF_WAT  ! roof water conductance
REAL, DIMENSION(:), INTENT(OUT)   :: PAC_RD_WAT  ! roof water conductance
!
REAL, DIMENSION(:), INTENT(IN)    :: PLW_WA_TO_WB      ! LW contrib. wall A (orB) -> wall B (or A)
REAL, DIMENSION(:), INTENT(IN)    :: PLW_WA_TO_R         ! LW contrib. wall       -> road
REAL, DIMENSION(:), INTENT(IN)    :: PLW_WB_TO_R         ! LW contrib. wall       -> road
REAL, DIMENSION(:), INTENT(IN)    :: PLW_WA_TO_NR        ! LW contrib. wall       -> road(snow)
REAL, DIMENSION(:), INTENT(IN)    :: PLW_WB_TO_NR        ! LW contrib. wall       -> road(snow)
REAL, DIMENSION(:), INTENT(IN)    :: PLW_R_TO_WA         ! LW contrib. road       -> wall
REAL, DIMENSION(:), INTENT(IN)    :: PLW_R_TO_WB         ! LW contrib. road       -> wall
REAL, DIMENSION(:), INTENT(IN)    :: PLW_G_TO_WA         ! LW contrib. GARDEN     -> wall
REAL, DIMENSION(:), INTENT(IN)    :: PLW_G_TO_WB         ! LW contrib. GARDEN     -> wall
REAL, DIMENSION(:), INTENT(IN)    :: PLW_NR_TO_WA        ! LW contrib. road(snow) -> wall
REAL, DIMENSION(:), INTENT(IN)    :: PLW_NR_TO_WB        ! LW contrib. road(snow) -> wall
REAL, DIMENSION(:), INTENT(IN)    :: PLW_S_TO_WA         ! LW contrib. sky        -> wall
REAL, DIMENSION(:), INTENT(IN)    :: PLW_S_TO_WB         ! LW contrib. sky        -> wall
REAL, DIMENSION(:), INTENT(IN)    :: PLW_S_TO_R          ! LW contrib. sky        -> road
REAL, DIMENSION(:), INTENT(IN)    :: PLW_S_TO_NR         ! LW contrib. sky        -> road(snow)
!
! new arguments after BEM
!
INTEGER,            INTENT(IN)     :: KDAY         ! Simulation day
REAL, DIMENSION(:), INTENT(IN)   :: PLW_WA_TO_WIN ! Radiative heat trasfer coeff wall-window 
                                                  ! [W K-1 m-2] 
REAL, DIMENSION(:), INTENT(IN)   :: PLW_WB_TO_WIN ! Radiative heat trasfer coeff wall-window 
                                                  ! [W K-1 m-2] 
REAL, DIMENSION(:), INTENT(IN)   :: PLW_G_TO_WIN  ! Radiative heat trasfer coeff garden-window 
                                                  ! [W K-1 m-2]
REAL, DIMENSION(:), INTENT(IN)   :: PLW_R_TO_WIN  ! Radiative heat trasfer coeff road-window 
                                                  ! [W K-1 m-2] 
REAL, DIMENSION(:), INTENT(IN)   :: PLW_S_TO_WIN ! Radiative heat trasfer coeff window-sky 
                                                 ! [W K-1 m-2]
REAL, DIMENSION(:), INTENT(IN)   :: PLW_WIN_TO_WA! Radiative heat trasfer coeff window-wall
                                                 ! [W K-1 m-2] 
REAL, DIMENSION(:), INTENT(IN)   :: PLW_WIN_TO_WB! Radiative heat trasfer coeff window-wall
                                                 ! [W K-1 m-2] 
REAL, DIMENSION(:), INTENT(IN)   :: PLW_WIN_TO_R ! Radiative heat trasfer coeff window-road 
                                                 ! [W K-1 m-2]
REAL, DIMENSION(:), INTENT(IN)   :: PLW_NR_TO_WIN! Radiative heat trasfer coeff road(snow)-win 
                                                 ! [W K-1 m-2]
REAL, DIMENSION(:), INTENT(IN)   :: PLW_WIN_TO_NR! Radiative heat trasfer coeff win-road(snow) 
                                                 ! [W K-1 m-2]
 !new argument for PET calculation
REAL, DIMENSION(:), INTENT(OUT) :: PEMIT_LW_RD ! LW fluxes emitted by road (W/m2 surf road)
REAL, DIMENSION(:), INTENT(OUT) :: PEMIT_LW_FAC  ! LW fluxes emitted by wall (W/m2 surf wall)
REAL, DIMENSION(:), INTENT(OUT) :: PT_RAD_IND    ! Indoor mean radiant temperature [K]
REAL, DIMENSION(:), INTENT(OUT) :: PHU_BLD       ! Indoor relative humidity 0 < (-) < 1
REAL,                INTENT(IN)  :: PTIME        ! current time since midnight (UTC, s)
REAL, DIMENSION(:), INTENT(IN)  :: PE_SHADING    !energy not ref., nor absorbed, nor
                                                 !trans. by glazing [Wm-2(win)]
!
!*      0.2    Declarations of local variables
!
REAL, DIMENSION(SIZE(PTA)) :: ZVMOD          ! wind
REAL, DIMENSION(SIZE(PTA)) :: ZWS_RF_MAX   ! maximum deepness of roof
REAL, DIMENSION(SIZE(PTA)) :: ZWS_RD_MAX   ! and road water reservoirs
!
REAL, DIMENSION(SIZE(PTA)) :: ZAC_BLD        ! surface conductance inside the building itself in DEF building model
REAL, DIMENSION(SIZE(PTA)) :: ZTA            ! air temperature extrapolated at roof level
REAL, DIMENSION(SIZE(PTA)) :: ZQA            ! air humidity extrapolated at roof level
!
REAL, DIMENSION(SIZE(PTA)) :: ZDQS_RD      ! heat storage inside road
REAL, DIMENSION(SIZE(PTA)) :: ZDQS_RF      ! heat storage inside roof
REAL, DIMENSION(SIZE(PTA)) :: ZDQS_WL_A    ! heat storage inside wall
REAL, DIMENSION(SIZE(PTA)) :: ZDQS_WL_B    ! heat storage inside wall
REAL, DIMENSION(SIZE(PTA)) :: ZFLX_BLD_RF  !heat flux from inside through roof
REAL, DIMENSION(SIZE(PTA)) :: ZFLX_BLD_WL_A!heat flux from inside through wall
REAL, DIMENSION(SIZE(PTA)) :: ZFLX_BLD_WL_B!heat flux from inside through wall
REAL, DIMENSION(SIZE(PTA)) :: ZFLX_BLD_FL !heat flux from inside through floor
REAL, DIMENSION(SIZE(PTA)) :: ZFLX_BLD_MA  !heat flux from inside through mass
!
REAL, DIMENSION(SIZE(PTA)) :: ZDQS_SN_RF ! heat storage inside roof snowpack
REAL, DIMENSION(SIZE(PTA)) :: ZDQS_SN_RD ! heat storage inside road snowpack
REAL, DIMENSION(SIZE(PTA)) :: ZMELT_BLT      ! Snow melt for built & impervious part
!
! coefficients for LW computations over snow (from previous time-step)
!
REAL, DIMENSION(SIZE(PTA)) :: ZTSSN_RD ! road snow temperature
!                                          ! at previous time-step
! new local variables after BEM
!
REAL, DIMENSION(SIZE(PTA)) :: ZIMB_RF      ! residual energy imbalance
                                             ! of the roof for
                                             ! verification
REAL, DIMENSION(SIZE(PTA)) :: ZIMB_RD      ! road residual energy imbalance 
                                             ! for verification [W m-2]
REAL, DIMENSION(SIZE(PTA)) :: ZIMB_WL      ! wall residual energy imbalance 
                                             ! for verification [W m-2]
REAL, DIMENSION(SIZE(PTA)) :: ZTS_RD       ! road surface temperature 
!                                            ! at previous time-step
REAL, DIMENSION(SIZE(PTA)) :: ZTS_WL_A     ! wall A surface temperature 
!                                            ! at previous time-step
REAL, DIMENSION(SIZE(PTA)) :: ZTS_WL_B     ! wall B surface temperature 
!                                            ! at previous time-step
REAL, DIMENSION(SIZE(PTA)) :: ZTS_WL       ! averaged wall surface temperature 
!                                            ! at previous time-step
REAL, DIMENSION(SIZE(PTA)) :: ZTS_RF       ! roof surface temperature 
!                                            ! at previous time-step
REAL, DIMENSION(SIZE(PTA),SIZE(T%XT_WALL_A,2)) :: ZT_WL ! averaged wall surface temperature 
!
INTEGER :: IWL, IRF                      ! number of wall, roof layer
REAL, DIMENSION(SIZE(PTA)) :: ZRADHT_IN     ! Indoor radiant heat transfer coefficient
                                                    ! [W K-1 m-2]
REAL, DIMENSION(SIZE(PTA)) :: ZTS_FL       ! floor surface temperature [K]
REAL, DIMENSION(SIZE(PTA)) :: ZRAD_RF_WL  ! rad. flux from roof to averaged wall [W m-2(roof)]
REAL, DIMENSION(SIZE(PTA)) :: ZRAD_RF_WIN   ! rad. flux from roof to window [W m-2(roof)]
REAL, DIMENSION(SIZE(PTA)) :: ZRAD_RF_FL ! rad. flux from roof to floor [W m-2(roof)]
REAL, DIMENSION(SIZE(PTA)) :: ZRAD_RF_MA  ! rad. flux from roof to mass [W m-2(roof)]
REAL, DIMENSION(SIZE(PTA)) :: ZCONV_RF_BLD  ! rad. flux from roof to bld [W m-2(roof)]
REAL, DIMENSION(SIZE(PTA)) :: ZRAD_WL_FL ! rad. flux from averaged wall to floor [W m-2(wall)]
REAL, DIMENSION(SIZE(PTA)) :: ZRAD_WL_MA  ! rad. flux from averaged wall to mass [W m-2(wall)]
REAL, DIMENSION(SIZE(PTA)) :: ZRAD_WIN_FL  ! rad. flux from averaged wall to floor [W m-2(win)]
REAL, DIMENSION(SIZE(PTA)) :: ZRAD_WIN_MA   ! rad. flux from averaged wall to mass [W m-2(win)]
REAL, DIMENSION(SIZE(PTA)) :: ZCONV_WL_BLD  ! rad. flux from roof to bld [W m-2(wall)]
REAL, DIMENSION(SIZE(PTA)) :: ZCONV_WIN_BLD   ! rad. flux from roof to bld [W m-2(win)]
REAL, DIMENSION(SIZE(PTA)) :: ZAC_WIN         ! window aerodynamic conductance

REAL, DIMENSION(SIZE(PTA)) :: ZLOAD_IN_RF   ! indoor load on roof W/m2[roof]
REAL, DIMENSION(SIZE(PTA)) :: ZLOAD_IN_FL   ! indoor load on floor W/m2[floor]
REAL, DIMENSION(SIZE(PTA)) :: ZLOAD_IN_WL   ! indoor load on wall W/m2[wall]
REAL, DIMENSION(SIZE(PTA)) :: ZLOAD_IN_WIN   ! indoor load on win W/m2[win]
REAL, DIMENSION(SIZE(PTA)) :: ZLOAD_IN_MA   ! indoor load on mass W/m2[mass]
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE                                             
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('TEB',0,zhook_handle)
!
!*      1.     Initializations
!              ---------------
!
!*      1.1    Water reservoirs
!              ----------------
!
zws_rf_max =  1. ! (1mm) maximum deepness of roof water reservoir
zws_rd_max =  1. ! (1mm) maximum deepness of road water reservoir
!
!*      1.2    radiative snow variables at previous time-step
!              ----------------------------------------------
!
ztssn_rd(:) = t%TSNOW_ROAD%TS(:)
!
!
!*      1.3    indoor aerodynamique conductance for DEF case
!              ----------------------------------------------
!
zac_bld(:) = xundef
IF (top%CBEM=='DEF') zac_bld=1. / 0.123 / (xcpd * prhoa(:)) !* (normalized by rho Cp for convenience)
!-------------------------------------------------------------------------------
!
!*      1.3    number of roof/wall layer
!              -------------------------
!
iwl = SIZE(t%XT_WALL_A,2)
irf = SIZE(t%XT_ROOF,2)
!
zts_wl_a(:)=t%XT_WALL_A   (:,1)
zts_wl_b(:)=t%XT_WALL_B   (:,1)
zts_wl(:)=0.5 * (zts_wl_a(:)+zts_wl_b(:))
zts_rd(:)=t%XT_ROAD     (:,1)
zts_rf(:)=t%XT_ROOF     (:,1)
!
!
!*      1.4    load on indoor walls
!              -------------------------
!
IF (top%CBEM=='BEM') THEN
  !
  zload_in_rf = b%XF_FLOOR_WIN * dmt%XTR_SW_WIN + dmt%XQIN * b%XN_FLOOR * (1-b%XQIN_FLAT) * b%XQIN_FRAD  &
           / (2 + t%XWALL_O_BLD + b%XGLAZ_O_BLD + b%XMASS_O_BLD ) ! W/m2 [ROOF]
  zload_in_fl = b%XF_FLOOR_WIN * dmt%XTR_SW_WIN + dmt%XQIN * b%XN_FLOOR * (1-b%XQIN_FLAT) * b%XQIN_FRAD  &
           / (2 + t%XWALL_O_BLD + b%XGLAZ_O_BLD + b%XMASS_O_BLD )
  zload_in_ma = b%XF_MASS_WIN * dmt%XTR_SW_WIN + dmt%XQIN * b%XN_FLOOR * (1-b%XQIN_FLAT) * b%XQIN_FRAD  &
           / (2 + t%XWALL_O_BLD + b%XGLAZ_O_BLD + b%XMASS_O_BLD )
  zload_in_wl = b%XF_WALL_WIN * dmt%XTR_SW_WIN + dmt%XQIN * b%XN_FLOOR * (1-b%XQIN_FLAT) * b%XQIN_FRAD  &
           / (2 + t%XWALL_O_BLD + b%XGLAZ_O_BLD + b%XMASS_O_BLD )
  zload_in_win = b%XF_WIN_WIN * dmt%XTR_SW_WIN + dmt%XQIN * b%XN_FLOOR * (1-b%XQIN_FLAT) * b%XQIN_FRAD  &
           / (2 + t%XWALL_O_BLD + b%XGLAZ_O_BLD + b%XMASS_O_BLD )
ELSE
  zload_in_rf = 0.
  zload_in_fl = 0.
  zload_in_ma = 0.
  zload_in_wl = 0.
  zload_in_win = 0.
ENDIF
!-------------------------------------------------------------------------------
!
!*      2.     Snow-covered surfaces relative effects
!              --------------------------------------
!
!*      2.1    Effects on water reservoirs
!              ---------------------------
!
zws_rf_max(:) = zws_rf_max(:) * pdf_rf(:)
zws_rd_max(:) = zws_rd_max(:) * pdf_rd(:)
!
!-------------------------------------------------------------------------------
!
!*      3.     Surface drag
!              ------------
!
 CALL urban_drag(top, t, b, himplicit_wind, ptstep, pt_canyon, pq_canyon, &
                 pu_canyon, pt_lowcan, pq_lowcan, pu_lowcan, pz_lowcan, &
                 zts_rf, zts_rd, zts_wl, pts_garden, pdn_rf, pdn_rd,    &
                 pexns, pexna, pta, pqa, pps, prhoa, pzref, puref,      &
                 pvmod, zws_rf_max, zws_rd_max, ppew_a_coef,            &
                 ppew_b_coef, ppew_a_coef_lowcan, ppew_b_coef_lowcan,   &
                 pqsat_rf, pqsat_rd, pdelt_rf, pdelt_rd, pcd, pcdn,     &
                 pac_rf, pac_rf_wat, pac_wl, pac_rd, pac_rd_wat,        &
                 pac_top, pac_garden, pri_twn, puw_rd, puw_rf,          &
                 pduwdu_rd, pduwdu_rf, pustar_twn, zac_win    )
!
!* area-averaged heat transfer coefficient
!
zvmod(:) = wind_threshold(pvmod(:),puref(:))
!
pch_twn(:) = (t%XBLD(:) * pac_rf(:) + (1.-t%XBLD(:)) * pac_top(:)) / zvmod(:)
!
!* aggregation of momentum fluxes for roofs (=> derivate of flux also recalculated)
!
puw_rf(:) = (1-t%XGREENROOF(:)) * puw_rf(:) + t%XGREENROOF(:) * puw_gr(:)
WHERE (pvmod(:)/=0.) pduwdu_rf(:) = 2. * puw_rf(:) / pvmod(:)
!
!-------------------------------------------------------------------------------
!
!*      4.     Extrapolation of atmospheric T and q at roof level (for fluxes computation)
!              --------------------------------------------------
!
zta(:) = pta(:) * pexns(:) / pexna(:)
zqa(:) = pqa(:) * qsat(pta(:),pps(:)) / qsat(zta(:),ppa(:))
!
!-------------------------------------------------------------------------------
!
!*      5.     Snow mantel model
!              -----------------
!
 CALL urban_snow_evol(t, b, pt_lowcan, pq_lowcan, pu_lowcan, zts_rf, zts_rd, zts_wl_a,    &
                      zts_wl_b, pps, zta, zqa, prhoa, plw_rad, psr, pzref, puref, pvmod,  &
                      ptstep,  pz_lowcan, pdn_rf, dmt%XABS_SW_SNOW_ROOF,                 &
                      dmt%XABS_LW_SNOW_ROOF, pdn_rd, dmt%XABS_SW_SNOW_ROAD,             &
                      dmt%XABS_LW_SNOW_ROAD, prnsn_rf, phsn_rf, plesn_rf, pgsn_rf,       &
                      pmelt_rf, prnsn_rd, phsn_rd, plesn_rd, pgsn_rd, pmelt_rd,           &
                      plw_wa_to_nr, plw_wb_to_nr, plw_s_to_nr, plw_win_to_nr, zdqs_sn_rf, &
                      zdqs_sn_rd      )
!
!-------------------------------------------------------------------------------
!
!*      6.    LW properties
!              -------------
!
pdf_rd(:) = 1. - pdn_rd(:)
!
!-------------------------------------------------------------------------------
!
!*      7.    Indoor radiative temperature
!              ---------------------------
!
! uses the averaged temperature of both walls for the building energy balance
zt_wl(:,:)=0.5 * (t%XT_WALL_A(:,:)+t%XT_WALL_B(:,:))
!
SELECT CASE(top%CBEM)
   CASE("DEF")
      zts_fl(:) = 19. + xtt
      pt_rad_ind(:) = ( t%XWALL_O_HOR(:) / t%XBLD(:) * zt_wl(:,iwl) + &
                  t%XT_ROOF(:,irf) + zts_fl(:) ) / (t%XWALL_O_HOR(:) / t%XBLD(:) + 1. + 1.) 
      zradht_in(:) = xundef
   CASE("BEM")
      zts_fl(:) = b%XT_FLOOR(:,1)
      pt_rad_ind(:)  = (b%XT_MASS(:,1)*b%XMASS_O_BLD(:) + zt_wl(:,iwl)*t%XWALL_O_BLD(:)     &
                   + zts_fl(:) + t%XT_ROOF(:,irf) + b%XT_WIN2(:) * b%XGLAZ_O_BLD(:)) &
                   /(b%XMASS_O_BLD(:) + t%XWALL_O_BLD(:) + 1. + 1. + b%XGLAZ_O_BLD(:))
      !             Assuming indoor surface emissivities of 0.9
      zradht_in(:)   = 0.9 * 0.9 * 4 * xstefan * pt_rad_ind(:)**3          
END SELECT
!
!
!*      7.    Roof Ts computation
!              -------------------
!
!* ts_roof and qsat_roof are updated
!

 CALL roof_layer_e_budget(top, t, b, pqsat_rf, zac_bld, ptstep, pdn_rf, prhoa,    &
                          pac_rf, pac_rf_wat, plw_rad, pps, pdelt_rf, zta, zqa,   &
                          pexna, pexns, dmt%XABS_SW_ROOF, pgsn_rf,  zflx_bld_rf, &
                          zdqs_rf, dmt%XABS_LW_ROOF, dmt%XH_ROOF, plew_rf, zimb_rf, &
                          dmt%XG_GREENROOF_ROOF, zradht_in, zts_fl, zt_wl(:,iwl),&
                          zrad_rf_wl, zrad_rf_win, zrad_rf_fl, zrad_rf_ma, zconv_rf_bld, &
                          prr, & !modif to add heating/cooling of rain
                          zload_in_rf )
!
!-------------------------------------------------------------------------------
!
!*      8.    Road Ts computations
!              -----------------------------
!
!* Road watering

 CALL teb_irrig(tir%LPAR_RD_IRRIG, ptstep, top%TTIME%TDATE%MONTH, ptsun,   &
               tir%XRD_START_MONTH, tir%XRD_END_MONTH, tir%XRD_START_HOUR, &
               tir%XRD_END_HOUR, tir%XRD_24H_IRRIG, dmt%XIRRIG_ROAD      )

!* ts_road, ts_wall, qsat_road, t_canyon and q_canyon are updated
!
 CALL road_layer_e_budget(t, b, ptstep, pdn_rd, prhoa, pac_rd, pac_rd_wat, &
                          plw_rad, pps, pqsat_rd, pdelt_rd, pexns,         &
                          dmt%XABS_SW_ROAD, pgsn_rd, pq_lowcan, pt_lowcan,&
                          zts_wl_a, zts_wl_b, ztssn_rd,  pts_garden,       &
                          plw_wa_to_r, plw_wb_to_r, plw_s_to_r,            &
                          plw_win_to_r, pemit_lw_rd, zdqs_rd, dmt%XABS_LW_ROAD,  &
                          dmt%XH_ROAD, plew_rd, zimb_rd, prr+dmt%XIRRIG_ROAD    )
!
!-------------------------------------------------------------------------------
!
!*      8.     Wall Ts computations
!              -----------------------------
!
 CALL facade_e_budget(top, t, b, dmt, ptstep, pdn_rd, prhoa, pac_wl, zac_bld,   &
                      plw_rad, pps, pexns, pt_canyon, zts_rd, ztssn_rd, pts_garden, &
                      zts_fl, plw_wa_to_wb, plw_r_to_wa, plw_r_to_wb,      &
                      plw_g_to_wa, plw_g_to_wb, plw_s_to_wa, plw_s_to_wb,  &
                      plw_nr_to_wa, plw_nr_to_wb, plw_win_to_wa,           &
                      plw_win_to_wb, plw_s_to_win, plw_wa_to_win,          &
                      plw_wb_to_win, plw_r_to_win, plw_g_to_win,           &
                      plw_nr_to_win, zflx_bld_wl_a, zdqs_wl_a,             &
                      zflx_bld_wl_b, zdqs_wl_b, pemit_lw_fac, zimb_wl,     &
                      zradht_in, zrad_rf_wl, zrad_rf_win, zrad_wl_fl,      &
                      zrad_wl_ma, zrad_win_fl, zrad_win_ma, zconv_wl_bld,  &
                      zconv_win_bld, zac_win, zload_in_wl, zload_in_win   )
!
!-------------------------------------------------------------------------------
!
!*      9.     Evolution of interior building air temperature
!              ----------------------------------------------
!
! uses the averaged temperature of both walls for the building energy balance
zt_wl(:,:)=0.5 * (t%XT_WALL_A(:,:)+t%XT_WALL_B(:,:))
!
SELECT CASE(top%CBEM)
CASE("DEF")
!
   CALL bld_e_budget(.true., ptstep, t%XBLD, t%XWALL_O_HOR,        &
                     prhoa, t%XT_ROOF, zt_wl, b%XTI_BLD, zts_fl(:) )

   !variables that needs to be computed apart
   b%XQI_BLD = 0.5 * qsat(b%XTI_BLD, pps)
   !variables that need to be set 0 for calculation
   zflx_bld_fl(:) = 0.
   zflx_bld_ma(:) = 0.
   !other variables
   phu_bld(:)     = xundef

CASE("BEM")
  CALL bem(bop, t, b, dmt, ptstep, ptsun, kday, pps, prhoa, pt_canyon, &
           pq_canyon, pu_canyon, phu_bld, pt_rad_ind, zflx_bld_fl,&
           zflx_bld_ma, zradht_in, zrad_rf_ma, zrad_rf_fl,        &
           zrad_wl_ma, zrad_wl_fl, zrad_win_ma, zrad_win_fl,      &
           zconv_rf_bld, zconv_wl_bld, zconv_win_bld, zload_in_fl,&
           zload_in_ma                                 )

   dmt%XH_WASTE  = dmt%XH_WASTE  * t%XBLD
   dmt%XLE_WASTE = dmt%XLE_WASTE * t%XBLD
END SELECT
!
!-------------------------------------------------------------------------------
!
!*      10.    Fluxes over built surfaces
!              --------------------------
!
 CALL urban_fluxes   (top, t, b, dmt, himplicit_wind, pt_canyon, ppew_a_coef, ppew_b_coef,      &
                      pexns, prhoa, pvmod, ph_traffic, ple_traffic,pac_wl, pcd, pdf_rf,         &
                      pdn_rf, pdf_rd, pdn_rd, prnsn_rf, phsn_rf, plesn_rf, pgsn_rf,             &
                      prnsn_rd, phsn_rd, plesn_rd, pgsn_rd, pmelt_rf, zdqs_rf, pmelt_rd,        &
                      zdqs_rd, zdqs_wl_a, zdqs_wl_b, zflx_bld_rf, zflx_bld_wl_a,                &
                      zflx_bld_wl_b, zflx_bld_fl, zflx_bld_ma, pe_shading, plew_rf,             &
                      prn_gr, ph_gr, ple_gr, pgflux_gr,                                         &
                      plew_rd, ple_wl_a, ple_wl_b, zmelt_blt, pustar_twn                        )
!
!
! Water transfer from snow reservoir to water reservoir in case of snow melt
!
WHERE (pmelt_rf(:) .GT. 0.)
  t%XWS_ROOF(:) = min(zws_rf_max,t%XWS_ROOF(:) + pmelt_rf(:)*ptstep)
ENDWHERE
!
WHERE (pmelt_rd(:) .GT. 0.)
  t%XWS_ROAD(:) = min(zws_rd_max,t%XWS_ROAD(:) + pmelt_rd(:)*ptstep)
ENDWHERE
!
!-------------------------------------------------------------------------------
!
!*      11.    Roof ans road reservoirs evolution
!              ----------------------------------
!
 CALL urban_hydro(zws_rf_max, zws_rd_max, t%XWS_ROOF, t%XWS_ROAD, prr,          &
                  dmt%XIRRIG_ROAD, ptstep, t%XBLD, dmt%XLE_ROOF, dmt%XLE_ROAD,  &
                  dmt%XRUNOFF_STRLROOF, dmt%XRUNOFF_ROAD   )
!
IF (top%LGREENROOF) THEN
  dmt%XRUNOFF_ROOF(:) =  (1.-t%XGREENROOF(:)) * dmt%XRUNOFF_STRLROOF(:) &
                        + t%XGREENROOF(:) * (prunoff_gr(:) + pdrain_gr(:))
ELSE
  dmt%XRUNOFF_ROOF(:) =  dmt%XRUNOFF_STRLROOF(:)
ENDIF                                                      
!
!-------------------------------------------------------------------------------
!
!*      19.    Compute aerodynamical resistance 
!              --------------------------------
!
presa_twn(:) = 1. / ( t%XBLD(:) * pac_rf(:)  + ( 1. - t%XBLD(:)) * pac_top(:))
!
IF (lhook) CALL dr_hook('TEB',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE teb