MODULE module_sf_bep_bem

#ifdef mpas
!USE mpas_dmpar, only : mpas_dmpar_global_abort
#define FATAL_ERROR(M) call mpas_dmpar_global_abort( M )
#else 
#define FATAL_ERROR(M) write(0,*) M ; stop
#endif

!USE module_model_constants
 USE module_sf_urban
 USE module_sf_bem

! SGClarke 09/11/2008
! Access urban_param.tbl values through calling urban_param_init in module_physics_init
! for CASE (BEPSCHEME) select sf_urban_physics
!
! -----------------------------------------------------------------------
!  Dimension for the array used in the BEP module
! -----------------------------------------------------------------------

      integer nurbm           ! Maximum number of urban classes    
      parameter (nurbm=3)

      integer ndm             ! Maximum number of street directions 
      parameter (ndm=2)

      integer nz_um           ! Maximum number of vertical levels in the urban grid
      parameter(nz_um=13)

      integer ng_u            ! Number of grid levels in the ground
      parameter (ng_u=10)
      integer nwr_u            ! Number of grid levels in the walls or roofs
      parameter (nwr_u=10)

      integer nf_u             !Number of grid levels in the floors (BEM)
      parameter (nf_u=10)

      integer ngb_u            !Number of grid levels in the ground below building (BEM)
      parameter (ngb_u=10)

      real dz_u                ! Urban grid resolution
      parameter (dz_u=5.) 

      integer nbui_max         !maximum number of types of buildings in an urban class
      parameter (nbui_max=4)   !must be less or equal than nz_um 

!---------------------------------------------------------------------------------
!Parameters of the windows. The glasses of windows are considered without films  -
!Read the paper of J.Karlsson and A.Roos(2000):"modelling the angular behaviour  -
!of the total solar energy transmittance of windows".Solar Energy Vol.69,No.4,   -
!pp. 321-329, for more details.                                                  -
!---------------------------------------------------------------------------------
      integer p_num            !number of panes in the windows (1,2 or 3)
      parameter (p_num=2)
      integer q_num            !category number for the windows (q_num= 4, standard glasses)
      parameter(q_num=4)       !Possible values 1,2,...,10


! The change of ng_u, nwr_u should be done in agreement with the block data
!  in the routine "surf_temp" 
! -----------------------------------------------------------------------
!  Constant used in the BEP module
! -----------------------------------------------------------------------
           
      real vk                 ! von Karman constant
      real g_u                ! Gravity acceleration
      real pi                 !
      real r                  ! Perfect gas constant
      real cp_u               ! Specific heat at constant pressure
      real rcp_u              !
      real sigma              !
      real p0                 ! Reference pressure at the sea level
      real cdrag              ! Drag force constant
      real latent             ! Latent heat of vaporization [J/kg] (used in BEM)
      parameter(vk=0.40,g_u=9.81,pi=3.141592653,r=287.,cp_u=1004.)        
      parameter(rcp_u=r/cp_u,sigma=5.67e-08,p0=1.e+5,cdrag=0.4,latent=2.45e+06)

! -----------------------------------------------------------------------     




   CONTAINS
 
      subroutine BEP_BEM(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy,      &
                      th_phy,rho,p_phy,swdown,glw,                                 &
                      gmt,julday,xlong,xlat,                                       &
                      declin_urb,cosz_urb2d,omg_urb2d,                             &
                      num_urban_layers,                                            &
                      trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d,                     &
                      tlev_urb3d,qlev_urb3d,tw1lev_urb3d,tw2lev_urb3d,             &
                      tglev_urb3d,tflev_urb3d,sf_ac_urb3d,lf_ac_urb3d,             &
                      cm_ac_urb3d,sfvent_urb3d,lfvent_urb3d,                       &
                      sfwin1_urb3d,sfwin2_urb3d,                                   &
                      sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d,                   &
                      a_u,a_v,a_t,a_e,b_u,b_v,                                     &
                      b_t,b_e,b_q,dlg,dl_u,sf,vl,                                  &
                      rl_up,rs_abs,emiss,grdflx_urb,qv_phy,                        &
                      ids,ide, jds,jde, kds,kde,                                   &
                      ims,ime, jms,jme, kms,kme,                                   &
                      its,ite, jts,jte, kts,kte)                    

      implicit none

!------------------------------------------------------------------------
!     Input
!------------------------------------------------------------------------
   INTEGER ::                       ids,ide, jds,jde, kds,kde,  &
                                    ims,ime, jms,jme, kms,kme,  &
                                    its,ite, jts,jte, kts,kte,  &
                                    itimestep
 

   REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   DZ8W
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   P_PHY
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   RHO
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   TH_PHY
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   T_PHY
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   U_PHY
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   V_PHY
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   U
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme )::   V
   REAL, DIMENSION( ims:ime , jms:jme )        ::   GLW
   REAL, DIMENSION( ims:ime , jms:jme )        ::   swdown
   REAL, DIMENSION( ims:ime, jms:jme )         ::   UST
   INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN )::   UTYPE_URB2D
   REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN )::   FRC_URB2D
   REAL, INTENT(IN  )   ::                                   GMT 
   INTEGER, INTENT(IN  ) ::                               JULDAY
   REAL, DIMENSION( ims:ime, jms:jme ),                           &
         INTENT(IN   )  ::                           XLAT, XLONG
   REAL, INTENT(IN) :: DECLIN_URB
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
   INTEGER, INTENT(IN  ) :: num_urban_layers
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: trb_urb4d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1_urb4d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2_urb4d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tgb_urb4d
!New variables used for BEM
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: qv_phy
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tlev_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: qlev_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1lev_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2lev_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tglev_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tflev_urb3d
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lf_ac_urb3d
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sf_ac_urb3d
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: cm_ac_urb3d
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfvent_urb3d
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lfvent_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin1_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin2_urb3d
!End variables
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw1_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw2_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfr_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfg_urb3d

   real z(ims:ime,kms:kme,jms:jme)            ! Vertical coordinates
   REAL, INTENT(IN )::   DT      ! Time step

!      
!------------------------------------------------------------------------
!     Output
!------------------------------------------------------------------------ 
!
!    Implicit and explicit components of the source and sink terms at each levels,
!     the fluxes can be computed as follow: FX = A*X + B   example: t_fluxes = a_t * pt + b_t
      real a_u(ims:ime,kms:kme,jms:jme)         ! Implicit component for the momemtum in X-direction (center)
      real a_v(ims:ime,kms:kme,jms:jme)         ! Implicit component for the momemtum in Y-direction (center)
      real a_t(ims:ime,kms:kme,jms:jme)         ! Implicit component for the temperature
      real a_e(ims:ime,kms:kme,jms:jme)         ! Implicit component for the TKE
      real b_u(ims:ime,kms:kme,jms:jme)         ! Explicit component for the momemtum in X-direction (center)
      real b_v(ims:ime,kms:kme,jms:jme)         ! Explicit component for the momemtum in Y-direction (center)
      real b_t(ims:ime,kms:kme,jms:jme)         ! Explicit component for the temperature
      real b_e(ims:ime,kms:kme,jms:jme)         ! Explicit component for the TKE
      real b_q(ims:ime,kms:kme,jms:jme)         ! Explicit component for the Humidity
      real dlg(ims:ime,kms:kme,jms:jme)         ! Height above ground (L_ground in formula (24) of the BLM paper). 
      real dl_u(ims:ime,kms:kme,jms:jme)        ! Length scale (lb in formula (22) ofthe BLM paper).
! urban surface and volumes        
      real sf(ims:ime,kms:kme,jms:jme)           ! surface of the urban grid cells
      real vl(ims:ime,kms:kme,jms:jme)             ! volume of the urban grid cells
! urban fluxes
      real rl_up(ims:ime,jms:jme) ! upward long wave radiation
      real rs_abs(ims:ime,jms:jme) ! absorbed short wave radiation
      real emiss(ims:ime,jms:jme)  ! emissivity averaged for urban surfaces
      real grdflx_urb(ims:ime,jms:jme)  ! ground heat flux for urban areas
!------------------------------------------------------------------------
!     Local
!------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    Building parameters      
      real alag_u(nurbm)                      ! Ground thermal diffusivity [m^2 s^-1]
      real alaw_u(nurbm)                      ! Wall thermal diffusivity [m^2 s^-1]
      real alar_u(nurbm)                      ! Roof thermal diffusivity [m^2 s^-1]
      real csg_u(nurbm)                       ! Specific heat of the ground material [J m^3 K^-1]
      real csw_u(nurbm)                       ! Specific heat of the wall material [J m^3 K^-1]
      real csr_u(nurbm)                       ! Specific heat of the roof material [J m^3 K^-1]
      real twini_u(nurbm)                     ! Initial temperature inside the building's wall [K]
      real trini_u(nurbm)                     ! Initial temperature inside the building's roof [K]
      real tgini_u(nurbm)                     ! Initial road temperature
!
! for twini_u, and trini_u the initial value at the deepest level is kept constant during the simulation
!
!    Radiation paramters
      real albg_u(nurbm)                      ! Albedo of the ground
      real albw_u(nurbm)                      ! Albedo of the wall
      real albr_u(nurbm)                      ! Albedo of the roof
      real albwin_u(nurbm)                    ! Albedo of the windows
      real emwind_u(nurbm)                    ! Emissivity of windows
      real emg_u(nurbm)                       ! Emissivity of ground
      real emw_u(nurbm)                       ! Emissivity of wall
      real emr_u(nurbm)                       ! Emissivity of roof

!   fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
!   and the short wave radation. 
      real fww(nz_um,nz_um,ndm,nurbm)         !  from wall to wall
      real fwg(nz_um,ndm,nurbm)               !  from wall to ground
      real fgw(nz_um,ndm,nurbm)               !  from ground to wall
      real fsw(nz_um,ndm,nurbm)               !  from sky to wall
      real fws(nz_um,ndm,nurbm)               !  from sky to wall
      real fsg(ndm,nurbm)                     !  from sky to ground

!    Roughness parameters
      real z0g_u(nurbm)         ! The ground's roughness length
      real z0r_u(nurbm)         ! The roof's roughness length

!    Street parameters
      integer nd_u(nurbm)       ! Number of street direction for each urban class 
      real strd_u(ndm,nurbm)    ! Street length (fix to greater value to the horizontal length of the cells)
      real drst_u(ndm,nurbm)    ! Street direction
      real ws_u(ndm,nurbm)      ! Street width
      real bs_u(ndm,nurbm)      ! Building width
      real h_b(nz_um,nurbm)     ! Bulding's heights
      real d_b(nz_um,nurbm)     ! Probability that a building has an height h_b
      real ss_u(nz_um,nurbm)  ! Probability that a building has an height equal to z
      real pb_u(nz_um,nurbm)  ! Probability that a building has an height greater or equal to z


!    Grid parameters
      integer nz_u(nurbm)       ! Number of layer in the urban grid
      
      real z_u(nz_um)         ! Height of the urban grid levels

!   MT
      real cop_u(nurbm)
      real pwin_u(nurbm)
      real beta_u(nurbm)
      integer sw_cond_u(nurbm)
      real time_on_u(nurbm)
      real time_off_u(nurbm)
      real targtemp_u(nurbm)
      real gaptemp_u(nurbm)
      real targhum_u(nurbm)
      real gaphum_u(nurbm)
      real perflo_u(nurbm)
      real hsesf_u(nurbm)
      real hsequip(24)

! 1D array used for the input and output of the routine "urban"

      real z1D(kms:kme)               ! vertical coordinates
      real ua1D(kms:kme)                ! wind speed in the x directions
      real va1D(kms:kme)                ! wind speed in the y directions
      real pt1D(kms:kme)                ! potential temperature
      real da1D(kms:kme)                ! air density
      real pr1D(kms:kme)                ! air pressure
      real pt01D(kms:kme)               ! reference potential temperature
      real zr1D                    ! zenith angle
      real deltar1D                ! declination of the sun
      real ah1D                    ! hour angle (it should come from the radiation routine)
      real rs1D                    ! solar radiation
      real rld1D                   ! downward flux of the longwave radiation

      real tw1D(2*ndm,nz_um,nwr_u,nbui_max) ! temperature in each layer of the wall
      real tg1D(ndm,ng_u)                   ! temperature in each layer of the ground
      real tr1D(ndm,nz_um,nwr_u)   ! temperature in each layer of the roof
!
!New variable for BEM
!
      real tlev1D(nz_um,nbui_max)            ! temperature in each floor and in each different type of building
      real qlev1D(nz_um,nbui_max)            ! specific humidity in each floor and in each different type of building
      real twlev1D(2*ndm,nz_um,nbui_max)     ! temperature in each window in each floor in each different type of building
      real tglev1D(ndm,ngb_u,nbui_max)       ! temperature in each layer of the ground below of a type of building
      real tflev1D(ndm,nf_u,nz_um-1,nbui_max)! temperature in each layer of the floors in each building
      real lflev1D(nz_um,nz_um)           ! latent heat flux due to the air conditioning systems
      real sflev1D(nz_um,nz_um)           ! sensible heat flux due to the air conditioning systems
      real lfvlev1D(nz_um,nz_um)          ! latent heat flux due to ventilation
      real sfvlev1D(nz_um,nz_um)          ! sensible heat flux due to ventilation
      real sfwin1D(2*ndm,nz_um,nbui_max)     ! sensible heat flux from windows
      real consumlev1D(nz_um,nz_um)       ! consumption due to the air conditioning systems
      real qv1D(kms:kme)                  ! specific humidity
      real meso_urb                       ! constant to link meso and urban scales [mÂ¯2]
      real d_urb(nz_um)    
      real sf_ac
      integer ibui,nbui
      integer nlev(nz_um)
!
!End new variables
!

      real sfw1D(2*ndm,nz_um,nbui_max)      ! sensible heat flux from walls
      real sfg1D(ndm)              ! sensible heat flux from ground (road)
      real sfr1D(ndm,nz_um)      ! sensible heat flux from roofs
      real sf1D(kms:kme)              ! surface of the urban grid cells
      real vl1D(kms:kme)                ! volume of the urban grid cells
      real a_u1D(kms:kme)               ! Implicit component of the momentum sources or sinks in the X-direction
      real a_v1D(kms:kme)               ! Implicit component of the momentum sources or sinks in the Y-direction
      real a_t1D(kms:kme)               ! Implicit component of the heat sources or sinks
      real a_e1D(kms:kme)               ! Implicit component of the TKE sources or sinks
      real b_u1D(kms:kme)               ! Explicit component of the momentum sources or sinks in the X-direction
      real b_v1D(kms:kme)               ! Explicit component of the momentum sources or sinks in the Y-direction
      real b_t1D(kms:kme)               ! Explicit component of the heat sources or sinks
      real b_ac1D(kms:kme)
      real b_e1D(kms:kme)               ! Explicit component of the TKE sources or sinks
      real b_q1D(kms:kme)               ! Explicit component of the Humidity sources or sinks
      real dlg1D(kms:kme)               ! Height above ground (L_ground in formula (24) of the BLM paper). 
      real dl_u1D(kms:kme)              ! Length scale (lb in formula (22) ofthe BLM paper)
      real time_bep
! arrays used to collapse indexes
      integer ind_zwd(nbui_max,nz_um,nwr_u,ndm)
      integer ind_gd(ng_u,ndm)
      integer ind_zd(nbui_max,nz_um,ndm)
      integer ind_zdf(nz_um,ndm)
      integer ind_zrd(nz_um,nwr_u,ndm)
!
      integer ind_bd(nbui_max,nz_um)
      integer ind_wd(nbui_max,nz_um,ndm)
      integer ind_gbd(nbui_max,ngb_u,ndm)  
      integer ind_fbd(nbui_max,nf_u,nz_um-1,ndm)
!
      integer ix,iy,iz,iurb,id,iz_u,iw,ig,ir,ix1,iy1,k
      integer it, nint
      integer iii
      real tempo
      logical first
      character(len=80) :: text
      data first/.true./
      save first,time_bep 

      save alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,                      &
          albg_u,albw_u,albr_u,emg_u,emw_u,emr_u,fww,fwg,fgw,fsw,fws,fsg, &
          z0g_u,z0r_u, nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,  &
          nz_u,z_u,albwin_u,emwind_u ,                                  &
          cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, targtemp_u,  &
          gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip

!------------------------------------------------------------------------
!    Calculation of the momentum, heat and turbulent kinetic fluxes
!     produced by builgings
!
! Reference:
! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
! 261-304
!
! F. Salamanca and A. Martilli, 2009: 'A new Building Energy Model coupled 
! with an Urban Canopy Parameterization for urban climate simulations_part II. 
! Validation with one dimension off-line simulations'. Theor Appl Climatol
! DOI 10.1007/s00704-009-0143-8 
!------------------------------------------------------------------------
!prepare the arrays to collapse indexes
      if(num_urban_layers.lt.nbui_max*nz_um*ndm*max(nwr_u,ng_u))then
        write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*nz_um*ndm*max(nwr_u,ng_u)
        stop
      endif
!
!New conditions for BEM
!
      if(num_urban_layers.lt.nbui_max*nz_um)then !limit for indoor temperature and indoor humidity
        write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*nz_um
        stop
      endif

      if(num_urban_layers.lt.nbui_max*nz_um*ndm)then !limit for window temperature
        write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*nz_um*ndm
        stop
      endif

      if(num_urban_layers.lt.nbui_max*ndm*ngb_u)then !limit for ground temperature below a building
        write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*ndm*ngb_u
        stop
      endif

      if(num_urban_layers.lt.(nz_um-1)*nbui_max*ndm*nf_u)then !limit for floor temperature 
        write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*ndm*nf_u*(nz_um-1),num_urban_layers
        stop
      endif

      if (ndm.ne.2)then
         write(*,*) 'number of directions is not correct',ndm
         stop
      endif

!
!End of new conditions
!
!
!Initialize collapse indexes
!
      ind_zwd=0       
      ind_gd=0
      ind_zd=0
      ind_zdf=0
      ind_zrd=0
      ind_bd=0
      ind_wd=0
      ind_gbd=0
      ind_fbd=0
!
!End initialization indexes
!

      iii=0
      do ibui=1,nbui_max
      do iz_u=1,nz_um
      do iw=1,nwr_u
      do id=1,ndm
       iii=iii+1
       ind_zwd(ibui,iz_u,iw,id)=iii
      enddo
      enddo
      enddo
      enddo
      
      iii=0
      do ig=1,ng_u
      do id=1,ndm
       iii=iii+1
       ind_gd(ig,id)=iii
      enddo
      enddo
      
      iii=0
      do ibui=1,nbui_max
      do iz_u=1,nz_um
      do id=1,ndm
       iii=iii+1
       ind_zd(ibui,iz_u,id)=iii
      enddo
      enddo
      enddo
      
      iii=0
      do iz_u=1,nz_um
      do iw=1,nwr_u
      do id=1,ndm
       iii=iii+1
       ind_zrd(iz_u,iw,id)=iii
      enddo
      enddo
      enddo
      
!
!New indexes for BEM
!
      iii=0
      do iz_u=1,nz_um
      do id=1,ndm
         iii=iii+1
         ind_zdf(iz_u,id)=iii
      enddo ! id
      enddo ! iz_u
      
      iii=0
      do ibui=1,nbui_max  !Type of building
      do iz_u=1,nz_um     !vertical levels
         iii=iii+1
         ind_bd(ibui,iz_u)=iii
      enddo !iz_u
      enddo !ibui
      
      
      iii=0
      do ibui=1,nbui_max !type of building
      do iz_u=1,nz_um !vertical levels
      do id=1,ndm !direction
         iii=iii+1
         ind_wd(ibui,iz_u,id)=iii
      enddo !id
      enddo !iz_u
      enddo !ibui
      
      iii=0
      do ibui=1,nbui_max!type of building
      do iw=1,ngb_u !layers in the wall (ground below a building)
      do id=1,ndm !direction
         iii=iii+1
         ind_gbd(ibui,iw,id)=iii  
      enddo !id
      enddo !iw 
      enddo !ibui    
      
      iii=0
      do ibui=1,nbui_max !type of building
      do iw=1,nf_u !layers in the wall (floor)
      do iz_u=1,nz_um-1 !vertical levels
      do id=1,ndm  !direction
         iii=iii+1
         ind_fbd(ibui,iw,iz_u,id)=iii
      enddo !id
      enddo !iz_u
      enddo !iw
      enddo !ibui
      
!
!End of new indexes
!      

      do ix=its,ite
      do iy=jts,jte
       z(ix,kts,iy)=0.
       do iz=kts+1,kte+1
        z(ix,iz,iy)=z(ix,iz-1,iy)+dz8w(ix,iz-1,iy)
       enddo
      enddo
      enddo

      if (first) then                           ! True only on first call
         call init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
                twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,&
                emr_u,emwind_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,&
                cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, &
                targtemp_u, gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip)

!Initialisation of the urban parameters and calculation of the view factor
       call icBEP(fww,fwg,fgw,fsw,fws,fsg,                         & 
                 z0g_u,z0r_u,                                      & 
                 nd_u,strd_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,          & 
                 nz_u,z_u)                                  

       first=.false.

      endif ! first
      
      do ix=its,ite
      do iy=jts,jte
        if (FRC_URB2D(ix,iy).gt.0.) then    ! Calling BEP only for existing urban classes.
	
         iurb=UTYPE_URB2D(ix,iy)
         ibui=0
         nlev=0
         nbui=0
         d_urb=0.
         do iz=1,nz_um		   
         if(ss_u(iz,iurb).gt.0) then		
           ibui=ibui+1		                
           nlev(ibui)=iz-1
           d_urb(ibui)=ss_u(iz,iurb)
           nbui=ibui
	 endif	  
         end do  !iz
         if (nbui.gt.nbui_max) then
            write (*,*) 'nbui_max must be increased to',nbui
            stop
         endif

       do iz= kts,kte
          ua1D(iz)=u_phy(ix,iz,iy)
          va1D(iz)=v_phy(ix,iz,iy)
	  pt1D(iz)=th_phy(ix,iz,iy)
	  da1D(iz)=rho(ix,iz,iy)
	  pr1D(iz)=p_phy(ix,iz,iy)
!	  pt01D(iz)=th_phy(ix,iz,iy)
	  pt01D(iz)=300.
	  z1D(iz)=z(ix,iz,iy)
          qv1D(iz)=qv_phy(ix,iz,iy)
          a_u1D(iz)=0.
          a_v1D(iz)=0.
          a_t1D(iz)=0.
          a_e1D(iz)=0.
          b_u1D(iz)=0.
          b_v1D(iz)=0.
          b_t1D(iz)=0.
          b_ac1D(iz)=0.
          b_e1D(iz)=0.           
         enddo
	 z1D(kte+1)=z(ix,kte+1,iy)

         do id=1,ndm
         do iz_u=1,nz_um
         do iw=1,nwr_u
         do ibui=1,nbui_max
          tw1D(2*id-1,iz_u,iw,ibui)=tw1_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)
          tw1D(2*id,iz_u,iw,ibui)=tw2_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)
         enddo
         enddo
         enddo
         enddo
	
         do id=1,ndm
          do ig=1,ng_u
            tg1D(id,ig)=tgb_urb4d(ix,ind_gd(ig,id),iy)
          enddo
          do iz_u=1,nz_um
          do ir=1,nwr_u
            tr1D(id,iz_u,ir)=trb_urb4d(ix,ind_zrd(iz_u,ir,id),iy)
          enddo
          enddo
         enddo
!
!Initialize variables for BEM
!
         tlev1D=0.  !Indoor temperature
         qlev1D=0.  !Indoor humidity

         twlev1D=0. !Window temperature
         tglev1D=0. !Ground temperature
         tflev1D=0. !Floor temperature

         sflev1D=0.    !Sensible heat flux from the a.c.
         lflev1D=0.    !latent heat flux from the a.c.
         consumlev1D=0.!consumption of the a.c.
         sfvlev1D=0.   !Sensible heat flux from natural ventilation
         lfvlev1D=0.   !Latent heat flux from natural ventilation
         sfwin1D=0.    !Sensible heat flux from windows
         sfw1D=0.      !Sensible heat flux from walls         

         do iz_u=1,nz_um    !vertical levels
         do ibui=1,nbui_max !Type of building
            tlev1D(iz_u,ibui)= tlev_urb3d(ix,ind_bd(ibui,iz_u),iy)  
            qlev1D(iz_u,ibui)= qlev_urb3d(ix,ind_bd(ibui,iz_u),iy)  
         enddo !ibui
         enddo !iz_u

         do id=1,ndm  !direction
            do iz_u=1,nz_um !vertical levels
               do ibui=1,nbui_max !type of building
                  twlev1D(2*id-1,iz_u,ibui)=tw1lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
                  twlev1D(2*id,iz_u,ibui)=tw2lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
                  sfwin1D(2*id-1,iz_u,ibui)=sfwin1_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
                  sfwin1D(2*id,iz_u,ibui)=sfwin2_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
               enddo !ibui  
            enddo !iz_u
         enddo !id

         do id=1,ndm !direction
            do iw=1,ngb_u !layer in the wall
               do ibui=1,nbui_max !type of building
                  tglev1D(id,iw,ibui)=tglev_urb3d(ix,ind_gbd(ibui,iw,id),iy)
               enddo !ibui
            enddo !iw
         enddo !id
       
         do id=1,ndm !direction
            do iw=1,nf_u !layer in the walls
               do iz_u=1,nz_um-1 !verticals levels
                  do ibui=1,nbui_max !type of building
                     tflev1D(id,iw,iz_u,ibui)=tflev_urb3d(ix,ind_fbd(ibui,iw,iz_u,id),iy)
                  enddo !ibui
               enddo ! iz_u
             enddo !iw
         enddo !id

!
!End initialization for BEM
!         

         do id=1,ndm
	 do iz=1,nz_um
         do ibui=1,nbui_max !type of building
	  sfw1D(2*id-1,iz,ibui)=sfw1_urb3d(ix,ind_zd(ibui,iz,id),iy)
	  sfw1D(2*id,iz,ibui)=sfw2_urb3d(ix,ind_zd(ibui,iz,id),iy)
         enddo
	 enddo
	 enddo
	 
	 do id=1,ndm
	  sfg1D(id)=sfg_urb3d(ix,id,iy)
	 enddo
	 
	 do id=1,ndm
	 do iz=1,nz_um
	  sfr1D(id,iz)=sfr_urb3d(ix,ind_zdf(iz,id),iy)
	 enddo
	 enddo
         
         rs1D=swdown(ix,iy)
         rld1D=glw(ix,iy)

         zr1D=acos(COSZ_URB2D(ix,iy))
         deltar1D=DECLIN_URB
         ah1D=OMG_URB2D(ix,iy)       

         call BEP1D(iurb,kms,kme,kts,kte,z1D,dt,ua1D,va1D,pt1D,da1D,pr1D,pt01D,  &
                   zr1D,deltar1D,ah1D,rs1D,rld1D,                   & 
                   alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,          & 
                   albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,emr_u, & 
                   emwind_u,fww,fwg,fgw,fsw,fws,fsg,                & 
                   z0g_u,z0r_u,                                     & 
                   nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,  & 
                   nz_u,z_u,                                        & 
                   cop_u,pwin_u,beta_u,sw_cond_u,time_on_u,         &
                   time_off_u,targtemp_u,gaptemp_u,targhum_u,       &
                   gaphum_u, perflo_u, hsesf_u, hsequip,            &
                   tw1D,tg1D,tr1D,sfw1D,sfg1D,sfr1D,                & 
                   a_u1D,a_v1D,a_t1D,a_e1D,                         & 
                   b_u1D,b_v1D,b_t1D,b_ac1D,b_e1D,b_q1D,            & 
                   dlg1D,dl_u1D,sf1D,vl1D,rl_up(ix,iy),             &
                   rs_abs(ix,iy),emiss(ix,iy),grdflx_urb(ix,iy),    &
                   qv1D,tlev1D,qlev1D,sflev1D,lflev1D,consumlev1D,  &
                   sfvlev1D,lfvlev1D,twlev1D,tglev1D,tflev1D,sfwin1D,&
                   ix,iy)                            

         do id=1,ndm ! direction 
	    do iz=1,nz_um   !vertical levels
               do ibui=1,nbui_max !type of building
	          sfw1_urb3d(ix,ind_zd(ibui,iz,id),iy)=sfw1D(2*id-1,iz,ibui) 
	          sfw2_urb3d(ix,ind_zd(ibui,iz,id),iy)=sfw1D(2*id,iz,ibui) 
	       enddo
	    enddo
         enddo
	 do id=1,ndm
	  sfg_urb3d(ix,id,iy)=sfg1D(id) 
	 enddo
         
	 do id=1,ndm
	 do iz=1,nz_um
	  sfr_urb3d(ix,ind_zdf(iz,id),iy)=sfr1D(id,iz)
	 enddo
	 enddo
         
         do id=1,ndm
         do iz_u=1,nz_um
         do iw=1,nwr_u
         do ibui=1,nbui_max
          tw1_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)=tw1D(2*id-1,iz_u,iw,ibui)
          tw2_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)=tw1D(2*id,iz_u,iw,ibui)
         enddo
         enddo
         enddo
         enddo
         
         do id=1,ndm
            do ig=1,ng_u
               tgb_urb4d(ix,ind_gd(ig,id),iy)=tg1D(id,ig)
            enddo
            do iz_u=1,nz_um
               do ir=1,nwr_u
                  trb_urb4d(ix,ind_zrd(iz_u,ir,id),iy)=tr1D(id,iz_u,ir)
               enddo
            enddo
         enddo
!
!Outputs of BEM
!
         do ibui=1,nbui_max !type of building
         do iz_u=1,nz_um !vertical levels
            tlev_urb3d(ix,ind_bd(ibui,iz_u),iy)=tlev1D(iz_u,ibui)  
            qlev_urb3d(ix,ind_bd(ibui,iz_u),iy)=qlev1D(iz_u,ibui)  
         enddo !iz_u
         enddo !ibui
         do id=1,ndm !direction
         do iz_u=1,nz_um !vertical levels
         do ibui=1,nbui_max !type of building
               tw1lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=twlev1D(2*id-1,iz_u,ibui)
               tw2lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=twlev1D(2*id,iz_u,ibui)
               sfwin1_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=sfwin1D(2*id-1,iz_u,ibui)
               sfwin2_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=sfwin1D(2*id,iz_u,ibui)
            enddo !ibui  
         enddo !iz_u
         enddo !id
        
         do id=1,ndm !direction
            do iw=1,ngb_u !layers in the walls
               do ibui=1,nbui_max  !type of building
                  tglev_urb3d(ix,ind_gbd(ibui,iw,id),iy)=tglev1D(id,iw,ibui)
               enddo !ibui
            enddo !iw
         enddo !id

         do id=1,ndm   !direction 
            do iw=1,nf_u !layer in the walls
               do iz_u=1,nz_um-1 !verticals levels
                  do ibui=1,nbui_max  !type of building
                     tflev_urb3d(ix,ind_fbd(ibui,iw,iz_u,id),iy)=tflev1D(id,iw,iz_u,ibui)
                  enddo !ibui
               enddo !iz_u
             enddo !iw
         enddo !id
         sf_ac_urb3d(ix,iy)=0.
         lf_ac_urb3d(ix,iy)=0.
         cm_ac_urb3d(ix,iy)=0.
         sfvent_urb3d(ix,iy)=0.
         lfvent_urb3d(ix,iy)=0.

         meso_urb=(1./4.)*FRC_URB2D(ix,iy)/((bs_u(1,iurb)+ws_u(1,iurb))*bs_u(2,iurb))+ &
                  (1./4.)*FRC_URB2D(ix,iy)/((bs_u(2,iurb)+ws_u(2,iurb))*bs_u(1,iurb))

                  
         ibui=0
         nlev=0
         nbui=0
         d_urb=0.
         do iz=1,nz_um		   
         if(ss_u(iz,iurb).gt.0) then		
           ibui=ibui+1		                
           nlev(ibui)=iz-1
           d_urb(ibui)=ss_u(iz,iurb)
           nbui=ibui
	 endif	  
         end do  !iz

         do ibui=1,nbui       !type of building   
         do iz_u=1,nlev(ibui) !vertical levels                    
               sf_ac_urb3d(ix,iy)=sf_ac_urb3d(ix,iy)+meso_urb*d_urb(ibui)*sflev1D(iz_u,ibui)
               lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)+meso_urb*d_urb(ibui)*lflev1D(iz_u,ibui)
               cm_ac_urb3d(ix,iy)=cm_ac_urb3d(ix,iy)+meso_urb*d_urb(ibui)*consumlev1D(iz_u,ibui)
               sfvent_urb3d(ix,iy)=sfvent_urb3d(ix,iy)+meso_urb*d_urb(ibui)*sfvlev1D(iz_u,ibui)
               lfvent_urb3d(ix,iy)=lfvent_urb3d(ix,iy)+meso_urb*d_urb(ibui)*lfvlev1D(iz_u,ibui)            
         enddo !iz_u
         enddo !ibui
!
!Add the latent heat exchanged throughout the ventilation in the lf_ac_urb3d output variable. 
!it is only a print variable
!
!        lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)+lfvent_urb3d(ix,iy)
!

         lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)-lfvent_urb3d(ix,iy)        

!
!End outputs of bem
!

        sf_ac=0.
        do iz= kts,kte
          sf(ix,iz,iy)=sf1D(iz)
          vl(ix,iz,iy)=vl1D(iz)
          a_u(ix,iz,iy)=a_u1D(iz)
          a_v(ix,iz,iy)=a_v1D(iz)
          a_t(ix,iz,iy)=a_t1D(iz)
          a_e(ix,iz,iy)=a_e1D(iz)
          b_u(ix,iz,iy)=b_u1D(iz)
          b_v(ix,iz,iy)=b_v1D(iz)
          b_t(ix,iz,iy)=b_t1D(iz)
          sf_ac=sf_ac+b_ac1D(iz)*da1D(iz)*cp_u*dz8w(ix,iz,iy)*vl1D(iz)*FRC_URB2D(ix,iy)
          b_e(ix,iz,iy)=b_e1D(iz)
          b_q(ix,iz,iy)=b_q1D(iz)
          dlg(ix,iz,iy)=dlg1D(iz)
          dl_u(ix,iz,iy)=dl_u1D(iz)
         enddo
         sf(ix,kte+1,iy)=sf1D(kte+1)

         endif ! FRC_URB2D
   
      enddo  ! iy
      enddo  ! ix

        time_bep=time_bep+dt

      return
      end subroutine BEP_BEM
            
! ===6=8===============================================================72

      subroutine BEP1D(iurb,kms,kme,kts,kte,z,dt,ua,va,pt,da,pr,pt0,   &  
                      zr,deltar,ah,rs,rld,                             & 
                      alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,          & 
                      albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,emr_u, & 
                      emwind_u,fww,fwg,fgw,fsw,fws,fsg,                & 
                      z0g_u,z0r_u,                                     & 
                      nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,  & 
                      nz_u,z_u,                                        & 
                      cop_u,pwin_u,beta_u,sw_cond_u,time_on_u,         &
                      time_off_u,targtemp_u,gaptemp_u,targhum_u,       &
                      gaphum_u, perflo_u, hsesf_u, hsequip,            &
                      tw,tg,tr,sfw,sfg,sfr,                            & 
                      a_u,a_v,a_t,a_e,                                 &
                      b_u,b_v,b_t,b_ac,b_e,b_q,                        & 
                      dlg,dl_u,sf,vl,rl_up,rs_abs,emiss,grdflx_urb,    &
                      qv,tlev,qlev,sflev,lflev,consumlev,              &
                      sfvlev,lfvlev,twlev,tglev,tflev,sfwin,ix,iy)                             

! ----------------------------------------------------------------------
! This routine computes the effects of buildings on momentum, heat and
!  TKE (turbulent kinetic energy) sources or sinks and on the mixing length.
! It provides momentum, heat and TKE sources or sinks at different levels of a
!  mesoscale grid defined by the altitude of its cell interfaces "z" and
!  its number of levels "nz".
! The meteorological input parameters (wind, temperature, solar radiation)
!  are specified on the "mesoscale grid".
! The inputs concerning the building and street charateristics are defined
!  on a "urban grid". The "urban grid" is defined with its number of levels
!  "nz_u" and its space step "dz_u".
! The input parameters are interpolated on the "urban grid". The sources or sinks
!  are calculated on the "urban grid". Finally the sources or sinks are 
!  interpolated on the "mesoscale grid".
 

!  Mesoscale grid            Urban grid             Mesoscale grid
!  
! z(4)    ---                                               ---
!          |                                                 |
!          |                                                 |
!          |   Interpolation                  Interpolation  |
!          |            Sources or sinks calculation         |
! z(3)    ---                                               ---
!          | ua               ua_u  ---  uv_a         a_u    |
!          | va               va_u   |   uv_b         b_u    |
!          | pt               pt_u  ---  uh_b         a_v    |
! z(2)    ---                        |    etc...      etc...---
!          |                 z_u(1) ---                      |
!          |                         |                       |
! z(1) ------------------------------------------------------------

!     
! Reference:
! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
! 261-304
 
! ----------------------------------------------------------------------

      implicit none

! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------

! Data relative to the "mesoscale grid"

      integer kms,kme,kts,kte
      real z(kms:kme)               ! Altitude above the ground of the cell interfaces.
      real ua(kms:kme)                ! Wind speed in the x direction
      real va(kms:kme)                ! Wind speed in the y direction
      real pt(kms:kme)                ! Potential temperature
      real da(kms:kme)                ! Air density
      real pr(kms:kme)                ! Air pressure
      real pt0(kms:kme)               ! Reference potential temperature (could be equal to "pt")
      real qv(kms:kme)              ! Specific humidity
      real dt                    ! Time step
      real zr                    ! Zenith angle
      real deltar                ! Declination of the sun
      real ah                    ! Hour angle
      real rs                    ! Solar radiation
      real rld                   ! Downward flux of the longwave radiation

! Data relative to the "urban grid"

      integer iurb               ! Current urban class

!    Building parameters      
      real alag_u(nurbm)         ! Ground thermal diffusivity [m^2 s^-1]
      real alaw_u(nurbm)         ! Wall thermal diffusivity [m^2 s^-1]
      real alar_u(nurbm)         ! Roof thermal diffusivity [m^2 s^-1]
      real csg_u(nurbm)          ! Specific heat of the ground material [J m^3 K^-1]
      real csw_u(nurbm)          ! Specific heat of the wall material [J m^3 K^-1]
      real csr_u(nurbm)          ! Specific heat of the roof material [J m^3 K^-1]

!    Radiation parameters
      real albg_u(nurbm)         ! Albedo of the ground
      real albw_u(nurbm)         ! Albedo of the wall
      real albr_u(nurbm)         ! Albedo of the roof
      real albwin_u(nurbm)       ! Albedo of the windows
      real emwind_u(nurbm)       ! Emissivity of windows
      real emg_u(nurbm)          ! Emissivity of ground
      real emw_u(nurbm)          ! Emissivity of wall
      real emr_u(nurbm)          ! Emissivity of roof

!    fww,fwg,fgw,fsw,fsg are the view factors used to compute the long and 
!    short wave radation. 
!    The calculation of these factor is explained in the Appendix A of the BLM paper
      real fww(nz_um,nz_um,ndm,nurbm)  !  from wall to wall
      real fwg(nz_um,ndm,nurbm)        !  from wall to ground
      real fgw(nz_um,ndm,nurbm)        !  from ground to wall
      real fsw(nz_um,ndm,nurbm)        !  from sky to wall
      real fws(nz_um,ndm,nurbm)        !  from wall to sky
      real fsg(ndm,nurbm)              !  from sky to ground

!    Roughness parameters
      real z0g_u(nurbm)            ! The ground's roughness length
      real z0r_u(nurbm)            ! The roof's roughness length
      
!    Street parameters
      integer nd_u(nurbm)          ! Number of street direction for each urban class 
      real strd_u(ndm,nurbm)       ! Street length (set to a greater value then the horizontal length of the cells)
      real drst_u(ndm,nurbm)       ! Street direction
      real ws_u(ndm,nurbm)         ! Street width
      real bs_u(ndm,nurbm)         ! Building width
      real h_b(nz_um,nurbm)        ! Bulding's heights
      real d_b(nz_um,nurbm)        ! The probability that a building has an height "h_b"
      real ss_u(nz_um,nurbm)     ! The probability that a building has an height equal to "z"
      real pb_u(nz_um,nurbm)     ! The probability that a building has an height greater or equal to "z"
        
!    Grid parameters
      integer nz_u(nurbm)     ! Number of layer in the urban grid
      real z_u(nz_um)       ! Height of the urban grid levels

!   MT
      real cop_u(nurbm)
      real pwin_u(nurbm)
      real beta_u(nurbm)
      integer sw_cond_u(nurbm)
      real time_on_u(nurbm)
      real time_off_u(nurbm)
      real targtemp_u(nurbm)
      real gaptemp_u(nurbm)
      real targhum_u(nurbm)
      real gaphum_u(nurbm)
      real perflo_u(nurbm)
      real hsesf_u(nurbm)
      real hsequip(24)

! ----------------------------------------------------------------------
! INPUT-OUTPUT
! ----------------------------------------------------------------------

! Data relative to the "urban grid" which should be stored from the current time step to the next one

      real tw(2*ndm,nz_um,nwr_u,nbui_max)  ! Temperature in each layer of the wall [K]
      real tr(ndm,nz_um,nwr_u)  ! Temperature in each layer of the roof [K]
      real tg(ndm,ng_u)          ! Temperature in each layer of the ground [K]
      real sfw(2*ndm,nz_um,nbui_max)      ! Sensible heat flux from walls
      real sfg(ndm)              ! Sensible heat flux from ground (road)
      real sfr(ndm,nz_um)      ! Sensible heat flux from roofs
      real gfg(ndm)             ! Heat flux transferred from the surface of the ground (road) towards the interior
      real gfr(ndm,nz_um)     ! Heat flux transferred from the surface of the roof towards the interior
      real gfw(2*ndm,nz_um,nbui_max)     ! Heat flux transfered from the surface of the walls towards the interior
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------

! Data relative to the "mesoscale grid"

      real sf(kms:kme)             ! Surface of the "mesoscale grid" cells taking into account the buildings
      real vl(kms:kme)               ! Volume of the "mesoscale grid" cells taking into account the buildings
     
!    Implicit and explicit components of the source and sink terms at each levels,
!     the fluxes can be computed as follow: FX = A*X + B   example: Heat fluxes = a_t * pt + b_t
      real a_u(kms:kme)              ! Implicit component of the momentum sources or sinks in the X-direction
      real a_v(kms:kme)              ! Implicit component of the momentum sources or sinks in the Y-direction
      real a_t(kms:kme)              ! Implicit component of the heat sources or sinks
      real a_e(kms:kme)              ! Implicit component of the TKE sources or sinks
      real b_u(kms:kme)              ! Explicit component of the momentum sources or sinks in the X-direction
      real b_v(kms:kme)              ! Explicit component of the momentum sources or sinks in the Y-direction
      real b_t(kms:kme)              ! Explicit component of the heat sources or sinks
      real b_ac(kms:kme)
      real b_e(kms:kme)              ! Explicit component of the TKE sources or sinks
      real b_q(kms:kme)              ! Explicit component of the humidity sources or sinks
      real dlg(kms:kme)              ! Height above ground (L_ground in formula (24) of the BLM paper). 
      real dl_u(kms:kme)             ! Length scale (lb in formula (22) ofthe BLM paper).
      
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------

      real dz(kms:kme)               ! vertical space steps of the "mesoscale grid"

! Data interpolated from the "mesoscale grid" to the "urban grid"

      real ua_u(nz_um)          ! Wind speed in the x direction
      real va_u(nz_um)          ! Wind speed in the y direction
      real pt_u(nz_um)          ! Potential temperature
      real da_u(nz_um)          ! Air density
      real pt0_u(nz_um)         ! Reference potential temperature
      real pr_u(nz_um)          ! Air pressure
      real qv_u(nz_um)          !Specific humidity

! Data defining the building and street charateristics

      real alag(ng_u)           ! Ground thermal diffusivity for the current urban class [m^2 s^-1] 
      
      real csg(ng_u)            ! Specific heat of the ground material of the current urban class [J m^3 K^-1]
      real csr(nwr_u)            ! Specific heat of the roof material for the current urban class [J m^3 K^-1]
      real csw(nwr_u)            ! Specific heat of the wall material for the current urban class [J m^3 K^-1]

      real z0(ndm,nz_um)      ! Roughness lengths "profiles"
      real ws(ndm)              ! Street widths of the current urban class
      real bs(ndm)              ! Building widths of the current urban class
      real strd(ndm)            ! Street lengths for the current urban class
      real drst(ndm)            ! Street directions for the current urban class
      real ss(nz_um)          ! Probability to have a building with height h
      real pb(nz_um)          ! Probability to have a building with an height equal

! Solar radiation at each level of the "urban grid"

      real rsg(ndm)             ! Short wave radiation from the ground
      real rsw(2*ndm,nz_um)     ! Short wave radiation from the walls
      real rlg(ndm)             ! Long wave radiation from the ground
      real rlw(2*ndm,nz_um)     ! Long wave radiation from the walls

! Potential temperature of the surfaces at each level of the "urban grid"

      real ptg(ndm)             ! Ground potential temperatures 
      real ptr(ndm,nz_um)     ! Roof potential temperatures 
      real ptw(2*ndm,nz_um,nbui_max)     ! Walls potential temperatures 

 
! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
! vertical surfaces (walls) ans horizontal surfaces (roofs and street)
! The fluxes can be computed as follow: Fluxes of X = A*X + B
! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u

      real uhb_u(ndm,nz_um)   ! U (wind component) Horizontal surfaces, B (explicit) term
      real uva_u(2*ndm,nz_um)   ! U (wind component)   Vertical surfaces, A (implicit) term
      real uvb_u(2*ndm,nz_um)   ! U (wind component)   Vertical surfaces, B (explicit) term
      real vhb_u(ndm,nz_um)   ! V (wind component) Horizontal surfaces, B (explicit) term
      real vva_u(2*ndm,nz_um)   ! V (wind component)   Vertical surfaces, A (implicit) term
      real vvb_u(2*ndm,nz_um)   ! V (wind component)   Vertical surfaces, B (explicit) term
      real thb_u(ndm,nz_um)   ! Temperature        Horizontal surfaces, B (explicit) term
      real tva_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, A (implicit) term
      real tvb_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, B (explicit) term
      real tvb_ac(2*ndm,nz_um)
      real ehb_u(ndm,nz_um)   ! Energy (TKE)       Horizontal surfaces, B (explicit) term
      real evb_u(2*ndm,nz_um)   ! Energy (TKE)         Vertical surfaces, B (explicit) term
      real qhb_u(ndm,nz_um)     ! Humidity      Horizontal surfaces, B (explicit) term
      real qvb_u(2*ndm,nz_um)   ! Humidity      Vertical surfaces, B (explicit) term      
!
      real rs_abs ! solar radiation absorbed by urban surfaces 
      real rl_up ! longwave radiation emitted by urban surface to the atmosphere 
      real emiss ! mean emissivity of the urban surface
      real grdflx_urb ! ground heat flux
      real dt_int ! internal time step
      integer nt_int ! number of internal time step
      integer iz,id, it_int
      integer iw,ix,iy

!---------------------------------------
!New variables uses in BEM
!----------------------------------------
   
      real tmp_u(nz_um)     !Air Temperature [K]

      real dzw(nwr_u)       !Layer sizes in the walls
      real dzr(nwr_u)       !Layer sizes in the roofs
      real dzf(nf_u)        !Layer sizes in the floors
      real dzgb(ngb_u)      !Layer sizes in the ground below the buildings

      real csgb(ngb_u)      !Specific heat of the ground material below the buildings 
                            !of the current urban class at each ground levels[J m^3 K^-1]
      real csf(nf_u)        !Specific heat of the floors materials in the buildings 
                            !of the current urban class at each levels[J m^3 K^-1]
      real alar(nwr_u+1)    ! Roof thermal diffusivity for the current urban class [W/m K]
      real alaw(nwr_u+1)    ! Walls thermal diffusivity for the current urban class [W/m K] 
      real alaf(nf_u+1)     ! Floor thermal diffusivity at each wall layers [W/m K]     
      real alagb(ngb_u+1)   ! Ground thermal diffusivity below the building at each wall layer [W/m K] 

      real sfrb(ndm,nbui_max)        ! Sensible heat flux from roofs [W/mÂ²]
      real gfrb(ndm,nbui_max)        ! Heat flux flowing inside the roofs [W/mÂ²]
      real sfwb1D(2*ndm,nz_um)    !Sensible heat flux from the walls [W/mÂ²] 
      real sfwin(2*ndm,nz_um,nbui_max)!Sensible heat flux from windows [W/mÂ²]
      real sfwinb1D(2*ndm,nz_um)  !Sensible heat flux from windows [W/mÂ²]
      real gfwb1D(2*ndm,nz_um)    !Heat flux flowing inside the walls [W/mÂ²]

      real qlev(nz_um,nbui_max)      !specific humidity [kg/kg]
      real qlevb1D(nz_um)         !specific humidity [kg/kg] 
      real tlev(nz_um,nbui_max)      !Indoor temperature [K]
      real tlevb1D(nz_um)         !Indoor temperature [K]
      real twb1D(2*ndm,nwr_u,nz_um)     !Wall temperature in BEM [K]      
      real twlev(2*ndm,nz_um,nbui_max)     !Window temperature in BEM [K]
      real twlevb1D(2*ndm,nz_um)        !Window temperature in BEM [K]
      real tglev(ndm,ngb_u,nbui_max)        !Ground temperature below a building in BEM [K]
      real tglevb1D(ngb_u)               !Ground temperature below a building in BEM [K]
      real tflev(ndm,nf_u,nz_um-1,nbui_max)!Floor temperature in BEM[K]
      real tflevb1D(nf_u,nz_um-1)       !Floor temperature in BEM[K]
      real trb(ndm,nwr_u,nbui_max)         !Roof temperature in BEM [K]
      real trb1D(nwr_u)                 !Roof temperature in BEM [K]
      
      real sflev(nz_um,nz_um)     ! sensible heat flux due to the air conditioning systems [W]
      real lflev(nz_um,nz_um)     ! latent heat flux due to the air conditioning systems  [W]
      real consumlev(nz_um,nz_um) ! consumption due to the air conditioning systems [W]
      real sflev1D(nz_um)         ! sensible heat flux due to the air conditioning systems [W]
      real lflev1D(nz_um)         ! latent heat flux due to the air conditioning systems  [W]
      real consumlev1D(nz_um)     ! consumption due to the air conditioning systems [W]
      real sfvlev(nz_um,nz_um)    ! sensible heat flux due to ventilation [W]
      real lfvlev(nz_um,nz_um)    ! latent heat flux due to ventilation [W]
      real sfvlev1D(nz_um)        ! sensible heat flux due to ventilation [W]
      real lfvlev1D(nz_um)        ! Latent heat flux due to ventilation [W]

      real ptwin(2*ndm,nz_um,nbui_max)  ! window potential temperature
      real tw_av(2*ndm,nz_um)        ! Averaged temperature of the wall surfaces
      real twlev_av(2*ndm,nz_um)     ! Averaged temperature of the windows
      real sfw_av(2*ndm,nz_um)       ! Averaged sensible heat from walls
      real sfwind_av(2*ndm,nz_um)    ! Averaged sensible heat from windows

      integer nbui                !Total number of different type of buildings in an urban class
      integer nlev(nz_um)         !Number of levels in each different type of buildings in an urban class
      integer ibui,ily  
      real :: nhourday   ! Number of hours from midnight, local time
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------

! Fix some usefull parameters for the computation of the sources or sinks
!
!initialize inside param
!
!      ss=0.
!      pb=0.

      do iz=kts,kte
         dz(iz)=z(iz+1)-z(iz)
      end do
      call param(iurb,nz_u(iurb),nd_u(iurb),                          &
                       csg_u,csg,alag_u,alag,csr_u,csr,               &
                       alar_u,alar,csw_u,csw,alaw_u,alaw,             &
                       ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0,                &
                       strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb,       &
                       dzw,dzr,dzf,csf,alaf,dzgb,csgb,alagb)

! Interpolation on the "urban grid"

      call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,ua,ua_u)
      call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,va,va_u)
      call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt,pt_u)
      call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt0,pt0_u)
      call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pr,pr_u)
      call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,da,da_u)
      call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,qv,qv_u)
                   
! Compute the modification of the radiation due to the buildings

      call averaging_temp(tw,twlev,ss,pb,tw_av,twlev_av, &
                          sfw_av,sfwind_av,sfw,sfwin)
      call modif_rad(iurb,nd_u(iurb),nz_u(iurb),z_u,ws,           &
                    drst,strd,ss,pb,                              &
                    tw_av,tg,twlev_av,albg_u(iurb),albw_u(iurb),  &
                    emw_u(iurb),emg_u(iurb),pwin_u(iurb),albwin_u(iurb),  &
                    emwind_u(iurb),fww,fwg,fgw,fsw,fsg,           &
                    zr,deltar,ah,                                 &
                    rs,rld,rsw,rsg,rlw,rlg)                       

! calculation of the urban albedo and the upward long wave radiation

       call upward_rad(nd_u(iurb),nz_u(iurb),ws,bs,sigma,pb,ss,                &
                       tg,emg_u(iurb),albg_u(iurb),rlg,rsg,sfg,                & 
                       tw_av,emw_u(iurb),albw_u(iurb),rlw,rsw,sfw_av,          & 
                       tr,emr_u(iurb),albr_u(iurb),emwind_u(iurb),             &
                       albwin_u(iurb),twlev_av,pwin_u(iurb),sfwind_av,rld,rs,sfr,      & 
                       rs_abs,rl_up,emiss,grdflx_urb)               

! Compute the surface temperatures
      
      call surf_temp(nd_u(iurb),pr_u,dt,              & 
                    rld,rsg,rlg,                      &
                    tg,alag,csg,emg_u(iurb),albg_u(iurb),ptg,sfg,gfg)


! Call the BEM (Building Energy Model) routine 
       
       do iz=1,nz_um !Compute the outdoor temperature 
	 tmp_u(iz)=pt_u(iz)*(pr_u(iz)/p0)**(rcp_u) 
       end do

       ibui=0
       nlev=0
       nbui=0
     
       sfrb=0.     !Sensible heat flux from roof
       gfrb=0.     !Heat flux flowing inside the roof
       sfwb1D=0.   !Sensible heat flux from walls
       sfwinb1D=0. !Sensible heat flux from windows
       gfwb1D=0.   !Heat flux flowing inside the walls[W/mÂ²]


       twb1D=0.    !Wall temperature
       twlevb1D=0. !Window temperature
       tglevb1D=0. !Ground temperature below a building
       tflevb1D=0. !Floor temperature     
       trb=0.      !Roof temperature
       trb1D=0.    !Roof temperature

       qlevb1D=0. !Indoor humidity
       tlevb1D=0. !indoor temperature

       sflev1D=0.    !Sensible heat flux from the a.c.
       lflev1D=0.    !Latent heat flux from the a.c.
       consumlev1D=0.!Consumption from the a.c.
       sfvlev1D=0.   !Sensible heat flux from the natural ventilation
       lfvlev1D=0.   !Latent heat flux from natural ventilation

       ptw=0.        !Wall potential temperature
       ptwin=0.      !Window potential temperature
       ptr=0.        !Roof potential temperature
       
       do iz=1,nz_um		   
         if(ss(iz).gt.0) then		
           ibui=ibui+1		                
           nlev(ibui)=iz-1
           nbui=ibui
           do id=1,ndm
              sfrb(id,ibui)=sfr(id,iz)
              do ily=1,nwr_u
                 trb(id,ily,ibui)=tr(id,iz,ily)
              enddo
           enddo
	 endif	  
       end do  !iz    
!--------------------------------------------------------------------------------
!Loop over BEM  -----------------------------------------------------------------
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
   
       nhourday=ah/PI*180./15.+12.
       if (nhourday >= 24) nhourday = nhourday - 24
       if (nhourday < 0)  nhourday = nhourday + 24

        do ibui=1,nbui
         
          do iz=1,nz_um
             qlevb1D(iz)=qlev(iz,ibui)
             tlevb1D(iz)=tlev(iz,ibui) 
          enddo
             
          do id=1,ndm

             do ily=1,nwr_u
                trb1D(ily)=trb(id,ily,ibui)
             enddo
             do ily=1,ngb_u
                tglevb1D(ily)=tglev(id,ily,ibui) 
             enddo
              
             do ily=1,nf_u
             do iz=1,nz_um-1
               tflevb1D(ily,iz)=tflev(id,ily,iz,ibui)
             enddo
             enddo
            
             do iz=1,nz_um
                sfwinb1D(2*id-1,iz)=sfwin(2*id-1,iz,ibui)
                sfwinb1D(2*id,iz)=sfwin(2*id,iz,ibui)
             enddo
            
             do iz=1,nz_um
                do ily=1,nwr_u
                   twb1D(2*id-1,ily,iz)=tw(2*id-1,iz,ily,ibui)
                   twb1D(2*id,ily,iz)=tw(2*id,iz,ily,ibui)
                enddo
                sfwb1D(2*id-1,iz)=sfw(2*id-1,iz,ibui)
                sfwb1D(2*id,iz)=sfw(2*id,iz,ibui)
                twlevb1D(2*id-1,iz)=twlev(2*id-1,iz,ibui)
                twlevb1D(2*id,iz)=twlev(2*id,iz,ibui)
             enddo
          enddo   

          call BEM(nz_um,nlev(ibui),nhourday,dt,bs_u(1,iurb),                   &
                   bs_u(2,iurb),dz_u,nwr_u,nf_u,nwr_u,ngb_u,sfwb1D,gfwb1D,      &
                   sfwinb1D,sfrb(1,ibui),gfrb(1,ibui),                          &
                   latent,sigma,albw_u(iurb),albwin_u(iurb),albr_u(iurb),       &
     	           emr_u(iurb),emw_u(iurb),emwind_u(iurb),rsw,rlw,r,cp_u,       &
     	           da_u,tmp_u,qv_u,pr_u,rs,rld,dzw,csw,alaw,pwin_u(iurb),       &
                   cop_u(iurb),beta_u(iurb),sw_cond_u(iurb),time_on_u(iurb),    &
                   time_off_u(iurb),targtemp_u(iurb),gaptemp_u(iurb),           &
                   targhum_u(iurb),gaphum_u(iurb),perflo_u(iurb),hsesf_u(iurb), &
                   hsequip,                                                     &
     	           dzf,csf,alaf,dzgb,csgb,alagb,dzr,csr,                        &
     	           alar,tlevb1D,qlevb1D,twb1D,twlevb1D,tflevb1D,tglevb1D,       &
     	           trb1D,sflev1D,lflev1D,consumlev1D,sfvlev1D,lfvlev1D)
                    
        
!
!Temporal modifications
!
         sfrb(2,ibui)=sfrb(1,ibui)
         gfrb(2,ibui)=gfrb(1,ibui)
!
!End temporal modifications  
!        
           do iz=1,nz_um
             qlev(iz,ibui)=qlevb1D(iz)
             tlev(iz,ibui)=tlevb1D(iz)
             sflev(iz,ibui)=sflev1D(iz)
             lflev(iz,ibui)=lflev1D(iz)
             consumlev(iz,ibui)=consumlev1D(iz)
             sfvlev(iz,ibui)=sfvlev1D(iz)
             lfvlev(iz,ibui)=lfvlev1D(iz)
           enddo
 
           do id=1,ndm
              do ily=1,nwr_u
                 trb(id,ily,ibui)=trb1D(ily)
              enddo   
              do ily=1,ngb_u
                 tglev(id,ily,ibui)=tglevb1D(ily) 
              enddo

              do ily=1,nf_u
              do iz=1,nz_um-1
                 tflev(id,ily,iz,ibui)=tflevb1D(ily,iz)
              enddo
              enddo
           

             do iz=1,nz_um
                do ily=1,nwr_u
                   tw(2*id-1,iz,ily,ibui)=twb1D(2*id-1,ily,iz)
                   tw(2*id,iz,ily,ibui)=twb1D(2*id,ily,iz)
                enddo
                gfw(2*id-1,iz,ibui)=gfwb1D(2*id-1,iz)
                gfw(2*id,iz,ibui)=gfwb1D(2*id,iz)
                twlev(2*id-1,iz,ibui)=twlevb1D(2*id-1,iz)
                twlev(2*id,iz,ibui)=twlevb1D(2*id,iz)
             enddo
           enddo       

        enddo !ibui     

!-----------------------------------------------------------------------------
!End loop over BEM -----------------------------------------------------------
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------

       ibui=0
       do iz=1,nz_um		   
         if(ss(iz).gt.0) then		
           ibui=ibui+1	
           do id=1,ndm	
              gfr(id,iz)=gfrb(id,ibui)
              sfr(id,iz)=sfrb(id,ibui)
              do ily=1,nwr_u
                 tr(id,iz,ily)=trb(id,ily,ibui)
              enddo
              ptr(id,iz)=tr(id,iz,nwr_u)*(pr_u(iz)/p0)**(-rcp_u)
           enddo
         endif
       enddo !iz

!Compute the potential temperature for the vertical surfaces of the buildings
      
       do id=1,ndm
          do iz=1,nz_um
             do ibui=1,nbui
                ptw(2*id-1,iz,ibui)=tw(2*id-1,iz,nwr_u,ibui)*(pr_u(iz)/p0)**(-rcp_u) 
                ptw(2*id,iz,ibui)=tw(2*id,iz,nwr_u,ibui)*(pr_u(iz)/p0)**(-rcp_u) 
                ptwin(2*id-1,iz,ibui)=twlev(2*id-1,iz,ibui)*(pr_u(iz)/p0)**(-rcp_u) 
                ptwin(2*id,iz,ibui)=twlev(2*id,iz,ibui)*(pr_u(iz)/p0)**(-rcp_u) 
             enddo
          enddo
       enddo
             
! Compute the implicit and explicit components of the sources or sinks on the "urban grid"
     
      call buildings(iurb,nd_u(iurb),nz_u(iurb),z0,ua_u,va_u,                 & 
                     pt_u,pt0_u,ptg,ptr,da_u,ptw,ptwin,pwin_u(iurb),drst,     &                      
                     uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u,qvb_u,qhb_u,   & 
                     uhb_u,vhb_u,thb_u,ehb_u,ss,dt,sfw,sfg,sfr,               &
                     sfwin,pb,bs_u,dz_u,sflev,lflev,sfvlev,lfvlev,tvb_ac)  
     
! Calculation of the sensible heat fluxes for the ground, the wall and roof
! Sensible Heat Flux = density * Cp_U * ( A* potential temperature + B )
! where A and B are the implicit and explicit components of the heat sources or sinks.

      
! Interpolation on the "mesoscale grid"

      call urban_meso(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z,dz,z_u,pb,ss,bs,ws,sf, & 
                     vl,uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u,                   &
                     uhb_u,vhb_u,thb_u,ehb_u,qhb_u,qvb_u,                            &
                     a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e,b_q,tvb_ac,b_ac)                    
       

! Calculation of the length scale taking into account the buildings effects

      call interp_length(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z_u,z,ss,ws,bs,dlg,dl_u)

      return
      end subroutine BEP1D

! ===6=8===============================================================72
! ===6=8===============================================================72

       subroutine param(iurb,nz,nd,                                   &
                       csg_u,csg,alag_u,alag,csr_u,csr,               &
                       alar_u,alar,csw_u,csw,alaw_u,alaw,             &
                       ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0,                &  
                       strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb,       &
                       dzw,dzr,dzf,csf,alaf,dzgb,csgb,alagb)        

! ----------------------------------------------------------------------
!    This routine prepare some usefull parameters       
! ----------------------------------------------------------------------

      implicit none

  
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      integer iurb                 ! Current urban class
      integer nz                   ! Number of vertical urban levels in the current class
      integer nd                   ! Number of street direction for the current urban class
      real alag_u(nurbm)           ! Ground thermal diffusivity [m^2 s^-1]
      real alar_u(nurbm)           ! Roof thermal diffusivity [m^2 s^-1]
      real alaw_u(nurbm)           ! Wall thermal diffusivity [m^2 s^-1]
      real bs_u(ndm,nurbm)         ! Building width
      real csg_u(nurbm)            ! Specific heat of the ground material [J m^3 K^-1]
      real csr_u(nurbm)            ! Specific heat of the roof material [J m^3 K^-1]
      real csw_u(nurbm)            ! Specific heat of the wall material [J m^3 K^-1]
      real drst_u(ndm,nurbm)       ! Street direction
      real strd_u(ndm,nurbm)       ! Street length 
      real ws_u(ndm,nurbm)         ! Street width
      real z0g_u(nurbm)            ! The ground's roughness length
      real z0r_u(nurbm)            ! The roof's roughness length
      real ss_u(nz_um,nurbm)       ! The probability that a building has an height equal to "z"
      real pb_u(nz_um,nurbm)       ! The probability that a building has an height greater or equal to "z"
     

! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
      real alag(ng_u)           ! Ground thermal diffusivity at each ground levels
      real csg(ng_u)            ! Specific heat of the ground material at each ground levels
      real bs(ndm)              ! Building width for the current urban class
      real drst(ndm)            ! street directions for the current urban class
      real strd(ndm)            ! Street lengths for the current urban class
      real ws(ndm)              ! Street widths of the current urban class
      real z0(ndm,nz_um)      ! Roughness lengths "profiles"
      real ss(nz_um)          ! Probability to have a building with height h
      real pb(nz_um)          ! Probability to have a building with an height greater or equal to "z"

!-----------------------------------------------------------------------------
!INPUT/OUTPUT
!-----------------------------------------------------------------------------

      real dzw(nwr_u)       !Layer sizes in the walls [m]
      real dzr(nwr_u)       !Layer sizes in the roofs [m]
      real dzf(nf_u)        !Layer sizes in the floors [m]
      real dzgb(ngb_u)      !layer sizes in the ground below the buildings [m]

      real csr(nwr_u)       ! Specific heat of the roof material at each roof levels
      real csw(nwr_u)       ! Specific heat of the wall material at each wall levels

      real csf(nf_u)        !Specific heat of the floors materials in the buildings 
                            !of the current urban class [J m^3 K^-1]
      real csgb(ngb_u)      !Specific heat of the ground material below the buildings 
                            !of the current urban class [J m^3 K^-1]
      real alar(nwr_u+1)    ! Roof thermal diffusivity at each roof levels [W/ m K]
      real alaw(nwr_u+1)    ! Wall thermal diffusivity at each wall levels [W/ m K]
      real alaf(nf_u+1)     ! Floor thermal diffusivity at each wall levels [W/m K]
      real alagb(ngb_u+1)   ! Ground thermal diffusivity below the building at each wall levels [W/m K]
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      integer id,ig,ir,iw,iz,iflo
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------  
!Define the layer sizes in the walls
        
      ss=0.
      pb=0.
      csg=0.
      alag=0.
      csr=0.
      alar=0.
      csw=0.
      alaw=0.
      z0=0.
      ws=0.
      bs=0.
      strd=0.
      drst=0.
      csgb=0.
      alagb=0.
      csf=0.
      alaf=0.
 
      dzgb=(/0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/)
      dzr=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/)   
      dzw=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/)
      dzf=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02/)   
      
      do ig=1,ngb_u
        csgb(ig) = csg_u(iurb)
        alagb(ig)= csg_u(iurb)*alag_u(iurb)
      enddo
      alagb(ngb_u+1)= csg_u(iurb)*alag_u(iurb)
      
      do iflo=1,nf_u
        csf(iflo) = csw_u(iurb)
        alaf(iflo)= csw_u(iurb)*alaw_u(iurb) 
      enddo
      alaf(nf_u+1)= csw_u(iurb)*alaw_u(iurb) 
     
      do ir=1,nwr_u
        csr(ir) = csr_u(iurb)
        alar(ir)= csr_u(iurb)*alar_u(iurb)
      enddo
      alar(nwr_u+1)= csr_u(iurb)*alar_u(iurb)

      do iw=1,nwr_u
        csw(iw) = csw_u(iurb)
        alaw(iw)= csw_u(iurb)*alaw_u(iurb)
      enddo
      alaw(nwr_u+1)=csw_u(iurb)*alaw_u(iurb) 

!------------------------------------------------------------------------  
        
      do iz=1,nz+1
         ss(iz)=ss_u(iz,iurb)
         pb(iz)=pb_u(iz,iurb)
      end do
                 
       do ig=1,ng_u
        csg(ig)=csg_u(iurb)
        alag(ig)=alag_u(iurb)
       enddo
       
       do id=1,nd
        z0(id,1)=z0g_u(iurb)
        do iz=2,nz+1
         z0(id,iz)=z0r_u(iurb)
        enddo
       enddo
      
       do id=1,nd
        ws(id)=ws_u(id,iurb)
        bs(id)=bs_u(id,iurb)
        strd(id)=strd_u(id,iurb)
        drst(id)=drst_u(id,iurb)     
       enddo

       
       return
       end subroutine param
       
! ===6=8===============================================================72
! ===6=8===============================================================72

      subroutine interpol(kms,kme,kts,kte,nz_u,z,z_u,c,c_u)

! ----------------------------------------------------------------------
!  This routine interpolate para
!  meters from the "mesoscale grid" to
!  the "urban grid".
!  See p300 Appendix B.1 of the BLM paper.
! ----------------------------------------------------------------------

      implicit none

! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
! Data relative to the "mesoscale grid"
      integer kts,kte,kms,kme            
      real z(kms:kme)          ! Altitude of the cell interface
      real c(kms:kme)            ! Parameter which has to be interpolated
! Data relative to the "urban grid"
      integer nz_u          ! Number of levels
      real z_u(nz_um)      ! Altitude of the cell interface
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
      real c_u(nz_um)        ! Interpolated paramters in the "urban grid"
       
! LOCAL:
! ----------------------------------------------------------------------
      integer iz_u,iz
      real ctot,dz

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------

       do iz_u=1,nz_u
        ctot=0.
        do iz=kts,kte
         dz=max(min(z(iz+1),z_u(iz_u+1))-max(z(iz),z_u(iz_u)),0.)
         ctot=ctot+c(iz)*dz
        enddo
        c_u(iz_u)=ctot/(z_u(iz_u+1)-z_u(iz_u))
       enddo
       
       return
       end subroutine interpol
         
! ===6=8===============================================================72       
! ===6=8===============================================================72    

      subroutine  averaging_temp(tw,twlev,ss,pb,tw_av,twlev_av,&
                                 sfw_av,sfwind_av,sfw,sfwin) 

      implicit none
!
!INPUT VARIABLES
!
      real tw(2*ndm,nz_um,nwr_u,nbui_max)        ! Temperature in each layer of the wall [K]
      real twlev(2*ndm,nz_um,nbui_max)     ! Window temperature in BEM [K]
      real pb(nz_um)                    ! Probability to have a building with an height equal or greater h
      real ss(nz_um)                    ! Probability to have a building with height h
      real sfw(2*ndm,nz_um,nbui_max)             ! Surface fluxes from the walls
      real sfwin(2*ndm,nz_um,nbui_max)     ! Surface fluxes from the windows
!
!OUTPUT VARIABLES
!
      real tw_av(2*ndm,nz_um)           ! Averaged temperature of the wall surfaces
      real twlev_av(2*ndm,nz_um)        ! Averaged temperature of the windows
      real sfw_av(2*ndm,nz_um)          ! Averaged sensible heat from walls
      real sfwind_av(2*ndm,nz_um)       ! Averaged sensible heat from windows
!
!LOCAL VARIABLES
!
      real d_urb(nz_um)    
      integer nlev(nz_um)            
      integer id,iz
      integer nbui,ibui
!
!Initialize Variables
!
      tw_av=0.
      twlev_av=0.
      sfw_av=0.
      sfwind_av=0.
      ibui=0
      nbui=0
      nlev=0
      d_urb=0.

      do iz=1,nz_um		   
         if(ss(iz).gt.0) then		
           ibui=ibui+1
           d_urb(ibui)=ss(iz)
           nlev(ibui)=iz-1
           nbui=ibui		               
         endif
      enddo

      do id=1,ndm
         do iz=1,nz_um-1
            if (pb(iz+1).gt.0) then
                do ibui=1,nbui
                   if (iz.le.nlev(ibui)) then
                      tw_av(2*id-1,iz)=tw_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*&
                                       tw(2*id-1,iz,nwr_u,ibui)**4
                      tw_av(2*id,iz)=tw_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*&
                                     tw(2*id,iz,nwr_u,ibui)**4
                      twlev_av(2*id-1,iz)=twlev_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*&
                                          twlev(2*id-1,iz,ibui)**4
                      twlev_av(2*id,iz)=twlev_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*&
                                        twlev(2*id,iz,ibui)**4
                      sfw_av(2*id-1,iz)=sfw_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*sfw(2*id-1,iz,ibui)
                      sfw_av(2*id,iz)=sfw_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*sfw(2*id,iz,ibui)
                      sfwind_av(2*id-1,iz)=sfwind_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*sfwin(2*id-1,iz,ibui)
                      sfwind_av(2*id,iz)=sfwind_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*sfwin(2*id,iz,ibui)
                   endif
                enddo
                tw_av(2*id-1,iz)=tw_av(2*id-1,iz)**(1./4.)
                tw_av(2*id,iz)=tw_av(2*id,iz)**(1./4.)
                twlev_av(2*id-1,iz)=twlev_av(2*id-1,iz)**(1./4.)
                twlev_av(2*id,iz)=twlev_av(2*id,iz)**(1./4.)
            endif
         enddo !iz         
      enddo !id
      return
      end subroutine averaging_temp
! ===6=8===============================================================72       
! ===6=8===============================================================72    

      subroutine modif_rad(iurb,nd,nz_u,z,ws,drst,strd,ss,pb,    &
                          tw,tg,twlev,albg,albw,emw,emg,pwin,albwin,   &
                          emwin,fww,fwg,fgw,fsw,fsg,             &
                          zr,deltar,ah,                          &    
                          rs,rl,rsw,rsg,rlw,rlg)                       
 
! ----------------------------------------------------------------------
! This routine computes the modification of the short wave and 
!  long wave radiation due to the buildings.
! ----------------------------------------------------------------------

      implicit none
 
 
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      integer iurb              ! current urban class
      integer nd                ! Number of street direction for the current urban class
      integer nz_u              ! Number of layer in the urban grid
      real z(nz_um)           ! Height of the urban grid levels
      real ws(ndm)              ! Street widths of the current urban class
      real drst(ndm)            ! street directions for the current urban class
      real strd(ndm)            ! Street lengths for the current urban class
      real ss(nz_um)          ! probability to have a building with height h
      real pb(nz_um)          ! probability to have a building with an height equal
      real tw(2*ndm,nz_um)    ! Temperature in each layer of the wall [K]
      real tg(ndm,ng_u)         ! Temperature in each layer of the ground [K]
      real albg                 ! Albedo of the ground for the current urban class
      real albw                 ! Albedo of the wall for the current urban class
      real emg                  ! Emissivity of ground for the current urban class
      real emw                  ! Emissivity of wall for the current urban class
      real fgw(nz_um,ndm,nurbm)       ! View factors from ground to wall
      real fsg(ndm,nurbm)             ! View factors from sky to ground
      real fsw(nz_um,ndm,nurbm)       ! View factors from sky to wall
      real fws(nz_um,ndm,nurbm)       ! View factors from wall to sky
      real fwg(nz_um,ndm,nurbm)       ! View factors from wall to ground
      real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
      real ah                   ! Hour angle (it should come from the radiation routine)
      real zr                   ! zenith angle
      real deltar               ! Declination of the sun
      real rs                   ! solar radiation
      real rl                   ! downward flux of the longwave radiation
!
!New variables BEM
!
      real twlev(2*ndm,nz_um)         ! Window temperature in BEM [K]
      real pwin                       ! Coverage area fraction of windows in the walls of the buildings 
      real albwin                     ! Albedo of the windows for the current urban class
      real emwin                      ! Emissivity of the windows for the current urban class
      real alb_av                     ! Averaged albedo (window and wall)
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
      real rlg(ndm)             ! Long wave radiation at the ground
      real rlw(2*ndm,nz_um)     ! Long wave radiation at the walls
      real rsg(ndm)             ! Short wave radiation at the ground
      real rsw(2*ndm,nz_um)     ! Short wave radiation at the walls

! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------

      integer id,iz

!  Calculation of the shadow effects

      call shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z,  &
                      rs,rsw,rsg)

! Calculation of the reflection effects          
      do id=1,nd
         call long_rad(iurb,nz_u,id,emw,emg,emwin,pwin,twlev,      &
                      fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)

         alb_av=pwin*albwin+(1.-pwin)*albw
         
         call short_rad(iurb,nz_u,id,alb_av,albg,fwg,fww,fgw,rsg,rsw,pb)
  
      enddo
      return
      end subroutine modif_rad


! ===6=8===============================================================72  
! ===6=8===============================================================72     

      subroutine surf_temp(nd,pr,dt,rl,rsg,rlg,  &
                           tg,alag,csg,emg,albg,ptg,sfg,gfg) 

! ----------------------------------------------------------------------
! Computation of the surface temperatures for walls, ground and roofs 
! ----------------------------------------------------------------------

      implicit none
                  
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      
      integer nd                ! Number of street direction for the current urban class
      real alag(ng_u)           ! Ground thermal diffusivity for the current urban class [m^2 s^-1] 

      real albg                 ! Albedo of the ground for the current urban class

      real csg(ng_u)            ! Specific heat of the ground material of the current urban class [J m^3 K^-1]

      real dt                   ! Time step
      real emg                  ! Emissivity of ground for the current urban class

      real pr(nz_um)            ! Air pressure
      real rl                   ! Downward flux of the longwave radiation
      real rlg(ndm)             ! Long wave radiation at the ground
      real rsg(ndm)             ! Short wave radiation at the ground
      real sfg(ndm)             ! Sensible heat flux from ground (road)

      real gfg(ndm)             ! Heat flux transferred from the surface of the ground (road) toward the interior

      real tg(ndm,ng_u)         ! Temperature in each layer of the ground [K]

      

! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
      real ptg(ndm)             ! Ground potential temperatures 


! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      integer id,ig,ir,iw,iz

      real rtg(ndm)             ! Total radiation at ground(road) surface (solar+incoming long+outgoing long)

      real tg_tmp(ng_u)


      real dzg_u(ng_u)          ! Layer sizes in the ground
      
      data dzg_u /0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------

        
   
      do id=1,nd

!      Calculation for the ground surfaces
       do ig=1,ng_u
        tg_tmp(ig)=tg(id,ig)
       end do
!	        
       call soil_temp(ng_u,dzg_u,tg_tmp,ptg(id),alag,csg,      &
                     rsg(id),rlg(id),pr(1),                    &
                     dt,emg,albg,                              &
                     rtg(id),sfg(id),gfg(id))    
       do ig=1,ng_u
        tg(id,ig)=tg_tmp(ig)
       end do
	
      end do !id
      
      return
      end subroutine surf_temp
     
! ===6=8===============================================================72     
! ===6=8===============================================================72  

      subroutine buildings(iurb,nd,nz,z0,ua_u,va_u,pt_u,pt0_u,       &
                        ptg,ptr,da_u,ptw,ptwin,pwin,                 &
                        drst,uva_u,vva_u,uvb_u,vvb_u,                &
                        tva_u,tvb_u,evb_u,qvb_u,qhb_u,               &
                        uhb_u,vhb_u,thb_u,ehb_u,ss,dt,sfw,sfg,sfr,   &
                        sfwin,pb,bs_u,dz_u,sflev,lflev,sfvlev,lfvlev,tvb_ac)                  

! ----------------------------------------------------------------------
! This routine computes the sources or sinks of the different quantities 
! on the urban grid. The actual calculation is done in the subroutines 
! called flux_wall, and flux_flat.
! ----------------------------------------------------------------------

      implicit none

        
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      integer nd                ! Number of street direction for the current urban class
      integer nz                ! number of vertical space steps
      real ua_u(nz_um)          ! Wind speed in the x direction on the urban grid
      real va_u(nz_um)          ! Wind speed in the y direction on the urban grid
      real da_u(nz_um)          ! air density on the urban grid
      real drst(ndm)            ! Street directions for the current urban class
      real dz
      real pt_u(nz_um)          ! Potential temperature on the urban grid
      real pt0_u(nz_um)         ! reference potential temperature on the urban grid
      real ptg(ndm)             ! Ground potential temperatures 
      real ptr(ndm,nz_um)       ! Roof potential temperatures 
      real ptw(2*ndm,nz_um,nbui_max)     ! Walls potential temperatures 
      real ss(nz_um)            ! probability to have a building with height h
      real pb(nz_um)
      real z0(ndm,nz_um)        ! Roughness lengths "profiles"
      real dt ! time step
      integer iurb              !Urban class

!
!New variables (BEM)
!
      real bs_u(ndm,nurbm)    ! Building width [m]
      real dz_u               ! Urban grid resolution
      real sflev(nz_um,nz_um)     ! sensible heat flux due to the air conditioning systems  [W]
      real lflev(nz_um,nz_um)     ! latent heat flux due to the air conditioning systems  [W]
      real sfvlev(nz_um,nz_um)    ! sensible heat flux due to ventilation [W]
      real lfvlev(nz_um,nz_um)    ! latent heat flux due to ventilation [W]
      real qvb_u(2*ndm,nz_um)
      real qhb_u(ndm,nz_um)
      real ptwin(2*ndm,nz_um,nbui_max)  ! window potential temperature
      real pwin
      real tvb_ac(2*ndm,nz_um)
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
! vertical surfaces (walls) and horizontal surfaces (roofs and street)
! The fluxes can be computed as follow: Fluxes of X = A*X + B
!  Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u

      real uhb_u(ndm,nz_um)   ! U (wind component) Horizontal surfaces, B (explicit) term
      real uva_u(2*ndm,nz_um)   ! U (wind component)   Vertical surfaces, A (implicit) term
      real uvb_u(2*ndm,nz_um)   ! U (wind component)   Vertical surfaces, B (explicit) term
      real vhb_u(ndm,nz_um)   ! V (wind component) Horizontal surfaces, B (explicit) term
      real vva_u(2*ndm,nz_um)   ! V (wind component)   Vertical surfaces, A (implicit) term
      real vvb_u(2*ndm,nz_um)   ! V (wind component)   Vertical surfaces, B (explicit) term
      real thb_u(ndm,nz_um)   ! Temperature        Horizontal surfaces, B (explicit) term
      real tva_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, A (implicit) term
      real tvb_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, B (explicit) term
      real ehb_u(ndm,nz_um)   ! Energy (TKE)       Horizontal surfaces, B (explicit) term
      real evb_u(2*ndm,nz_um)   ! Energy (TKE)         Vertical surfaces, B (explicit) term
      real sfw(2*ndm,nz_um,nbui_max)   ! sensible heat flux from walls
      real sfwin(2*ndm,nz_um,nbui_max) ! sensible heat flux form windows
      real sfr(ndm,nz_um)           ! sensible heat flux from roof
      real sfg(ndm)                 ! sensible heat flux from street

! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      real d_urb(nz_um)
      real uva_tmp
      real vva_tmp
      real uvb_tmp
      real vvb_tmp 
      real evb_tmp     
      integer nlev(nz_um)
      integer id,iz,ibui,nbui

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
      dz=dz_u
      ibui=0
      nbui=0
      nlev=0
      d_urb=0.
      
      uva_u=0.
      uvb_u=0.
      vhb_u=0.
      vva_u=0.
      vvb_u=0.
      thb_u=0.
      tva_u=0.
      tvb_u=0.
      tvb_ac=0.
      ehb_u=0.
      evb_u=0.
      qvb_u=0.
      qhb_u=0.

      do iz=1,nz_um		   
         if(ss(iz).gt.0) then		
           ibui=ibui+1
           d_urb(ibui)=ss(iz)
           nlev(ibui)=iz-1
           nbui=ibui		               
         endif
      enddo
      if (nbui.gt.nbui_max) then
         write(*,*) 'nbui_max must be increased to',nbui
         stop
      endif
      do id=1,nd

!        Calculation at the ground surfaces

         call flux_flat(dz,z0(id,1),ua_u(1),va_u(1),pt_u(1),pt0_u(1),  &
                       ptg(id),uhb_u(id,1),                            & 
                       vhb_u(id,1),sfg(id),ehb_u(id,1),da_u(1))          
         thb_u(id,1)=- sfg(id)/(da_u(1)*cp_u)  
              
!        Calculation at the roof surfaces    

         do iz=2,nz
            if(ss(iz).gt.0)then
               call flux_flat(dz,z0(id,iz),ua_u(iz),                  &              
                       va_u(iz),pt_u(iz),pt0_u(iz),                   &   
                       ptr(id,iz),uhb_u(id,iz),                       &   
                       vhb_u(id,iz),sfr(id,iz),ehb_u(id,iz),da_u(iz)) 
               thb_u(id,iz)=- sfr(id,iz)/(da_u(iz)*cp_u) 
            else
               uhb_u(id,iz) = 0.0
               vhb_u(id,iz) = 0.0
               thb_u(id,iz) = 0.0
               ehb_u(id,iz) = 0.0
            endif
         end do

!        Calculation at the wall surfaces        
 
         do ibui=1,nbui
         do iz=1,nlev(ibui)  
                   
            call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz),             &  
                        ptw(2*id-1,iz,ibui),ptwin(2*id-1,iz,ibui),          &   
                        uva_tmp,vva_tmp,                                    &   
                        uvb_tmp,vvb_tmp,                                    &   
                        sfw(2*id-1,iz,ibui),sfwin(2*id-1,iz,ibui),          &   
                        evb_tmp,drst(id),dt)      
   
            if (pb(iz+1).gt.0.) then

                    uva_u(2*id-1,iz)=uva_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*uva_tmp
                    vva_u(2*id-1,iz)=vva_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*vva_tmp
                    uvb_u(2*id-1,iz)=uvb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*uvb_tmp
                    vvb_u(2*id-1,iz)=vvb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*vvb_tmp
                    evb_u(2*id-1,iz)=evb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*evb_tmp
                    tvb_u(2*id-1,iz)=tvb_u(2*id-1,iz)-(d_urb(ibui)/pb(iz+1))*                       &
                                    (sfw(2*id-1,iz,ibui)*(1.-pwin)+sfwin(2*id-1,iz,ibui)*pwin)/     &
                                    da_u(iz)/cp_u-(1./4.)*(d_urb(ibui)/pb(iz+1))*(sfvlev(iz,ibui)-sflev(iz,ibui))/&
                                    (dz*bs_u(id,iurb))/da_u(iz)/cp_u
                    tvb_ac(2*id-1,iz)=tvb_ac(2*id-1,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(-sflev(iz,ibui))/&
                                    (dz*bs_u(id,iurb))/da_u(iz)/cp_u
                    qvb_u(2*id-1,iz)=qvb_u(2*id-1,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(lfvlev(iz,ibui)-lflev(iz,ibui))/&
                                    (dz*bs_u(id,iurb))/da_u(iz)/latent
                                     
            endif

            call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz),    &   
                        ptw(2*id,iz,ibui),ptwin(2*id,iz,ibui),     &    
                        uva_tmp,vva_tmp,                           &    
                        uvb_tmp,vvb_tmp,                           &    
                        sfw(2*id,iz,ibui),sfwin(2*id,iz,ibui),     &   
                        evb_tmp,drst(id),dt) 

            if (pb(iz+1).gt.0.) then

                    uva_u(2*id,iz)=uva_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*uva_tmp
                    vva_u(2*id,iz)=vva_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*vva_tmp
                    uvb_u(2*id,iz)=uvb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*uvb_tmp
                    vvb_u(2*id,iz)=vvb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*vvb_tmp
                    evb_u(2*id,iz)=evb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*evb_tmp
                    tvb_u(2*id,iz)=tvb_u(2*id,iz)-(d_urb(ibui)/pb(iz+1))*                    &
                                    (sfw(2*id,iz,ibui)*(1.-pwin)+sfwin(2*id,iz,ibui)*pwin)/  &
                                     da_u(iz)/cp_u-(1./4.)*(d_urb(ibui)/pb(iz+1))*(sfvlev(iz,ibui)-sflev(iz,ibui))/&
                                    (dz*bs_u(id,iurb))/da_u(iz)/cp_u
                    tvb_ac(2*id,iz)=tvb_ac(2*id,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(-sflev(iz,ibui))/&
                                    (dz*bs_u(id,iurb))/da_u(iz)/cp_u
                    qvb_u(2*id,iz)=qvb_u(2*id,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(lfvlev(iz,ibui)-lflev(iz,ibui))/&
                                    (dz*bs_u(id,iurb))/da_u(iz)/latent

            endif
!
          enddo !iz
         enddo !ibui
         
      end do !id
                
      return
      end subroutine buildings
      

! ===6=8===============================================================72
! ===6=8===============================================================72

        subroutine urban_meso(nd,kms,kme,kts,kte,nz_u,z,dz,z_u,pb,ss,bs,ws,sf,vl,    &
                             uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &       
                             uhb_u,vhb_u,thb_u,ehb_u,qhb_u,qvb_u,       &      
                             a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e,b_q,tvb_ac,b_ac)           

! ----------------------------------------------------------------------
!  This routine interpolates the parameters from the "urban grid" to the
!  "mesoscale grid".
!  See p300-301 Appendix B.2 of the BLM paper.  
! ----------------------------------------------------------------------

      implicit none

! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
! Data relative to the "mesoscale grid"
      integer kms,kme,kts,kte               
      real z(kms:kme)              ! Altitude above the ground of the cell interface
      real dz(kms:kme)               ! Vertical space steps

! Data relative to the "uban grid"
      integer nz_u              ! Number of layer in the urban grid
      integer nd                ! Number of street direction for the current urban class
      real bs(ndm)              ! Building widths of the current urban class
      real ws(ndm)              ! Street widths of the current urban class
      real z_u(nz_um)         ! Height of the urban grid levels
      real pb(nz_um)          ! Probability to have a building with an height equal
      real ss(nz_um)          ! Probability to have a building with height h
      real uhb_u(ndm,nz_um)   ! U (x-wind component) Horizontal surfaces, B (explicit) term
      real uva_u(2*ndm,nz_um)   ! U (x-wind component) Vertical surfaces, A (implicit) term
      real uvb_u(2*ndm,nz_um)   ! U (x-wind component) Vertical surfaces, B (explicit) term
      real vhb_u(ndm,nz_um)   ! V (y-wind component) Horizontal surfaces, B (explicit) term
      real vva_u(2*ndm,nz_um)   ! V (y-wind component) Vertical surfaces, A (implicit) term
      real vvb_u(2*ndm,nz_um)   ! V (y-wind component) Vertical surfaces, B (explicit) term
      real thb_u(ndm,nz_um)   ! Temperature        Horizontal surfaces, B (explicit) term
      real tva_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, A (implicit) term
      real tvb_u(2*ndm,nz_um)   ! Temperature          Vertical surfaces, B (explicit) term
      real tvb_ac(2*ndm,nz_um)
      real ehb_u(ndm,nz_um)   ! Energy (TKE)       Horizontal surfaces, B (explicit) term
      real evb_u(2*ndm,nz_um)   ! Energy (TKE)         Vertical surfaces, B (explicit) term
!
!New variables for BEM
!
      real qhb_u(ndm,nz_um)
      real qvb_u(2*ndm,nz_um)
     

! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
! Data relative to the "mesoscale grid"
      real sf(kms:kme)             ! Surface of the "mesoscale grid" cells taking into account the buildings
      real vl(kms:kme)               ! Volume of the "mesoscale grid" cells taking into account the buildings
      real a_u(kms:kme)              ! Implicit component of the momentum sources or sinks in the X-direction
      real a_v(kms:kme)              ! Implicit component of the momentum sources or sinks in the Y-direction
      real a_t(kms:kme)              ! Implicit component of the heat sources or sinks
      real a_e(kms:kme)              ! Implicit component of the TKE sources or sinks
      real b_u(kms:kme)              ! Explicit component of the momentum sources or sinks in the X-direction
      real b_v(kms:kme)              ! Explicit component of the momentum sources or sinks in the Y-direction
      real b_t(kms:kme)              ! Explicit component of the heat sources or sinks
      real b_ac(kms:kme)
      real b_e(kms:kme)              ! Explicit component of the TKE sources or sinks
      real b_q(kms:kme)               ! Explicit component of the humidity sources or sinks
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      real dzz
      real fact
      integer id,iz,iz_u
      real se,sr,st,su,sv,sq
      real uet(kms:kme)                ! Contribution to TKE due to walls
      real veb,vta,vtb,vte,vtot,vua,vub,vva,vvb,vqb,vtb_ac


! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ---------------------------------------------------------------------- 

! initialisation

      do iz=kts,kte
         a_u(iz)=0.
         a_v(iz)=0.
         a_t(iz)=0.
         a_e(iz)=0.
         b_u(iz)=0.
         b_v(iz)=0.
         b_e(iz)=0.
         b_t(iz)=0.
         b_ac(iz)=0.
         uet(iz)=0.
         b_q(iz)=0.
      end do
            
! horizontal surfaces
      do iz=kts,kte
         sf(iz)=0.
         vl(iz)=0.
      enddo
      sf(kte+1)=0. 
      
      do id=1,nd      
         do iz=kts+1,kte+1
            sr=0.
            do iz_u=2,nz_u
               if(z(iz).lt.z_u(iz_u).and.z(iz).ge.z_u(iz_u-1))then
                  sr=pb(iz_u)
               endif
            enddo
            sf(iz)=sf(iz)+((ws(id)+(1.-sr)*bs(id))/(ws(id)+bs(id)))/nd
         enddo
      enddo

! volume      
      do iz=kts,kte
         do id=1,nd
            vtot=0.
            do iz_u=1,nz_u
               dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
               vtot=vtot+pb(iz_u+1)*dzz
            enddo
            vtot=vtot/(z(iz+1)-z(iz))
            vl(iz)=vl(iz)+(1.-vtot*bs(id)/(ws(id)+bs(id)))/nd
         enddo
      enddo
      
! horizontal surface impact  

      do id=1,nd
      
         fact=1./vl(kts)/dz(kts)*ws(id)/(ws(id)+bs(id))/nd
         b_t(kts)=b_t(kts)+thb_u(id,1)*fact
         b_u(kts)=b_u(kts)+uhb_u(id,1)*fact
         b_v(kts)=b_v(kts)+vhb_u(id,1)*fact 
         b_e(kts)=b_e(kts)+ehb_u(id,1)*fact*(z_u(2)-z_u(1))
         b_q(kts)=b_q(kts)+qhb_u(id,1)*fact         

         do iz=kts,kte
            st=0.
            su=0.
            sv=0.
            se=0.
            sq=0.
            do iz_u=2,nz_u
               if(z(iz).le.z_u(iz_u).and.z(iz+1).gt.z_u(iz_u))then
                  st=st+ss(iz_u)*thb_u(id,iz_u)
                  su=su+ss(iz_u)*uhb_u(id,iz_u)
                  sv=sv+ss(iz_u)*vhb_u(id,iz_u)          
                  se=se+ss(iz_u)*ehb_u(id,iz_u)*(z_u(iz_u+1)-z_u(iz_u))
                  sq=sq+ss(iz_u)*qhb_u(id,iz_u)
               endif
            enddo
      
            fact=bs(id)/(ws(id)+bs(id))/vl(iz)/dz(iz)/nd
            b_t(iz)=b_t(iz)+st*fact
            b_u(iz)=b_u(iz)+su*fact
            b_v(iz)=b_v(iz)+sv*fact
            b_e(iz)=b_e(iz)+se*fact
            b_q(iz)=b_q(iz)+sq*fact
         enddo
      enddo              

! vertical surface impact
           
      do iz=kts,kte 
         uet(iz)=0.
         do id=1,nd              
            vtb=0.
            vtb_ac=0.
            vta=0.
            vua=0.
            vub=0.
            vva=0.
            vvb=0.
            veb=0.
	    vte=0.
            vqb=0.
            do iz_u=1,nz_u
               dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
               fact=dzz/(ws(id)+bs(id))
               vtb=vtb+pb(iz_u+1)*                                  &        
                    (tvb_u(2*id-1,iz_u)+tvb_u(2*id,iz_u))*fact 
               vtb_ac=vtb_ac+pb(iz_u+1)*                            &        
                    (tvb_ac(2*id-1,iz_u)+tvb_ac(2*id,iz_u))*fact     
               vta=vta+pb(iz_u+1)*                                  &        
                   (tva_u(2*id-1,iz_u)+tva_u(2*id,iz_u))*fact
               vua=vua+pb(iz_u+1)*                                  &        
                    (uva_u(2*id-1,iz_u)+uva_u(2*id,iz_u))*fact
               vva=vva+pb(iz_u+1)*                                  &        
                    (vva_u(2*id-1,iz_u)+vva_u(2*id,iz_u))*fact
               vub=vub+pb(iz_u+1)*                                  &        
                    (uvb_u(2*id-1,iz_u)+uvb_u(2*id,iz_u))*fact
               vvb=vvb+pb(iz_u+1)*                                  &        
                    (vvb_u(2*id-1,iz_u)+vvb_u(2*id,iz_u))*fact
               veb=veb+pb(iz_u+1)*                                  &        
                    (evb_u(2*id-1,iz_u)+evb_u(2*id,iz_u))*fact
               vqb=vqb+pb(iz_u+1)*                                  &        
                    (qvb_u(2*id-1,iz_u)+qvb_u(2*id,iz_u))*fact   
            enddo
           
            fact=1./vl(iz)/dz(iz)/nd
            b_t(iz)=b_t(iz)+vtb*fact
            b_ac(iz)=b_ac(iz)+vtb_ac*fact
            a_t(iz)=a_t(iz)+vta*fact
            a_u(iz)=a_u(iz)+vua*fact
            a_v(iz)=a_v(iz)+vva*fact
            b_u(iz)=b_u(iz)+vub*fact
            b_v(iz)=b_v(iz)+vvb*fact
            b_e(iz)=b_e(iz)+veb*fact
            uet(iz)=uet(iz)+vte*fact
            b_q(iz)=b_q(iz)+vqb*fact
         enddo              
      enddo
      

      
      return
      end subroutine urban_meso

! ===6=8===============================================================72 
! ===6=8===============================================================72 

      subroutine interp_length(nd,kms,kme,kts,kte,nz_u,z_u,z,ss,ws,bs, &
                             dlg,dl_u)

! ----------------------------------------------------------------------     
!    Calculation of the length scales
!    See p272-274 formula (22) and (24) of the BLM paper    
! ----------------------------------------------------------------------     
     
      implicit none


! ----------------------------------------------------------------------     
! INPUT:
! ----------------------------------------------------------------------     
      integer kms,kme,kts,kte                
      real z(kms:kme)              ! Altitude above the ground of the cell interface
      integer nd                ! Number of street direction for the current urban class
      integer nz_u              ! Number of levels in the "urban grid"
      real z_u(nz_um)         ! Height of the urban grid levels
      real bs(ndm)              ! Building widths of the current urban class
      real ss(nz_um)          ! Probability to have a building with height h
      real ws(ndm)              ! Street widths of the current urban class


! ----------------------------------------------------------------------     
! OUTPUT:
! ----------------------------------------------------------------------     
      real dlg(kms:kme)              ! Height above ground (L_ground in formula (24) of the BLM paper). 
      real dl_u(kms:kme)             ! Length scale (lb in formula (22) ofthe BLM paper).

! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      real dlgtmp
      integer id,iz,iz_u
      real sftot
      real ulu,ssl

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
   
        do iz=kts,kte
         ulu=0.
         ssl=0.
         do id=1,nd        
          do iz_u=2,nz_u
           if(z_u(iz_u).gt.z(iz))then
            ulu=ulu+ss(iz_u)/z_u(iz_u)/nd
            ssl=ssl+ss(iz_u)/nd
           endif
          enddo
         enddo

        if(ulu.ne.0)then
          dl_u(iz)=ssl/ulu
         else
          dl_u(iz)=0.
         endif
        enddo
       

        do iz=kts,kte
         dlg(iz)=0.
          do id=1,nd
           sftot=ws(id)  
           dlgtmp=ws(id)/((z(iz)+z(iz+1))/2.)
           do iz_u=1,nz_u
            if((z(iz)+z(iz+1))/2..gt.z_u(iz_u))then
             dlgtmp=dlgtmp+ss(iz_u)*bs(id)/                           &                
                    ((z(iz)+z(iz+1))/2.-z_u(iz_u))
             sftot=sftot+ss(iz_u)*bs(id)
            endif
           enddo
           dlg(iz)=dlg(iz)+dlgtmp/sftot/nd
         enddo
         dlg(iz)=1./dlg(iz)
        enddo
        
       return
       end subroutine interp_length

! ===6=8===============================================================72
! ===6=8===============================================================72   

      subroutine shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z,    &
                           rs,rsw,rsg)
        
! ----------------------------------------------------------------------
!         Modification of short wave radiation to take into account
!         the shadow produced by the buildings
! ----------------------------------------------------------------------

      implicit none
     
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      integer nd                ! Number of street direction for the current urban class
      integer nz_u              ! number of vertical layers defined in the urban grid
      real ah                   ! Hour angle (it should come from the radiation routine)
      real deltar               ! Declination of the sun
      real drst(ndm)            ! street directions for the current urban class
      real rs                   ! solar radiation
      real ss(nz_um)          ! probability to have a building with height h
      real pb(nz_um)          ! Probability that a building has an height greater or equal to h
      real ws(ndm)              ! Street width of the current urban class
      real z(nz_um)           ! Height of the urban grid levels
      real zr                   ! zenith angle

! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
      real rsg(ndm)             ! Short wave radiation at the ground
      real rsw(2*ndm,nz_um)     ! Short wave radiation at the walls

! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      integer id,iz,jz
      real aae,aaw,bbb,phix,rd,rtot,wsd
      
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------

      if(rs.eq.0.or.sin(zr).eq.1)then
         do id=1,nd
            rsg(id)=0.
            do iz=1,nz_u
               rsw(2*id-1,iz)=0.
               rsw(2*id,iz)=0.
            enddo
         enddo
      else            
!test              
         if(abs(sin(zr)).gt.1.e-10)then
          if(cos(deltar)*sin(ah)/sin(zr).ge.1)then
           bbb=pi/2.
          elseif(cos(deltar)*sin(ah)/sin(zr).le.-1)then
           bbb=-pi/2.
          else
           bbb=asin(cos(deltar)*sin(ah)/sin(zr))
          endif
         else
          if(cos(deltar)*sin(ah).ge.0)then
           bbb=pi/2.
          elseif(cos(deltar)*sin(ah).lt.0)then
           bbb=-pi/2.
          endif
         endif

         phix=zr
           
         do id=1,nd
        
            rsg(id)=0.
           
            aae=bbb-drst(id)
            aaw=bbb-drst(id)+pi
                    
            do iz=1,nz_u
               rsw(2*id-1,iz)=0.
               rsw(2*id,iz)=0.          
            if (pb(iz+1).gt.0.) then          
               do jz=1,nz_u                    
                if(abs(sin(aae)).gt.1.e-10)then
                  call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aae,   &    
                      ws(id),rd)                  
                 rsw(2*id-1,iz)=rsw(2*id-1,iz)+rs*rd*ss(jz+1)/pb(iz+1)
                endif
              
                if(abs(sin(aaw)).gt.1.e-10)then
                  call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aaw,   &    
                      ws(id),rd)
                 rsw(2*id,iz)=rsw(2*id,iz)+rs*rd*ss(jz+1)/pb(iz+1)                  
                endif
               enddo             
            endif  
            enddo
        if(abs(sin(aae)).gt.1.e-10)then
            wsd=abs(ws(id)/sin(aae))
              
            do jz=1,nz_u           
               rd=max(0.,wsd-z(jz+1)*tan(phix))
               rsg(id)=rsg(id)+rs*rd*ss(jz+1)/wsd          
            enddo
            rtot=0.
           
            do iz=1,nz_u
               rtot=rtot+(rsw(2*id,iz)+rsw(2*id-1,iz))*            &
                         (z(iz+1)-z(iz))
            enddo
            rtot=rtot+rsg(id)*ws(id)
        else
            rsg(id)=rs
        endif
            
         enddo
      endif
         
      return
      end subroutine shadow_mas
         
! ===6=8===============================================================72     
! ===6=8===============================================================72     

      subroutine shade_wall(z1,z2,hu,phix,aa,ws,rd)

! ----------------------------------------------------------------------
! This routine computes the effects of a shadow induced by a building of 
! height hu, on a portion of wall between z1 and z2. See equation A10, 
! and correction described below formula A11, and figure A1. Basically rd
! is the ratio between the horizontal surface illuminated and the portion
! of wall. Referring to figure A1, multiplying radiation flux density on 
! a horizontal surface (rs) by x1-x2 we have the radiation energy per 
! unit time. Dividing this by z2-z1, we obtain the radiation flux 
! density reaching the portion of the wall between z2 and z1 
! (everything is assumed in 2D)
! ----------------------------------------------------------------------

      implicit none
      
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      real aa                   ! Angle between the sun direction and the face of the wall (A12)
      real hu                   ! Height of the building that generates the shadow
      real phix                 ! Solar zenith angle
      real ws                   ! Width of the street
      real z1                   ! Height of the level z(iz)
      real z2                   ! Height of the level z(iz+1)

! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
      real rd                   ! Ratio between (x1-x2)/(z2-z1), see Fig. 1A. 
                                ! Multiplying rd by rs (radiation flux 
                                ! density on a horizontal surface) gives 
                                ! the radiation flux density on the 
                                ! portion of wall between z1 and z2. 
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      real x1,x2                ! x1,x2 see Fig. A1.

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------

      x1=min((hu-z1)*tan(phix),max(0.,ws/sin(aa)))
      
      x2=max((hu-z2)*tan(phix),0.)

      rd=max(0.,sin(aa)*(max(0.,x1-x2))/(z2-z1))
      
      return
      end subroutine shade_wall

! ===6=8===============================================================72     
! ===6=8===============================================================72     

      subroutine long_rad(iurb,nz_u,id,emw,emg,emwin,pwin,twlev,&
                         fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)

! ----------------------------------------------------------------------
! This routine computes the effects of the reflections of long-wave 
! radiation in the street canyon by solving the system 
! of 2*nz_u+1 eqn. in 2*nz_u+1
! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
! The system is solved by solving A X= B,
! with A matrix B vector, and X solution. 
! ----------------------------------------------------------------------

      implicit none

  
      
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      real emg                        ! Emissivity of ground for the current urban class
      real emw                        ! Emissivity of wall for the current urban class
      real fgw(nz_um,ndm,nurbm)       ! View factors from ground to wall
      real fsg(ndm,nurbm)             ! View factors from sky to ground
      real fsw(nz_um,ndm,nurbm)       ! View factors from sky to wall
      real fwg(nz_um,ndm,nurbm)       ! View factors from wall to ground
      real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
      integer id                      ! Current street direction
      integer iurb                    ! Current urban class
      integer nz_u                    ! Number of layer in the urban grid
      real pb(nz_um)                  ! Probability to have a building with an height equal
      real rl                         ! Downward flux of the longwave radiation
      real tg(ndm,ng_u)               ! Temperature in each layer of the ground [K]
      real tw(2*ndm,nz_um)            ! Temperature in each layer of the wall [K]
!
!New Variables for BEM
!
      real twlev(2*ndm,nz_um)         ! Window temperature in BEM [K]
      real emwin                      ! Emissivity of windows
      real pwin                       ! Coverage area fraction of windows in the walls of the buildings (BEM)
      

! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
      real rlg(ndm)                   ! Long wave radiation at the ground
      real rlw(2*ndm,nz_um)           ! Long wave radiation at the walls

! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      integer i,j
      real aaa(2*nz_um+1,2*nz_um+1)   ! terms of the matrix
      real bbb(2*nz_um+1)             ! terms of the vector

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------


! west wall
       
      do i=1,nz_u
        
        do j=1,nz_u
         aaa(i,j)=0.
        enddo
        
        aaa(i,i)=1.        
       
        do j=nz_u+1,2*nz_u
         aaa(i,j)=-(1.-emw*(1.-pwin)-emwin*pwin)* &
                  fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
        enddo
        
        aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i,id,iurb)
        
        bbb(i)=fsw(i,id,iurb)*rl+emg*fgw(i,id,iurb)*sigma*tg(id,ng_u)**4
        do j=1,nz_u
           bbb(i)=bbb(i)+pb(j+1)*sigma*fww(j,i,id,iurb)* &
                 (emw*(1.-pwin)*tw(2*id,j)**4+emwin*pwin*twlev(2*id,j)**4)+ &
                 fww(j,i,id,iurb)*rl*(1.-pb(j+1))
        enddo
        
       enddo
       
! east wall

       do i=1+nz_u,2*nz_u
        
        do j=1,nz_u
         aaa(i,j)=-(1.-emw*(1.-pwin)-emwin*pwin)*fww(j,i-nz_u,id,iurb)*pb(j+1)
        enddo
        
        do j=1+nz_u,2*nz_u
         aaa(i,j)=0.
        enddo
        
        aaa(i,i)=1.
        
        aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i-nz_u,id,iurb)
        
        bbb(i)=fsw(i-nz_u,id,iurb)*rl+  &     
               emg*fgw(i-nz_u,id,iurb)*sigma*tg(id,ng_u)**4

        do j=1,nz_u
         bbb(i)=bbb(i)+pb(j+1)*sigma*fww(j,i-nz_u,id,iurb)*  &   
                (emw*(1.-pwin)*tw(2*id-1,j)**4+emwin*pwin*twlev(2*id-1,j)**4)+&   
                fww(j,i-nz_u,id,iurb)*rl*(1.-pb(j+1))
        enddo
       
       enddo

! ground
       do j=1,nz_u
        aaa(2*nz_u+1,j)=-(1.-emw*(1.-pwin)-emwin*pwin)* &
                         fwg(j,id,iurb)*pb(j+1)
       enddo
       
       do j=nz_u+1,2*nz_u
        aaa(2*nz_u+1,j)=-(1.-emw*(1.-pwin)-emwin*pwin)* &
                         fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
       enddo
       
       aaa(2*nz_u+1,2*nz_u+1)=1.
       
       bbb(2*nz_u+1)=fsg(id,iurb)*rl
       
       do i=1,nz_u
        bbb(2*nz_u+1)=bbb(2*nz_u+1)+sigma*fwg(i,id,iurb)*pb(i+1)*         &
                      (emw*(1.-pwin)*(tw(2*id-1,i)**4+tw(2*id,i)**4)+     &
                      emwin*pwin*(twlev(2*id-1,i)**4+twlev(2*id,i)**4))+  &
                      2.*fwg(i,id,iurb)*(1.-pb(i+1))*rl                  
       enddo
   

     
       call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)

       do i=1,nz_u
        rlw(2*id-1,i)=bbb(i)
       enddo
       
       do i=nz_u+1,2*nz_u
        rlw(2*id,i-nz_u)=bbb(i)
       enddo
       
       rlg(id)=bbb(2*nz_u+1)
  
       return
       end subroutine long_rad
             
! ===6=8===============================================================72
! ===6=8===============================================================72

       subroutine short_rad(iurb,nz_u,id,albw,                        & 
                           albg,fwg,fww,fgw,rsg,rsw,pb)

! ----------------------------------------------------------------------
! This routine computes the effects of the reflections of short-wave 
! (solar) radiation in the street canyon by solving the system 
! of 2*nz_u+1 eqn. in 2*nz_u+1
! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
! The system is solved by solving A X= B,
! with A matrix B vector, and X solution. 
! ----------------------------------------------------------------------

      implicit none

  
      
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      real albg                 ! Albedo of the ground for the current urban class
      real albw                 ! Albedo of the wall for the current urban class
      real fgw(nz_um,ndm,nurbm)       ! View factors from ground to wall
      real fwg(nz_um,ndm,nurbm)       ! View factors from wall to ground
      real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
      integer id                ! current street direction 
      integer iurb              ! current urban class
      integer nz_u              ! Number of layer in the urban grid
      real pb(nz_um)          ! probability to have a building with an height equal

! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
      real rsg(ndm)             ! Short wave radiation at the ground
      real rsw(2*ndm,nz_um)     ! Short wave radiation at the walls

! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      integer i,j
      real aaa(2*nz_um+1,2*nz_um+1)  ! terms of the matrix
      real bbb(2*nz_um+1)            ! terms of the vector

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------

      
! west wall
       
      do i=1,nz_u
         do j=1,nz_u
            aaa(i,j)=0.
         enddo
         
         aaa(i,i)=1.        
         
         do j=nz_u+1,2*nz_u
            aaa(i,j)=-albw*fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
         enddo
         
         aaa(i,2*nz_u+1)=-albg*fgw(i,id,iurb)
         bbb(i)=rsw(2*id-1,i)
         
      enddo
       
! east wall

      do i=1+nz_u,2*nz_u
         do j=1,nz_u
            aaa(i,j)=-albw*fww(j,i-nz_u,id,iurb)*pb(j+1)
         enddo
         
         do j=1+nz_u,2*nz_u
            aaa(i,j)=0.
         enddo
         
        aaa(i,i)=1.
        aaa(i,2*nz_u+1)=-albg*fgw(i-nz_u,id,iurb)
        bbb(i)=rsw(2*id,i-nz_u)
        
      enddo

! ground

      do j=1,nz_u
         aaa(2*nz_u+1,j)=-albw*fwg(j,id,iurb)*pb(j+1)
      enddo
       
      do j=nz_u+1,2*nz_u
         aaa(2*nz_u+1,j)=-albw*fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
      enddo
       
      aaa(2*nz_u+1,2*nz_u+1)=1.
      bbb(2*nz_u+1)=rsg(id)
       
      call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)

      do i=1,nz_u
         rsw(2*id-1,i)=bbb(i)
      enddo
       
      do i=nz_u+1,2*nz_u
         rsw(2*id,i-nz_u)=bbb(i) 
      enddo
       
      rsg(id)=bbb(2*nz_u+1)
       
      return
      end subroutine short_rad
             

! ===6=8===============================================================72     
! ===6=8===============================================================72     
      
      subroutine gaussj(a,n,b,np)

! ----------------------------------------------------------------------
! This routine solve a linear system of n equations of the form
!              A X = B
!  where  A is a matrix a(i,j)
!         B a vector and X the solution
! In output b is replaced by the solution     
! ----------------------------------------------------------------------

      implicit none

! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      integer np
      real a(np,np)

! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
      real b(np)

! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      integer nmax
      parameter (nmax=150)

      real big,dum
      integer i,icol,irow
      integer j,k,l,ll,n
      integer ipiv(nmax)
      real pivinv

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
       
      do j=1,n
         ipiv(j)=0.
      enddo
       
      do i=1,n
         big=0.
         do j=1,n
            if(ipiv(j).ne.1)then
               do k=1,n
                  if(ipiv(k).eq.0)then
                     if(abs(a(j,k)).ge.big)then
                        big=abs(a(j,k))
                        irow=j
                        icol=k
                     endif
                  elseif(ipiv(k).gt.1)then
                     FATAL_ERROR('singular matrix in gaussj')
                  endif
               enddo
            endif
         enddo
         
         ipiv(icol)=ipiv(icol)+1
         
         if(irow.ne.icol)then
            do l=1,n
               dum=a(irow,l)
               a(irow,l)=a(icol,l)
               a(icol,l)=dum
            enddo
            
            dum=b(irow)
            b(irow)=b(icol)
            b(icol)=dum
          
         endif
         
         if(a(icol,icol).eq.0) FATAL_ERROR('singular matrix in gaussj')
         
         pivinv=1./a(icol,icol)
         a(icol,icol)=1
         
         do l=1,n
            a(icol,l)=a(icol,l)*pivinv
         enddo
         
         b(icol)=b(icol)*pivinv
         
         do ll=1,n
            if(ll.ne.icol)then
               dum=a(ll,icol)
               a(ll,icol)=0.
               do l=1,n
                  a(ll,l)=a(ll,l)-a(icol,l)*dum
               enddo
               
               b(ll)=b(ll)-b(icol)*dum
               
            endif
         enddo
      enddo
      
      return
      end subroutine gaussj
         
! ===6=8===============================================================72     
! ===6=8===============================================================72     
       
      subroutine soil_temp(nz,dz,temp,pt,ala,cs,                       &
                          rs,rl,press,dt,em,alb,rt,sf,gf)

! ----------------------------------------------------------------------
! This routine solves the Fourier diffusion equation for heat in 
! the material (wall, roof, or ground). Resolution is done implicitely.
! Boundary conditions are: 
! - fixed temperature at the interior
! - energy budget at the surface
! ----------------------------------------------------------------------

      implicit none

     
                
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      integer nz                ! Number of layers
      real ala(nz)              ! Thermal diffusivity in each layers [m^2 s^-1] 
      real alb                  ! Albedo of the surface
      real cs(nz)               ! Specific heat of the material [J m^3 K^-1]
      real dt                   ! Time step
      real em                   ! Emissivity of the surface
      real press                ! Pressure at ground level
      real rl                   ! Downward flux of the longwave radiation
      real rs                   ! Solar radiation
      real sf                   ! Sensible heat flux at the surface
      real temp(nz)             ! Temperature in each layer [K]
      real dz(nz)               ! Layer sizes [m]


! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
      real gf                   ! Heat flux transferred from the surface toward the interior
      real pt                   ! Potential temperature at the surface
      real rt                   ! Total radiation at the surface (solar+incoming long+outgoing long)

! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      integer iz
      real a(nz,3)
      real alpha
      real c(nz)
      real cddz(nz+2)
      real tsig

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
       
      tsig=temp(nz)
      alpha=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf
! Compute cddz=2*cd/dz  
        
      cddz(1)=ala(1)/dz(1)
      do iz=2,nz
         cddz(iz)=2.*ala(iz)/(dz(iz)+dz(iz-1))
      enddo
       
      a(1,1)=0.
      a(1,2)=1.
      a(1,3)=0.
      c(1)=temp(1)
          
      do iz=2,nz-1
         a(iz,1)=-cddz(iz)*dt/dz(iz)
         a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dz(iz)          
         a(iz,3)=-cddz(iz+1)*dt/dz(iz)
         c(iz)=temp(iz)
      enddo          
                     
      a(nz,1)=-dt*cddz(nz)/dz(nz)
      a(nz,2)=1.+dt*cddz(nz)/dz(nz)
      a(nz,3)=0.
      c(nz)=temp(nz)+dt*alpha/cs(nz)/dz(nz)

      
      call invert(nz,a,c,temp)

           
      pt=temp(nz)*(press/1.e+5)**(-rcp_u)

      rt=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)
                        
       gf=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf                                   
      return
      end subroutine soil_temp

! ===6=8===============================================================72 
! ===6=8===============================================================72 

      subroutine invert(n,a,c,x)

! ----------------------------------------------------------------------
!        Inversion and resolution of a tridiagonal matrix
!                   A X = C
! ----------------------------------------------------------------------

      implicit none
                
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
       integer n
       real a(n,3)              !  a(*,1) lower diagonal (Ai,i-1)
                                !  a(*,2) principal diagonal (Ai,i)
                                !  a(*,3) upper diagonal (Ai,i+1)
       real c(n)

! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
       real x(n)    

! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
       integer i

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
                     
       do i=n-1,1,-1                 
          c(i)=c(i)-a(i,3)*c(i+1)/a(i+1,2)
          a(i,2)=a(i,2)-a(i,3)*a(i+1,1)/a(i+1,2)
       enddo
       
       do i=2,n        
          c(i)=c(i)-a(i,1)*c(i-1)/a(i-1,2)
       enddo
        
       do i=1,n
          x(i)=c(i)/a(i,2)
       enddo

       return
       end subroutine invert
  

! ===6=8===============================================================72  
! ===6=8===============================================================72
  
      subroutine flux_wall(ua,va,pt,da,ptw,ptwin,uva,vva,uvb,vvb,  &
                           sfw,sfwin,evb,drst,dt)         
       
! ----------------------------------------------------------------------
! This routine computes the surface sources or sinks of momentum, tke,
! and heat from vertical surfaces (walls).   
! ----------------------------------------------------------------------
      implicit none   
         
! INPUT:
! -----
      real drst                 ! street directions for the current urban class
      real da                   ! air density
      real pt                   ! potential temperature
      real ptw                  ! Walls potential temperatures 
      real ptwin                ! Windows potential temperatures
      real ua                   ! wind speed
      real va                   ! wind speed
      real dt                   !time step

! OUTPUT:
! ------
! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
! vertical surfaces (walls).
! The fluxes can be computed as follow: Fluxes of X = A*X + B
! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u

      real uva                  ! U (wind component)   Vertical surfaces, A (implicit) term
      real uvb                  ! U (wind component)   Vertical surfaces, B (explicit) term
      real vva                  ! V (wind component)   Vertical surfaces, A (implicit) term
      real vvb                  ! V (wind component)   Vertical surfaces, B (explicit) term
!     real tva                  ! Temperature          Vertical surfaces, A (implicit) term
!     real tvb                  ! Temperature          Vertical surfaces, B (explicit) term
      real evb                  ! Energy (TKE)         Vertical surfaces, B (explicit) term
      real sfw                  ! Surfaces fluxes from the walls
      real sfwin                ! Surfaces fluxes from the windows

! LOCAL:
! -----
      real hc
      real hcwin
      real u_ort
      real vett


! -------------------------
! END VARIABLES DEFINITIONS
! -------------------------

      vett=(ua**2+va**2)**.5         
         
      u_ort=abs((cos(drst)*ua-sin(drst)*va))
       
      uva=-cdrag*u_ort/2.*cos(drst)*cos(drst)
      vva=-cdrag*u_ort/2.*sin(drst)*sin(drst)
         
      uvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*va
      vvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*ua         

      if (vett.lt.4.88) then   
         hc=5.678*(1.09+0.23*(vett/0.3048))  
      else
         hc=5.678*0.53*((vett/0.3048)**0.78)
      endif 

      if (hc.gt.da*cp_u/dt)then
         hc=da*cp_u/dt
      endif

       if (vett.lt.4.88) then
          hcwin=5.678*(0.99+0.21*(vett/0.3048))
       else
          hcwin=5.678*0.50*((vett/0.3048)**0.78)
       endif

       if (hcwin.gt.da*cp_u/dt) then
           hcwin=da*cp_u/dt
       endif
         
!         tvb=hc*ptw/da/cp_u
!         tva=-hc/da/cp_u
!!!!!!!!!!!!!!!!!!!!
! explicit 

      sfw=hc*(pt-ptw)
      sfwin=hcwin*(pt-ptwin)  
       
         
      evb=cdrag*(abs(u_ort)**3.)/2.
              
      return
      end subroutine flux_wall
         
! ===6=8===============================================================72
! ===6=8===============================================================72

      subroutine flux_flat(dz,z0,ua,va,pt,pt0,ptg,                     &
                          uhb,vhb,sf,ehb,da)
                                
! ----------------------------------------------------------------------
!           Calculation of the flux at the ground 
!           Formulation of Louis (Louis, 1979)       
! ----------------------------------------------------------------------

      implicit none

      real dz                   ! first vertical level
      real pt                   ! potential temperature
      real pt0                  ! reference potential temperature
      real ptg                  ! ground potential temperature
      real ua                   ! wind speed
      real va                   ! wind speed
      real z0                   ! Roughness length
      real da                   ! air density
      
     

! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
! Explicit component of the momentum, temperature and TKE sources or sinks on horizontal 
!  surfaces (roofs and street)
! The fluxes can be computed as follow: Fluxes of X = B
!  Example: Momentum fluxes on horizontal surfaces =  uhb_u
      real uhb                  ! U (wind component) Horizontal surfaces, B (explicit) term
      real vhb                  ! V (wind component) Horizontal surfaces, B (explicit) term
!     real thb                  ! Temperature        Horizontal surfaces, B (explicit) term
!     real tva                  ! Temperature          Vertical surfaces, A (implicit) term
!     real tvb                  ! Temperature          Vertical surfaces, B (explicit) term
      real ehb                  ! Energy (TKE)       Horizontal surfaces, B (explicit) term
      real sf


! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      real aa
      real al
      real buu
      real c
      real fbuw
      real fbpt
      real fh
      real fm
      real ric
      real tstar
      real ustar
      real utot
      real wstar
      real zz
      
      real b,cm,ch,rr,tol
      parameter(b=9.4,cm=7.4,ch=5.3,rr=0.74,tol=.001)

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------


! computation of the ground temperature
         
      utot=(ua**2+va**2)**.5
        
      
!!!! Louis formulation
!
! compute the bulk Richardson Number

      zz=dz/2.
   

      utot=max(utot,0.01)
          
      ric=2.*g_u*zz*(pt-ptg)/((pt+ptg)*(utot**2))
              
      aa=vk/log(zz/z0)

! determine the parameters fm and fh for stable, neutral and unstable conditions

      if(ric.gt.0)then
         fm=1/(1+0.5*b*ric)**2
         fh=fm
      else
         c=b*cm*aa*aa*(zz/z0)**.5
         fm=1-b*ric/(1+c*(-ric)**.5)
         c=c*ch/cm
         fh=1-b*ric/(1+c*(-ric)**.5)
      endif
      
      fbuw=-aa*aa*utot*utot*fm
      fbpt=-aa*aa*utot*(pt-ptg)*fh/rr
                     
      ustar=(-fbuw)**.5
      tstar=-fbpt/ustar

      al=(vk*g_u*tstar)/(pt*ustar*ustar)                      
      
      buu=-g_u/pt0*ustar*tstar
       
      uhb=-ustar*ustar*ua/utot
      vhb=-ustar*ustar*va/utot 
      sf= ustar*tstar*da*cp_u   
       
!     thb= 0.      
      ehb=buu
!!!!!!!!!!!!!!!
         
      return
      end subroutine flux_flat

! ===6=8===============================================================72
! ===6=8===============================================================72

      subroutine icBEP (fww,fwg,fgw,fsw,fws,fsg,                       &
                       z0g_u,z0r_u,                                    &
                       nd_u,strd_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u,        &
                       nz_u,z_u)                               


      implicit none 
      
   
!    Building parameters      

!    Radiation parameters

!    Roughness parameters
      real z0g_u(nurbm)       ! The ground's roughness length      
      real z0r_u(nurbm)       ! The roof's roughness length     

!    Street parameters
      integer nd_u(nurbm)     ! Number of street direction for each urban class

      real strd_u(ndm,nurbm)  ! Street length (fix to greater value to the horizontal length of the cells)
      real ws_u(ndm,nurbm)    ! Street width [m]
      real bs_u(ndm,nurbm)    ! Building width [m]
      real h_b(nz_um,nurbm)   ! Bulding's heights [m]
      real d_b(nz_um,nurbm)   ! The probability that a building has an height h_b
! -----------------------------------------------------------------------
!     Output
!------------------------------------------------------------------------



!   fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
!   and the short wave radation. They are the part of radiation from a surface
!   or from the sky to another surface.
      real fww(nz_um,nz_um,ndm,nurbm)        !  from wall to wall
      real fwg(nz_um,ndm,nurbm)              !  from wall to ground
      real fgw(nz_um,ndm,nurbm)              !  from ground to wall
      real fsw(nz_um,ndm,nurbm)              !  from sky to wall
      real fws(nz_um,ndm,nurbm)              !  from wall to sky
      real fsg(ndm,nurbm)                    !  from sky to ground

      real ss_u(nz_um,nurbm)     ! The probability that a building has an height equal to z
      real pb_u(nz_um,nurbm)     ! The probability that a building has an height greater or equal to z
        
!    Grid parameters
      integer nz_u(nurbm)     ! Number of layer in the urban grid
      real z_u(nz_um)       ! Height of the urban grid levels


! -----------------------------------------------------------------------
!     Local
!------------------------------------------------------------------------

      integer iz_u,id,ilu,iurb

      real dtot
      real hbmax

!------------------------------------------------------------------------
     

! -----------------------------------------------------------------------
!     This routine initialise the urban paramters for the BEP module
!------------------------------------------------------------------------
!
!Initialize some variables
!
       nz_u=0
       z_u=0.    
       ss_u=0.
       pb_u=0.
       fww=0.
       fwg=0.
       fgw=0.
       fsw=0.
       fws=0.
       fsg=0.

! Computation of the urban levels height

      z_u(1)=0.
     
      do iz_u=1,nz_um-1
         z_u(iz_u+1)=z_u(iz_u)+dz_u
      enddo
      
! Normalisation of the building density

      do iurb=1,nurbm
         dtot=0.
         do ilu=1,nz_um
            dtot=dtot+d_b(ilu,iurb)
         enddo
         do ilu=1,nz_um
            d_b(ilu,iurb)=d_b(ilu,iurb)/dtot
         enddo
      enddo      

! Compute the view factors, pb and ss 
      
      do iurb=1,nurbm         
         hbmax=0.
         nz_u(iurb)=0
         do ilu=1,nz_um
            if(h_b(ilu,iurb).gt.hbmax)hbmax=h_b(ilu,iurb)
         enddo
         
         do iz_u=1,nz_um-1
            if(z_u(iz_u+1).gt.hbmax)go to 10
         enddo
         
 10      continue
         nz_u(iurb)=iz_u+1

         do id=1,nd_u(iurb)

            call view_factors(iurb,nz_u(iurb),id,strd_u(id,iurb),   &    
                            z_u,ws_u(id,iurb),                      &    
                            fww,fwg,fgw,fsg,fsw,fws) 

            do iz_u=1,nz_u(iurb)
               ss_u(iz_u,iurb)=0.
               do ilu=1,nz_um
                  if(z_u(iz_u).le.h_b(ilu,iurb)                      &    
                    .and.z_u(iz_u+1).gt.h_b(ilu,iurb))then            
                        ss_u(iz_u,iurb)=ss_u(iz_u,iurb)+d_b(ilu,iurb)
                  endif 
               enddo
            enddo

            pb_u(1,iurb)=1.
            do iz_u=1,nz_u(iurb)
               pb_u(iz_u+1,iurb)=max(0.,pb_u(iz_u,iurb)-ss_u(iz_u,iurb))
            enddo

         enddo
      end do
     
                  
      return       
      end subroutine icBEP

! ===6=8===============================================================72
! ===6=8===============================================================72

      subroutine view_factors(iurb,nz_u,id,dxy,z,ws,fww,fwg,fgw,fsg,fsw,fws) 
     
      implicit none

 

! -----------------------------------------------------------------------
!     Input
!------------------------------------------------------------------------

      integer iurb            ! Number of the urban class
      integer nz_u            ! Number of levels in the urban grid
      integer id              ! Street direction number
      real ws                 ! Street width
      real z(nz_um)         ! Height of the urban grid levels
      real dxy                ! Street lenght


! -----------------------------------------------------------------------
!     Output
!------------------------------------------------------------------------

!   fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
!   and the short wave radation. They are the part of radiation from a surface
!   or from the sky to another surface.

      real fww(nz_um,nz_um,ndm,nurbm)            !  from wall to wall
      real fwg(nz_um,ndm,nurbm)                  !  from wall to ground
      real fgw(nz_um,ndm,nurbm)                  !  from ground to wall
      real fsw(nz_um,ndm,nurbm)                  !  from sky to wall
      real fws(nz_um,ndm,nurbm)                  !  from wall to sky
      real fsg(ndm,nurbm)                        !  from sky to ground


! -----------------------------------------------------------------------
!     Local
!------------------------------------------------------------------------

      integer jz,iz

      real hut
      real f1,f2,f12,f23,f123,ftot
      real fprl,fnrm
      real a1,a2,a3,a4,a12,a23,a123

! -----------------------------------------------------------------------
!     This routine calculates the view factors
!------------------------------------------------------------------------
        
      hut=z(nz_u+1)
        
      do jz=1,nz_u      
      
! radiation from wall to wall
       
         do iz=1,nz_u
     
            call fprls (fprl,dxy,abs(z(jz+1)-z(iz  )),ws)
            f123=fprl
            call fprls (fprl,dxy,abs(z(jz+1)-z(iz+1)),ws)
            f23=fprl
            call fprls (fprl,dxy,abs(z(jz  )-z(iz  )),ws)
            f12=fprl
            call fprls (fprl,dxy,abs(z(jz  )-z(iz+1)),ws)
            f2 = fprl
       
            a123=dxy*(abs(z(jz+1)-z(iz  )))
            a12 =dxy*(abs(z(jz  )-z(iz  )))
            a23 =dxy*(abs(z(jz+1)-z(iz+1)))
            a1  =dxy*(abs(z(iz+1)-z(iz  )))
            a2  =dxy*(abs(z(jz  )-z(iz+1)))
            a3  =dxy*(abs(z(jz+1)-z(jz  )))
       
            ftot=0.5*(a123*f123-a23*f23-a12*f12+a2*f2)/a1
       
            fww(iz,jz,id,iurb)=ftot*a1/a3

         enddo 

! radiation from ground to wall
       
         call fnrms (fnrm,z(jz+1),dxy,ws)
         f12=fnrm
         call fnrms (fnrm,z(jz)  ,dxy,ws)
         f1=fnrm
       
         a1 = ws*dxy
         
         a12= ws*dxy
       
         a4=(z(jz+1)-z(jz))*dxy
       
         ftot=(a12*f12-a12*f1)/a1
                    
         fgw(jz,id,iurb)=ftot*a1/a4
     
!  radiation from sky to wall
     
         call fnrms(fnrm,hut-z(jz)  ,dxy,ws)
         f12 = fnrm
         call fnrms (fnrm,hut-z(jz+1),dxy,ws)
         f1 =fnrm
       
         a1 = ws*dxy
       
         a12= ws*dxy
              
         a4 = (z(jz+1)-z(jz))*dxy
       
         ftot=(a12*f12-a12*f1)/a1
        
         fsw(jz,id,iurb)=ftot*a1/a4       
      
      enddo

! radiation from wall to sky      
      do iz=1,nz_u
       call fnrms(fnrm,ws,dxy,hut-z(iz))
       f12=fnrm
       call fnrms(fnrm,ws,dxy,hut-z(iz+1))
       f1=fnrm
       a1 = (z(iz+1)-z(iz))*dxy
       a2 = (hut-z(iz+1))*dxy
       a12= (hut-z(iz))*dxy
       a4 = ws*dxy
       ftot=(a12*f12-a2*f1)/a1
       fws(iz,id,iurb)=ftot*a1/a4 
 
      enddo
!!!!!!!!!!!!!


       do iz=1,nz_u

! radiation from wall to ground
      
         call fnrms (fnrm,ws,dxy,z(iz+1))
         f12=fnrm
         call fnrms (fnrm,ws,dxy,z(iz  ))
         f1 =fnrm
         
         a1= (z(iz+1)-z(iz) )*dxy
       
         a2 = z(iz)*dxy
         a12= z(iz+1)*dxy
         a4 = ws*dxy

         ftot=(a12*f12-a2*f1)/a1        
                    
         fwg(iz,id,iurb)=ftot*a1/a4
        
      enddo

! radiation from sky to ground
      
      call fprls (fprl,dxy,ws,hut)
      fsg(id,iurb)=fprl

      return
      end subroutine view_factors

! ===6=8===============================================================72
! ===6=8===============================================================72

      SUBROUTINE fprls (fprl,a,b,c)

      implicit none

     
            
      real a,b,c
      real x,y
      real fprl


      x=a/c
      y=b/c
      
      if(a.eq.0.or.b.eq.0.)then
       fprl=0.
      else
       fprl=log( ( (1.+x**2)*(1.+y**2)/(1.+x**2+y**2) )**.5)+  &
           y*((1.+x**2)**.5)*atan(y/((1.+x**2)**.5))+          &  
           x*((1.+y**2)**.5)*atan(x/((1.+y**2)**.5))-          &   
           y*atan(y)-x*atan(x)
       fprl=fprl*2./(pi*x*y)
      endif
      
      return
      end subroutine fprls

! ===6=8===============================================================72     
! ===6=8===============================================================72

      SUBROUTINE fnrms (fnrm,a,b,c)

      implicit none



      real a,b,c
      real x,y,z,a1,a2,a3,a4,a5,a6
      real fnrm
      
      x=a/b
      y=c/b
      z=x**2+y**2
      
      if(y.eq.0.or.x.eq.0)then
       fnrm=0.
      else
       a1=log( (1.+x*x)*(1.+y*y)/(1.+z) )
       a2=y*y*log(y*y*(1.+z)/z/(1.+y*y) )
       a3=x*x*log(x*x*(1.+z)/z/(1.+x*x) )
       a4=y*atan(1./y)
       a5=x*atan(1./x)
       a6=sqrt(z)*atan(1./sqrt(z))
       fnrm=0.25*(a1+a2+a3)+a4+a5-a6
       fnrm=fnrm/(pi*y)
      endif
      
      return
      end subroutine fnrms
  ! ===6=8===============================================================72  
     
      SUBROUTINE init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
        twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,&
        emr_u,emwind_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b, &
        cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, &
        targtemp_u, gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip)

! initialization routine, where the variables from the table are read

      implicit none
      
      integer iurb            ! urban class number
!    Building parameters      
      real alag_u(nurbm)      ! Ground thermal diffusivity [m^2 s^-1]
      real alaw_u(nurbm)      ! Wall thermal diffusivity [m^2 s^-1]
      real alar_u(nurbm)      ! Roof thermal diffusivity [m^2 s^-1]
      real csg_u(nurbm)       ! Specific heat of the ground material [J m^3 K^-1]
      real csw_u(nurbm)       ! Specific heat of the wall material [J m^3 K^-1]
      real csr_u(nurbm)       ! Specific heat of the roof material [J m^3 K^-1]
      real twini_u(nurbm)     ! Temperature inside the buildings behind the wall [K]
      real trini_u(nurbm)     ! Temperature inside the buildings behind the roof [K]
      real tgini_u(nurbm)     ! Initial road temperature

!    Radiation parameters
      real albg_u(nurbm)      ! Albedo of the ground
      real albw_u(nurbm)      ! Albedo of the wall
      real albr_u(nurbm)      ! Albedo of the roof
      real albwin_u(nurbm)    ! Albedo of the window
      real emg_u(nurbm)       ! Emissivity of ground
      real emw_u(nurbm)       ! Emissivity of wall
      real emr_u(nurbm)       ! Emissivity of roof
      real emwind_u(nurbm)    ! Emissivity of windows

!    Roughness parameters
      real z0g_u(nurbm)       ! The ground's roughness length      
      real z0r_u(nurbm)       ! The roof's roughness length

!    Street parameters
      integer nd_u(nurbm)     ! Number of street direction for each urban class

      real strd_u(ndm,nurbm)  ! Street length (fix to greater value to the horizontal length of the cells)
      real drst_u(ndm,nurbm)  ! Street direction [degree]
      real ws_u(ndm,nurbm)    ! Street width [m]
      real bs_u(ndm,nurbm)    ! Building width [m]
      real h_b(nz_um,nurbm)   ! Bulding's heights [m]
      real d_b(nz_um,nurbm)   ! The probability that a building has an height h_b

      integer i,iu
      integer nurb ! number of urban classes used

      real, intent(out) :: cop_u(nurbm)
      real, intent(out) :: pwin_u(nurbm)
      real, intent(out) :: beta_u(nurbm)
      integer, intent(out) :: sw_cond_u(nurbm)
      real, intent(out) :: time_on_u(nurbm)
      real, intent(out) :: time_off_u(nurbm)
      real, intent(out) :: targtemp_u(nurbm)
      real, intent(out) :: gaptemp_u(nurbm)
      real, intent(out) :: targhum_u(nurbm)
      real, intent(out) :: gaphum_u(nurbm)
      real, intent(out) :: perflo_u(nurbm)
      real, intent(out) :: hsesf_u(nurbm)
      real, intent(out) :: hsequip(24)

!
!We initialize
!
        h_b=0.
        d_b=0.

       nurb=ICATE
       do iu=1,nurb                         
          nd_u(iu)=0
       enddo

       csw_u=CAPB_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
       csr_u=CAPR_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
       csg_u=CAPG_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
       do i=1,icate
         alaw_u(i)=AKSB_TBL(i) / csw_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
         alar_u(i)=AKSR_TBL(i) / csr_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
         alag_u(i)=AKSG_TBL(i) / csg_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
       enddo
       twini_u=TBLEND_TBL
       trini_u=TRLEND_TBL
       tgini_u=TGLEND_TBL
       albw_u=ALBB_TBL
       albr_u=ALBR_TBL
       albg_u=ALBG_TBL
       emw_u=EPSB_TBL
       emr_u=EPSR_TBL
       emg_u=EPSG_TBL
       z0r_u=Z0R_TBL
       z0g_u=Z0G_TBL
       nd_u=NUMDIR_TBL
!MT BEM
       cop_u = cop_tbl
       pwin_u = pwin_tbl
       beta_u = beta_tbl
       sw_cond_u = sw_cond_tbl
       time_on_u = time_on_tbl
       time_off_u = time_off_tbl
       targtemp_u = targtemp_tbl
       gaptemp_u = gaptemp_tbl
       targhum_u = targhum_tbl
       gaphum_u = gaphum_tbl
       perflo_u = perflo_tbl
       hsesf_u = hsesf_tbl
       hsequip = hsequip_tbl

       do iu=1,icate
              if(ndm.lt.nd_u(iu))then
                write(*,*)'ndm too small in module_sf_bep, please increase to at least ', nd_u(iu)
                write(*,*)'remember also that num_urban_layers should be equal or greater than nz_um*ndm*nwr-u!'
                stop
              endif
         do i=1,nd_u(iu)
           drst_u(i,iu)=STREET_DIRECTION_TBL(i,iu) * pi/180.
           ws_u(i,iu)=STREET_WIDTH_TBL(i,iu)
           bs_u(i,iu)=BUILDING_WIDTH_TBL(i,iu)
         enddo
       enddo
       do iu=1,ICATE
          if(nz_um.lt.numhgt_tbl(iu)+3)then
              write(*,*)'nz_um too small in module_sf_bep, please increase to at least ',numhgt_tbl(iu)+3
              write(*,*)'remember also that num_urban_layers should be equal or greater than nz_um*ndm*nwr-u!'
              stop
          endif
         do i=1,NUMHGT_TBL(iu)
           h_b(i,iu)=HEIGHT_BIN_TBL(i,iu)
           d_b(i,iu)=HPERCENT_BIN_TBL(i,iu)
         enddo
       enddo

       do i=1,ndm
        do iu=1,nurbm
         strd_u(i,iu)=100000.
        enddo
       enddo

       do iu=1,nurb  
          emwind_u(iu)=0.9                       
          call albwindow(albwin_u(iu))  
       enddo
       
       return
       end subroutine init_para
! ===6================================================================72
! ===6================================================================72
      subroutine upward_rad(nd_u,nz_u,ws,bs,sigma,pb,ss,                     &
                       tg,emg_u,albg_u,rlg,rsg,sfg,                          & 
                       tw,emw_u,albw_u,rlw,rsw,sfw,                          & 
                       tr,emr_u,albr_u,emwind,albwind,twlev,pwin,            &
                       sfwind,rld,rs, sfr,                                   & 
                       rs_abs,rl_up,emiss,grdflx_urb)
!
! IN this surboutine we compute the upward longwave flux, and the albedo
! needed for the radiation scheme
!
      implicit none

!
!INPUT VARIABLES
!
      real rsw(2*ndm,nz_um)        ! Short wave radiation at the wall for a given canyon direction [W/m2]
      real rlw(2*ndm,nz_um)         ! Long wave radiation at the walls for a given canyon direction [W/m2]
      real rsg(ndm)                   ! Short wave radiation at the canyon for a given canyon direction [W/m2]
      real rlg(ndm)                   ! Long wave radiation at the ground for a given canyon direction [W/m2]
      real rs                        ! Short wave radiation at the horizontal surface from the sun [W/mÂ²]  
      real sfw(2*ndm,nz_um)      ! Sensible heat flux from walls [W/mÂ²]
      real sfg(ndm)              ! Sensible heat flux from ground (road) [W/mÂ²]
      real sfr(ndm,nz_um)      ! Sensible heat flux from roofs [W/mÂ²]                      
      real rld                        ! Long wave radiation from the sky [W/mÂ²]
      real albg_u                    ! albedo of the ground/street
      real albw_u                    ! albedo of the walls
      real albr_u                    ! albedo of the roof 
      real ws(ndm)                        ! width of the street
      real bs(ndm)
                        ! building size
      real pb(nz_um)                ! Probability to have a building with an height equal or higher   
      integer nz_u
      real ss(nz_um)                ! Probability to have a building of a given height
      real sigma                       
      real emg_u                       ! emissivity of the street
      real emw_u                       ! emissivity of the wall
      real emr_u                       ! emissivity of the roof
      real tw(2*ndm,nz_um)  ! Temperature in each layer of the wall [K]
      real tr(ndm,nz_um,nwr_u)  ! Temperature in each layer of the roof [K]
      real tg(ndm,ng_u)          ! Temperature in each layer of the ground [K]
      integer id ! street direction
      integer nd_u ! number of street directions
!
!New variables BEM
!
      real emwind               !Emissivity of the windows
      real albwind              !Albedo of the windows
      real twlev(2*ndm,nz_um)   !Averaged Temperature of the windows 
      real pwin                 !Coverage area fraction of the windows
      real gflwin               !Heat stored for the windows
      real sfwind(2*ndm,nz_um)  !Sensible heat flux from windows [W/mÂ²]

!OUTPUT/INPUT
      real rs_abs  ! absrobed solar radiationfor this street direction
      real rl_up   ! upward longwave radiation for this street direction
      real emiss ! mean emissivity
      real grdflx_urb ! ground heat flux 
!LOCAL
      integer iz,iw
      real rl_inc,rl_emit
      real gfl
      integer ix,iy,iwrong

         iwrong=1
      do iz=1,nz_u+1
      do id=1,nd_u
      do iw=1,nwr_u
        if(tr(id,iz,iw).lt.100.)then
              write(*,*)'in upward_rad ',iz,id,iw,tr(id,iz,iw) 
              iwrong=0
        endif
      enddo
      enddo
      enddo
           if(iwrong.eq.0)stop

      rl_up=0.
 
      rs_abs=0.
      rl_inc=0.
      emiss=0.
      rl_emit=0.
      grdflx_urb=0.
      do id=1,nd_u          
       rl_emit=rl_emit-( emg_u*sigma*(tg(id,ng_u)**4.)+(1-emg_u)*rlg(id))*ws(id)/(ws(id)+bs(id))/nd_u
       rl_inc=rl_inc+rlg(id)*ws(id)/(ws(id)+bs(id))/nd_u       
       rs_abs=rs_abs+(1.-albg_u)*rsg(id)*ws(id)/(ws(id)+bs(id))/nd_u
         gfl=(1.-albg_u)*rsg(id)+emg_u*rlg(id)-emg_u*sigma*(tg(id,ng_u)**4.)+sfg(id)
         grdflx_urb=grdflx_urb-gfl*ws(id)/(ws(id)+bs(id))/nd_u  
 
          do iz=2,nz_u
            rl_emit=rl_emit-(emr_u*sigma*(tr(id,iz,nwr_u)**4.)+(1-emr_u)*rld)*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
            rl_inc=rl_inc+rld*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u            
            rs_abs=rs_abs+(1.-albr_u)*rs*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
            gfl=(1.-albr_u)*rs+emr_u*rld-emr_u*sigma*(tr(id,iz,nwr_u)**4.)+sfr(id,iz)
            grdflx_urb=grdflx_urb-gfl*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
         enddo
           
         do iz=1,nz_u 
           
            rl_emit=rl_emit-(emw_u*(1.-pwin)*sigma*(tw(2*id-1,iz)**4.+tw(2*id,iz)**4.)+ &
                            (emwind*pwin*sigma*(twlev(2*id-1,iz)**4.+twlev(2*id,iz)**4.))+ &
                ((1.-emw_u)*(1.-pwin)+pwin*(1.-emwind))*(rlw(2*id-1,iz)+rlw(2*id,iz)))* &
                            dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u

            rl_inc=rl_inc+((rlw(2*id-1,iz)+rlw(2*id,iz)))*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u

            rs_abs=rs_abs+(((1.-albw_u)*(1.-pwin)+(1.-albwind)*pwin)*(rsw(2*id-1,iz)+rsw(2*id,iz)))*&
                          dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u 

            gfl=(1.-albw_u)*(rsw(2*id-1,iz)+rsw(2*id,iz)) +emw_u*( rlw(2*id-1,iz)+rlw(2*id,iz) )   &
             -emw_u*sigma*( tw(2*id-1,iz)**4.+tw(2*id,iz)**4. )+(sfw(2*id-1,iz)+sfw(2*id,iz))   

            gflwin=(1.-albwind)*(rsw(2*id-1,iz)+rsw(2*id,iz)) +emwind*(rlw(2*id-1,iz)+rlw(2*id,iz))   &
             -emwind*sigma*( twlev(2*id-1,iz)**4.+twlev(2*id,iz)**4.)+(sfwind(2*id-1,iz)+sfwind(2*id,iz)) 
               
           
            grdflx_urb=grdflx_urb-(gfl*(1.-pwin)+pwin*gflwin)*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u

         enddo
          
      enddo
        emiss=(emg_u+emw_u+emr_u)/3.
        rl_up=(rl_inc+rl_emit)-rld
       
         
      return

      END SUBROUTINE upward_rad

!====6=8===============================================================72         
!====6=8===============================================================72 
! ===6================================================================72
! ===6================================================================72

         subroutine albwindow(albwin)
		
!-------------------------------------------------------------------
	 implicit none


! -------------------------------------------------------------------
!Based on the 
!paper of J.Karlsson and A.Roos(2000):"modelling the angular behaviour
!of the total solar energy transmittance of windows"
!Solar Energy Vol.69,No.4,pp. 321-329.      
! -------------------------------------------------------------------
!Input
!-----	
        
!Output
!------
         real albwin	        ! albedo of the window  
!Local
!-----
	 real a,b,c		!Polynomial coefficients
	 real alfa,delta,gama	!Polynomial powers
	 real g0	        !transmittance when the angle 
                                !of incidence is normal to the surface.
         real asup,ainf
	 real fonc

!Constants
!--------------------
         
         real epsilon              !accuracy of the integration
         parameter (epsilon=1.e-07) 
         real n1,n2                !Index of refraction for glasses and air
         parameter(n1=1.,n2=1.5)
         integer intg,k
!--------------------------------------------------------------------		
         if (q_num.eq.0) then
           write(*,*) 'Category parameter of the windows no valid'
           stop
         endif

         g0=4.*n1*n2/((n1+n2)*(n1+n2))
	 a=8.
	 b=0.25/q_num
         c=1.-a-b	
	 alfa =5.2 + (0.7*q_num)
	 delta =2.
	 gama =(5.26+0.06*p_num)+(0.73+0.04*p_num)*q_num

         intg=1
!----------------------------------------------------------------------


100      asup=0.
         ainf=0.

         do k=1,intg
          call foncs(fonc,(pi*k/intg),a,b,c,alfa,delta,gama)
          asup=asup+(pi/intg)*fonc
         enddo

         intg=intg+1

         do k=1,intg
          call foncs(fonc,(pi*k/intg),a,b,c,alfa,delta,gama)
          ainf=ainf+(pi/intg)*fonc
         enddo
	 
         if(abs(asup-ainf).lt.epsilon) then
           albwin=1-g0+(g0/2.)*asup
         else
           goto 100
         endif
        
!---------------------------------------------------------------------- 	
	return
	end subroutine albwindow
!====================================================================72
!====================================================================72

        subroutine foncs(fonc,x,aa,bb,cc,alf,delt,gam)

        implicit none
!
        real x,aa,bb,cc
        real alf,delt,gam
        real fonc
  
        fonc=(((aa*(x**alf))/(pi**alf))+   &
             ((bb*(x**delt))/(pi**delt))+  &
             ((cc*(x**gam))/(pi**gam)))*sin(x)
        
        return
	end subroutine foncs
!====================================================================72
!====================================================================72  
END MODULE module_sf_bep_bem
