MODULE module_sf_bep
#if defined(mpas)
use mpas_atmphys_utilities, only: physics_error_fatal
#define FATAL_ERROR(M) call physics_error_fatal( M )
#else
use module_wrf_error
#define FATAL_ERROR(M) call wrf_error_fatal( M)
!USE module_model_constants
#endif
 USE module_sf_urban

! 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=18)

      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)

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

! 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
      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)

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




   CONTAINS
 
      subroutine BEP(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_ndm,  urban_map_zrd,  urban_map_zwd, urban_map_gd, &
                       urban_map_zd,  urban_map_zdf,   urban_map_bd, urban_map_wd, &
                      urban_map_gbd,  urban_map_fbd,                               &
                       num_urban_hi,                                               &
                      trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d,        &
                      sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d,      &
                      lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d,           &
                      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,                  &
                      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_ndm
   INTEGER, INTENT(IN  ) :: urban_map_zrd
   INTEGER, INTENT(IN  ) :: urban_map_zwd
   INTEGER, INTENT(IN  ) :: urban_map_gd
   INTEGER, INTENT(IN  ) :: urban_map_zd
   INTEGER, INTENT(IN  ) :: urban_map_zdf
   INTEGER, INTENT(IN  ) :: urban_map_bd
   INTEGER, INTENT(IN  ) :: urban_map_wd
   INTEGER, INTENT(IN  ) :: urban_map_gbd
   INTEGER, INTENT(IN  ) :: urban_map_fbd
   INTEGER, INTENT(IN  ) :: num_urban_hi
   REAL, DIMENSION( ims:ime, 1:urban_map_zrd, jms:jme ), INTENT(INOUT) :: trb_urb4d
   REAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw1_urb4d
   REAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw2_urb4d
   REAL, DIMENSION( ims:ime, 1:urban_map_gd , jms:jme ), INTENT(INOUT) :: tgb_urb4d
   REAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfw1_urb3d
   REAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfw2_urb3d
   REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: sfr_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ), INTENT(INOUT) :: sfg_urb3d
   REAL, DIMENSION( ims:ime, 1:num_urban_hi, jms:jme ), INTENT(IN) :: hi_urb2d
   REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: lp_urb2d
   REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: lb_urb2d
   REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: hgt_urb2d

!      integer nx,ny,nz              ! Number of points in the mesocsale grid
      real z(ims:ime,kms:kme,jms:jme)            ! Vertical coordinates
      REAL, INTENT(IN )::   DT      ! Time step
!      real zr(ims:ime,jms:jme)                ! Solar zenith angle
!      real deltar(ims:ime,jms:jme)            ! Declination of the sun
!      real ah(ims:ime,jms:jme)                ! Hour angle
!      real rs(ims:ime,jms:jme)                ! Solar radiation
!------------------------------------------------------------------------
!     Output
!------------------------------------------------------------------------ 
!      real tsk(ims:ime,jms:jme)           ! Average of surface temperatures (roads and roofs)
!
!    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(its:ite,jts:jte) ! upward long wave radiation
      real rs_abs(its:ite,jts:jte) ! absorbed short wave radiation
      real emiss(its:ite,jts:jte)  ! emissivity averaged for urban surfaces
      real grdflx_urb(its:ite,jts:jte)  ! ground heat flux for urban areas
!------------------------------------------------------------------------
!     Local
!------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      real hi_urb(its:ite,1:nz_um,jts:jte)    ! Height histograms of buildings
      real hi_urb1D(nz_um)                    ! Height histograms of buildings
      real hb_u(nz_um)                        ! Bulding's heights
      real ss_urb(nz_um)  ! Probability that a building has an height equal to z
      real pb_urb(nz_um)  ! Probability that a building has an height greater or equal to z
      integer nz_urb(nurbm)                   ! Number of layer in the urban grid
      integer nzurban(nurbm)                  

!    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
!
!   Building materials
!
      real csg(ng_u)                          ! Specific heat of the ground material [J m^3 K^-1]
      real csr(nwr_u)                         ! Specific heat of the roof material [J m^3 K^-1]
      real csw(nwr_u)                         ! Specific heat of the wall material [J m^3 K^-1]
      real alag(ng_u)                         ! Ground thermal diffusivity [m^2 s^-1]
      real alaw(nwr_u)                        ! Wall thermal diffusivity [m^2 s^-1]
      real alar(nwr_u)                        ! Roof thermal diffusivity [m^2 s^-1]
!
! for twini_u, and trini_u the initial value at the deepest level is kept constant during the simulation
!
!    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 emg_u(nurbm)                       ! Emissivity of ground
      real emw_u(nurbm)                       ! Emissivity of wall
      real emr_u(nurbm)                       ! Emissivity of roof

!   fww_u,fwg_u,fgw_u,fsw_u,fsg_u are the view factors used to compute the long wave
!   and the short wave radation. 
      real fww_u(nz_um,nz_um,ndm,nurbm)         !  from wall to wall
      real fwg_u(nz_um,ndm,nurbm)               !  from wall to ground
      real fgw_u(nz_um,ndm,nurbm)               !  from ground to wall
      real fsw_u(nz_um,ndm,nurbm)               !  from sky to wall
      real fws_u(nz_um,ndm,nurbm)               !  from sky to wall
      real fsg_u(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

!    Roughness parameters
      real z0(ndm,nz_um)           ! Roughness lengths "profiles"

!    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
!
!    Street parameters
!
      real bs(ndm)                 ! Building width 
      real ws(ndm)                 ! Street width
      real drst(ndm)               ! street directions 
      real strd(ndm)               ! Street lengths 
      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

!    Grid parameters

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

! 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)  ! 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
      real sfw1D(2*ndm,nz_um)      ! 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_e1D(kms:kme)               ! Explicit component of the TKE 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(nz_um,nwr_u,ndm)
      integer ind_gd(ng_u,ndm)
      integer ind_zd(nz_um,ndm)
!
      integer ix,iy,iz,iurb,id,iz_u,iw,ig,ir,ix1,iy1,k
      integer it, nint
      integer iii
      real time_h,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,                      &
           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 

!------------------------------------------------------------------------
!    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
!------------------------------------------------------------------------
!prepare the arrays to collapse indexes

      if(urban_map_zrd.lt.nz_um*ndm*nwr_u)then
        write(*,*)'urban_map_zrd too small, please increase to at least ', nz_um*ndm*nwr_u
        stop
      endif
      iii=0
      do iz_u=1,nz_um
      do iw=1,nwr_u
      do id=1,ndm
       iii=iii+1
       ind_zwd(iz_u,iw,id)=iii
      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 iz_u=1,nz_um
      do id=1,ndm
       iii=iii+1
       ind_zd(iz_u,id)=iii
      enddo
      enddo

      if (num_urban_hi.ge.nz_um)then
          write(*,*)'nz_um too small, please increase to at least ', num_urban_hi+1
          stop         
      endif

      do ix=its,ite
          do iy=jts,jte
              do iz_u=1,nz_um
                  hi_urb(ix,iz_u,iy)=0.
              enddo
          enddo
      enddo
      
      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 !iz
       do iz_u=1,num_urban_hi
          hi_urb(ix,iz_u,iy)= hi_urb2d(ix,iz_u,iy)
       enddo !iz_u
      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,emg_u,emw_u,&
                emr_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b)

!      Initialisation of the urban parameters and calculation of the view factors
             
         call icBEP(nd_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)

           hi_urb1D=0.
           do iz_u=1,nz_um
              hi_urb1D(iz_u)=hi_urb(ix,iz_u,iy)
           enddo

           call icBEPHI_XY(hb_u,hi_urb1D,ss_urb,pb_urb,    &
                           nz_urb(iurb),z_u)
           
           call param(iurb,nz_u(iurb),nz_urb(iurb),nzurban(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_urb,ss,pb_u,   &
                      pb_urb,pb,lp_urb2d(ix,iy),                     &
                      lb_urb2d(ix,iy),hgt_urb2d(ix,iy),FRC_URB2D(ix,iy))  
!
!We compute the view factors in the icBEP_XY routine
!           
         
           call icBEP_XY(iurb,fww_u,fwg_u,fgw_u,fsw_u,fws_u,fsg_u,   &
                         nd_u(iurb),strd,ws,nzurban(iurb),z_u)   

       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)
          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_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
!          tw1D(2*id-1,iz_u,iw)=tw1_u(ix,iy,ind_zwd(iz_u,iw,id))
!          tw1D(2*id,iz_u,iw)=tw2_u(ix,iy,ind_zwd(iz_u,iw,id))
            if(ind_zwd(iz_u,iw,id).gt.urban_map_zwd)write(*,*)'ind_zwd too big w',ind_zwd(iz_u,iw,id)
          tw1D(2*id-1,iz_u,iw)=tw1_urb4d(ix,ind_zwd(iz_u,iw,id),iy)
          tw1D(2*id,iz_u,iw)=tw2_urb4d(ix,ind_zwd(iz_u,iw,id),iy)
         enddo
         enddo
         enddo
	
         do id=1,ndm
          do ig=1,ng_u
!           tg1D(id,ig)=tg_u(ix,iy,ind_gd(ig,id))
            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)=tr_u(ix,iy,ind_zwd(iz_u,ir,id))
            if(ind_zwd(iz_u,ir,id).gt.urban_map_zwd)write(*,*)'ind_zwd too big r',ind_zwd(iz_u,ir,id)
            tr1D(id,iz_u,ir)=trb_urb4d(ix,ind_zwd(iz_u,ir,id),iy)
          enddo
          enddo
         enddo

         do id=1,ndm
	 do iz=1,nz_um
!	  sfw1D(2*id-1,iz)=sfw1(ix,iy,ind_zd(iz,id))
!	  sfw1D(2*id,iz)=sfw2(ix,iy,ind_zd(iz,id))
	  sfw1D(2*id-1,iz)=sfw1_urb3d(ix,ind_zd(iz,id),iy)
	  sfw1D(2*id,iz)=sfw2_urb3d(ix,ind_zd(iz,id),iy)
	 enddo
	 enddo
	 
	 do id=1,ndm
!	  sfg1D(id)=sfg(ix,iy,id)
	  sfg1D(id)=sfg_urb3d(ix,id,iy)
	 enddo
	 
	 do id=1,ndm
	 do iz=1,nz_um
!	  sfr1D(id,iz)=sfr(ix,iy,ind_zd(iz,id))
	  sfr1D(id,iz)=sfr_urb3d(ix,ind_zd(iz,id),iy)
	 enddo
	 enddo

         
         rs1D=swdown(ix,iy)
         rld1D=glw(ix,iy)
         time_h=(itimestep*dt)/3600.+gmt

         zr1D=acos(COSZ_URB2D(ix,iy))
         deltar1D=DECLIN_URB
         ah1D=OMG_URB2D(ix,iy)
!	 call angle(xlong(ix,iy),xlat(ix,iy),julday,time_h,zr1D,deltar1D,ah1D)
!        write(*,*) 'entro en BEP1D'
         call BEP1D(iurb,kms,kme,kts,kte,z1D,dt,ua1D,va1D,pt1D,da1D,pr1D,pt01D,  &
                   zr1D,deltar1D,ah1D,rs1D,rld1D,                   & 
                   alag,alaw,alar,csg,csw,csr,                      & 
                   albg_u(iurb),albw_u(iurb),albr_u(iurb),          &
                   emg_u(iurb),emw_u(iurb),emr_u(iurb),             & 
                   fww_u,fwg_u,fgw_u,fsw_u,                         &                    
                   fws_u,fsg_u,z0,                                  &                               
                   nd_u(iurb),strd,drst,ws,bs,ss,pb,                & 
                   nzurban(iurb),z_u,                               & 
                   tw1D,tg1D,tr1D,sfw1D,sfg1D,sfr1D,                & 
                   a_u1D,a_v1D,a_t1D,a_e1D,                         & 
                   b_u1D,b_v1D,b_t1D,b_e1D,                         & 
                   dlg1D,dl_u1D,sf1D,vl1D,rl_up(ix,iy),             &
                   rs_abs(ix,iy),emiss(ix,iy),grdflx_urb(ix,iy))                            
!        write(*,*) 'salgo de BEP1D'
         do id=1,ndm
	 do iz=1,nz_um
	  sfw1_urb3d(ix,ind_zd(iz,id),iy)=sfw1D(2*id-1,iz) 
	  sfw2_urb3d(ix,ind_zd(iz,id),iy)=sfw1D(2*id,iz) 
	 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_zd(iz,id),iy)=sfr1D(id,iz)
	 enddo
	 enddo
!
         do id=1,ndm
         do iz_u=1,nz_um
         do iw=1,nwr_u
          tw1_urb4d(ix,ind_zwd(iz_u,iw,id),iy)=tw1D(2*id-1,iz_u,iw)
          tw2_urb4d(ix,ind_zwd(iz_u,iw,id),iy)=tw1D(2*id,iz_u,iw)
         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_zwd(iz_u,ir,id),iy)=tr1D(id,iz_u,ir)
          enddo
          enddo
         enddo
        
          sf(ix,kts:kte,iy)=0.
          vl(ix,kts:kte,iy)=0.
          a_u(ix,kts:kte,iy)=0.
          a_v(ix,kts:kte,iy)=0.
          a_t(ix,kts:kte,iy)=0.
          a_e(ix,kts:kte,iy)=0.
          b_u(ix,kts:kte,iy)=0.
          b_v(ix,kts:kte,iy)=0.
          b_t(ix,kts:kte,iy)=0.
          b_e(ix,kts:kte,iy)=0.
          b_q(ix,kts:kte,iy)=0.
          dlg(ix,kts:kte,iy)=0.
          dl_u(ix,kts:kte,iy)=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)
          b_e(ix,iz,iy)=b_e1D(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
            
! ===6=8===============================================================72

      subroutine BEP1D(iurb,kms,kme,kts,kte,z,dt,ua,va,pt,da,pr,pt0,  &  
                      zr,deltar,ah,rs,rld,                            & 
                      alag,alaw,alar,csg,csw,csr,                     & 
                      albg,albw,albr,emg,emw,emr,                     & 
                      fww,fwg,fgw,fsw,fws,fsg,z0,                     &                                             
                      ndu,strd,drst,ws,bs,ss,pb,                      & 
                      nzu,z_u,                                        & 
                      tw,tg,tr,sfw,sfg,sfr,                           & 
                      a_u,a_v,a_t,a_e,                                &
                      b_u,b_v,b_t,b_e,                                & 
                      dlg,dl_u,sf,vl,rl_up,rs_abs,emiss,grdflx_urb)                             

! ----------------------------------------------------------------------
! 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 nz                 ! Number of vertical levels
      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 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(ng_u)            ! Ground thermal diffusivity [m^2 s^-1]
      real alaw(nwr_u)           ! Wall thermal diffusivity [m^2 s^-1]
      real alar(nwr_u)           ! Roof thermal diffusivity [m^2 s^-1]
      real csg(ng_u)             ! Specific heat of the ground material [J m^3 K^-1]
      real csw(nwr_u)            ! Specific heat of the wall material [J m^3 K^-1]
      real csr(nwr_u)            ! Specific heat of the roof material [J m^3 K^-1]

!    Radiation parameters
      real albg                  ! Albedo of the ground
      real albw                  ! Albedo of the wall
      real albr                  ! Albedo of the roof
      real emg                   ! Emissivity of ground
      real emw                   ! Emissivity of wall
      real emr                   ! 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 z0(ndm,nz_um)           ! Roughness lengths "profiles"
      
!    Street parameters
      integer ndu                  ! Number of street direction for each urban class 
      real strd(ndm)               ! Street length (set to a greater value then the horizontal length of the cells)
      real drst(ndm)               ! Street direction
      real ws(ndm)                 ! Street width
      real bs(ndm)                 ! Building width
      real ss(nz_um)               ! The probability that a building has an height equal to "z"
      real pb(nz_um)               ! The probability that a building has an height greater or equal to "z"
        
!    Grid parameters
      integer nzu                  ! Number of layer in the urban grid
      real z_u(nz_um)              ! Height of the urban grid levels


! ----------------------------------------------------------------------
! 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)  ! 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)      ! 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)     ! 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_e(kms:kme)              ! Explicit component of the TKE 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

! 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)     ! 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 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 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
      integer iz,id
      integer iw,ix,iy

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

! Fix some usefull parameters for the computation of the sources or sinks

      do iz=kts,kte
         dz(iz)=z(iz+1)-z(iz)
      end do

! Interpolation on the "urban grid"
      call interpol(kms,kme,kts,kte,nzu,z,z_u,ua,ua_u)
      call interpol(kms,kme,kts,kte,nzu,z,z_u,va,va_u)
      call interpol(kms,kme,kts,kte,nzu,z,z_u,pt,pt_u)
      call interpol(kms,kme,kts,kte,nzu,z,z_u,pt0,pt0_u)
      call interpol(kms,kme,kts,kte,nzu,z,z_u,pr,pr_u)
      call interpol(kms,kme,kts,kte,nzu,z,z_u,da,da_u)

                   
! Compute the modification of the radiation due to the buildings

      call modif_rad(iurb,ndu,nzu,z_u,ws,               &
                    drst,strd,ss,pb,                    &
                    tw,tg,albg,albw,emw,emg,            &
                    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(ndu,nzu,ws,bs,                   &
                       sigma,pb,ss,                     &
                       tg,emg,albg,rlg,rsg,sfg,         & 
                       tw,emw,albw,rlw,rsw,sfw,         & 
                       tr,emr,albr,rld,rs,sfr,          & 
                       rs_abs,rl_up,emiss,grdflx_urb)               
        
! Compute the surface temperatures
     

      call surf_temp(nzu,ndu,pr_u,dt,ss,                & 
                    rs,rld,rsg,rlg,rsw,rlw,             &
                    tg,alag,csg,emg,albg,ptg,sfg,gfg,   &
                    tr,alar,csr,emr,albr,ptr,sfr,gfr,   &
                    tw,alaw,csw,emw,albw,ptw,sfw,gfw)  
      
      
! Compute the implicit and explicit components of the sources or sinks on the "urban grid"
       
      call buildings(ndu,nzu,z0,ua_u,va_u,                       & 
                     pt_u,pt0_u,ptg,ptr,da_u,ptw,drst,           &                      
                     uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u,  & 
                     uhb_u,vhb_u,thb_u,ehb_u,ss,dt)                        
 
         
! 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.
!  
! 
 
      do id=1,ndu         
         sfg(id)=-da_u(1)*cp_u*thb_u(id,1)
         do iz=2,nzu
            sfr(id,iz)=-da_u(iz)*cp_u*thb_u(id,iz)
         enddo
         
         do iz=1,nzu
            sfw(2*id-1,iz)=-da_u(iz)*cp_u*(tvb_u(2*id-1,iz)+     &
                tva_u(2*id-1,iz)*pt_u(iz))
            sfw(2*id,iz)=-da_u(iz)*cp_u*(tvb_u(2*id,iz)+         &
                tva_u(2*id,iz)*pt_u(iz)) 
         enddo
      enddo
      
! calculation of the urban albedo and the upward long wave radiation
          
!!       call upward_rad(ndu,nzu,ws,bs,                      &
!!                       sigma,pb,ss,                        &
!!                       tg,emg,albg,rlg,rsg,sfg,            & 
!!                       tw,emw,albw,rlw,rsw,sfw,            & 
!!                       tr,emr,albr,rld,rs,sfr,             & 
!!                       rs_abs,rl_up,emiss,grdflx_urb)             

! Interpolation on the "mesoscale grid"

      call urban_meso(ndu,kms,kme,kts,kte,nzu,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,                          &
                     a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e)                    
       

! computation of the mean road temperature tsk (this value could be used 
! to replace the surface temperature in the radiation routines, if needed).

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

      call interp_length(ndu,kms,kme,kts,kte,nzu,z_u,z,ss,ws,bs,dlg,dl_u)

      return
      end subroutine BEP1D

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

       subroutine param(iurb,nzu,nzurb,nzurban,ndu,                   &
                       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_urb,ss,pb_u,   &
                       pb_urb,pb,lp_urb,lb_urb,hgt_urb,frc_urb)        

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

      implicit none

  
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      integer iurb                 ! Current urban class
      integer nzu                  ! Number of vertical urban levels in the current class
      integer nzurb                ! Number of vertical urban levels in the current class
      integer ndu                  ! 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"
      real ss_urb(nz_um)     ! The probability that a building has an height equal to "z"
      real pb_urb(nz_um)     ! The probability that a building has an height greater or equal to "z"
      real lp_urb                ! Building plan area density
      real lb_urb                ! Building surface area to plan area ratio
      real hgt_urb               ! Average building height weighted by building plan area [m]
      real frc_urb               ! Urban fraction
     
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
      real alag(ng_u)           ! Ground thermal diffusivity at each ground levels
      real alar(nwr_u)           ! Roof thermal diffusivity at each roof levels
      real alaw(nwr_u)           ! Wall thermal diffusivity at each wall levels
      real csg(ng_u)            ! Specific heat of the ground material at each ground levels
      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 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 equal
      integer nzurban

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

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------     
!
!Initialize
!
      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.
      nzurban=0
      
      ihu=0

       do iz=1,nz_um
          if (ss_urb(iz)/=0.) then
             ihu=1
             exit
          else
             continue
          endif
       enddo
           
       if (ihu==1) then
          do iz=1,nzurb+1
             ss(iz)=ss_urb(iz)
             pb(iz)=pb_urb(iz)
          enddo
          nzurban=nzurb
       else
          do iz=1,nzu+1
             ss(iz)=ss_u(iz,iurb)
             pb(iz)=pb_u(iz,iurb)
          end do 
          nzurban=nzu
       endif

       do id=1,ndu
        z0(id,1)=z0g_u(iurb)
        do iz=2,nzurban+1
           z0(id,iz)=z0r_u(iurb)
        enddo
       enddo
                 
       do ig=1,ng_u
        csg(ig)=csg_u(iurb)
        alag(ig)=alag_u(iurb)
       enddo
       
       do ir=1,nwr_u
        csr(ir)=csr_u(iurb)
        alar(ir)=alar_u(iurb)
       enddo
       
       do iw=1,nwr_u
        csw(iw)=csw_u(iurb)
        alaw(iw)=alaw_u(iurb)
       enddo
      
       do id=1,ndu
        strd(id)=strd_u(id,iurb)
        drst(id)=drst_u(id,iurb)     
       enddo
    
       do id=1,ndu
          if ((hgt_urb<=0.).OR.(lp_urb<=0.).OR.(lb_urb<=0.)) then
             ws(id)=ws_u(id,iurb)
             bs(id)=bs_u(id,iurb)
          else if ((lp_urb/frc_urb<1.).and.(lp_urb<lb_urb)) then
                  bs(id)=2.*hgt_urb*lp_urb/(lb_urb-lp_urb)
                  ws(id)=2.*hgt_urb*lp_urb*((frc_urb/lp_urb)-1.)/(lb_urb-lp_urb)
               else
                  ws(id)=ws_u(id,iurb)
                  bs(id)=bs_u(id,iurb) 
          endif
       enddo
       do id=1,ndu
          if ((bs(id)<=1.).OR.(bs(id)>=150.)) then
!            write(*,*) 'WARNING, WIDTH OF THE BUILDING WRONG',id,bs(id)
!            write(*,*) 'WIDTH OF THE STREET',id,ws(id)
             bs(id)=bs_u(id,iurb)
             ws(id)=ws_u(id,iurb)
          endif
          if ((ws(id)<=1.).OR.(ws(id)>=150.)) then
!            write(*,*) 'WARNING, WIDTH OF THE STREET WRONG',id,ws(id)
!            write(*,*) 'WIDTH OF THE BUILDING',id,bs(id)
             bs(id)=bs_u(id,iurb)
             ws(id)=ws_u(id,iurb)
          endif
       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_u+1)      ! Altitude of the cell interface
      real z_u(nz_um)       ! Altitude of the cell interface
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
!!    real c_u(nz_u)        ! Interpolated paramters in the "urban grid"
      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 modif_rad(iurb,nd,nz_u,z,ws,drst,strd,ss,pb,    &
                          tw,tg,albg,albw,emw,emg,               &
                          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,nwr_u) ! 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
! ----------------------------------------------------------------------
! 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,                 &
                      fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
         
         call short_rad(iurb,nz_u,id,albw,albg,fwg,fww,fgw,rsg,rsw,pb)
  
      enddo
      
      return
      end subroutine modif_rad


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

      subroutine surf_temp(nz_u,nd,pr,dt,ss,rs,rl,rsg,rlg,rsw,rlw,     &
                          tg,alag,csg,emg,albg,ptg,sfg,gfg,             &
                          tr,alar,csr,emr,albr,ptr,sfr,gfr,             &
                          tw,alaw,csw,emw,albw,ptw,sfw,gfw)             
                 

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

      implicit none
      
      
      
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      integer nz_u              ! Number of vertical layers defined in the urban grid
      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 alar(nwr_u)           ! Roof thermal diffusivity for the current urban class [m^2 s^-1]  
      real alaw(nwr_u)           ! Wall thermal diffusivity for the current urban class [m^2 s^-1]  
      real albg                 ! Albedo of the ground for the current urban class
      real albr                 ! Albedo of the roof for the current urban class
      real albw                 ! Albedo of the wall 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 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 dt                   ! Time step
      real emg                  ! Emissivity of ground for the current urban class
      real emr                  ! Emissivity of roof for the current urban class
      real emw                  ! Emissivity of wall for the current urban class
      real pr(nz_um)            ! Air pressure
      real rs                   ! Solar radiation
      real rl                   ! Downward flux of the longwave radiation
      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
      real sfg(ndm)             ! Sensible heat flux from ground (road)
      real sfr(ndm,nz_um)     ! Sensible heat flux from roofs
      real sfw(2*ndm,nz_um)     ! Sensible heat flux from walls
      real gfg(ndm)             ! Heat flux transferred from the surface of the ground (road) toward the interior
      real gfr(ndm,nz_um)     ! Heat flux transferred from the surface of the roof toward the interior
      real gfw(2*ndm,nz_um)     ! Heat flux transfered from the surface of the walls toward the interior
      real ss(nz_um)          ! Probability to have a building with height h
      real tg(ndm,ng_u)         ! Temperature in each layer of the ground [K]
      real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
      real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
      

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

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

      real rtg(ndm)             ! Total radiation at ground(road) surface (solar+incoming long+outgoing long)
      real rtr(ndm,nz_um)     ! Total radiation at roof surface (solar+incoming long+outgoing long)
      real rtw(2*ndm,nz_um)     ! Radiation at walls surface (solar+incoming long+outgoing long)
      real tg_tmp(ng_u)
      real tr_tmp(nwr_u)
      real tw_tmp(nwr_u)

      real dzg_u(ng_u)          ! Layer sizes in the ground

      real dzr_u(nwr_u)          ! Layers sizes in the roof
         
      real dzw_u(nwr_u)          ! Layer sizes in the wall
      

      data dzg_u /0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/
      data dzr_u /0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/   
      data dzw_u /0.02,0.02,0.02,0.02,0.02,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

!      Calculation for the roofs surfaces
         
       do iz=2,nz_u
            
        if(ss(iz).gt.0.)then
         do ir=1,nwr_u        
          tr_tmp(ir)=tr(id,iz,ir)
         end do
        
         call soil_temp(nwr_u,dzr_u,tr_tmp,ptr(id,iz),          &
                       alar,csr,rs,rl,pr(iz),dt,emr,albr,    &
                       rtr(id,iz),sfr(id,iz),gfr(id,iz))     
         do ir=1,nwr_u        
          tr(id,iz,ir)=tr_tmp(ir)
         end do
       
        end if
            
       end do !iz

!      Calculation for the walls surfaces
         
       do iz=1,nz_u
            
        do iw=1,nwr_u        
         tw_tmp(iw)=tw(2*id-1,iz,iw)
        end do
        call soil_temp(nwr_u,dzw_u,tw_tmp,ptw(2*id-1,iz),alaw,          &
                      csw,                                           &     
                      rsw(2*id-1,iz),rlw(2*id-1,iz),                 &     
                      pr(iz),dt,emw,                                 &    
                albw,rtw(2*id-1,iz),sfw(2*id-1,iz),gfw(2*id-1,iz))   

        do iw=1,nwr_u        
         tw(2*id-1,iz,iw)=tw_tmp(iw)
        end do
            
        do iw=1,nwr_u        
         tw_tmp(iw)=tw(2*id,iz,iw)
        end do
            
        call soil_temp(nwr_u,dzw_u,tw_tmp,ptw(2*id,iz),alaw,          &      
                      csw,                                         &     
                      rsw(2*id,iz),rlw(2*id,iz),                   &     
                      pr(iz),dt,emw,                               &     
               albw,rtw(2*id,iz),sfw(2*id,iz),gfw(2*id,iz))        
         do iw=1,nwr_u        
          tw(2*id,iz,iw)=tw_tmp(iw)
         end do
                   
        end do !iz
	
      end do !id
      
      return
      end subroutine surf_temp
     
! ===6=8===============================================================72     
! ===6=8===============================================================72  

      subroutine buildings(nd,nz,z0,ua_u,va_u,pt_u,pt0_u,         &
                        ptg,ptr,da_u,ptw,                            &
                        drst,uva_u,vva_u,uvb_u,vvb_u,                &
                        tva_u,tvb_u,evb_u,                           &
                        uhb_u,vhb_u,thb_u,ehb_u,ss,dt)                  

! ----------------------------------------------------------------------
! 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)     ! Walls potential temperatures 
      real ss(nz_um)          ! probability to have a building with height h
      real z0(ndm,nz_um)      ! Roughness lengths "profiles"
      real dt ! time step


! ----------------------------------------------------------------------
! 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
  
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      integer id,iz

! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
       dz=dz_u

      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),thb_u(id,1),ehb_u(id,1))            

!        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),thb_u(id,iz),ehb_u(id,iz))        
            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 iz=1,nz         
            call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz),     &  
                        ptw(2*id-1,iz),                             &   
                        uva_u(2*id-1,iz),vva_u(2*id-1,iz),          &   
                        uvb_u(2*id-1,iz),vvb_u(2*id-1,iz),          &   
                        tva_u(2*id-1,iz),tvb_u(2*id-1,iz),          &   
                        evb_u(2*id-1,iz),drst(id),dt)                  
                    
            call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz),    &   
                        ptw(2*id,iz),                              &    
                        uva_u(2*id,iz),vva_u(2*id,iz),             &    
                        uvb_u(2*id,iz),vvb_u(2*id,iz),             &    
                        tva_u(2*id,iz),tvb_u(2*id,iz),             &   
                        evb_u(2*id,iz),drst(id),dt) 
!
       
         end do
         
      end do
                
      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,                   &      
                             a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e)           

! ----------------------------------------------------------------------
!  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 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

! ----------------------------------------------------------------------
! 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_e(kms:kme)              ! Explicit component of the TKE sources or sinks
      
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
      real dzz
      real fact
      integer id,iz,iz_u
      real se,sr,st,su,sv
      real uet(kms:kme)                ! Contribution to TKE due to walls
      real veb,vta,vtb,vte,vtot,vua,vub,vva,vvb


! ----------------------------------------------------------------------
! 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.
         uet(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))
         
         do iz=kts,kte
            st=0.
            su=0.
            sv=0.
            se=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))
               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
         enddo
      enddo              

! vertical surface impact
           
      do iz=kts,kte 
         uet(iz)=0.
         do id=1,nd              
            vtb=0.
            vta=0.
            vua=0.
            vub=0.
            vva=0.
            vvb=0.
            veb=0.
	    vte=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   
               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
            enddo
           
            fact=1./vl(iz)/dz(iz)/nd
            b_t(iz)=b_t(iz)+vtb*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
         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,                  &
                         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,nwr_u)       ! Temperature in each layer of the wall [K]
      

! ----------------------------------------------------------------------
! 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)*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)*pb(i+1)
        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)*emw*sigma*fww(j,i,id,iurb)*       &
               tw(2*id,j,nwr_u)**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)*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)*pb(i-nz_u+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)*emw*sigma*fww(j,i-nz_u,id,iurb)*  &   
                tw(2*id-1,j,nwr_u)**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)*fwg(j,id,iurb)*pb(j+1)
       enddo
       
       do j=nz_u+1,2*nz_u
        aaa(2*nz_u+1,j)=-(1.-emw)*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)+emw*sigma*fwg(i,id,iurb)*pb(i+1)*    &
                      (tw(2*id-1,i,nwr_u)**4+tw(2*id,i,nwr_u)**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
!        cddz(nz+1)=ala(nz+1)/dz(nz)
       
      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=-cddz(nz)*(temp(nz)-temp(nz-1))*cs(nz)
       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,uva,vva,uvb,vvb,            &
                                  tva,tvb,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 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

! LOCAL:
! -----
      real hc
      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
         
      hc=5.678*(1.09+0.23*(vett/0.3048))  

      if(hc.gt.da*cp_u/dt)then
        hc=da*cp_u/dt
      endif
         
!         tvb=hc*ptw/da/cp_u
!         tva=-hc/da/cp_u
!!!!!!!!!!!!!!!!!!!!
! explicit 
      tvb=hc*ptw/da/cp_u-hc/da/cp_u*pt !c
      tva = 0.                  !c
         
      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,thb,ehb)
                                
! ----------------------------------------------------------------------
!           Calculation of the flux at the ground 
!           Formulation of Louis (Louis, 1979)       
! ----------------------------------------------------------------------

      implicit none

     

! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
      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

! ----------------------------------------------------------------------
! 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


! ----------------------------------------------------------------------
! 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.
   
!        if(tstar.lt.0.)then
!         wstar=(-ustar*tstar*g*hii/pt)**(1./3.)
!        else
!         wstar=0.
!        endif
!        
!      if (utot.le.0.7*wstar) utot=max(0.7*wstar,0.00001)

      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 
      thb=-ustar*tstar       
!      thb= 0.      
      ehb=buu
!!!!!!!!!!!!!!!
         
      return
      end subroutine flux_flat

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

      subroutine icBEP (nd_u,h_b,d_b,ss_u,pb_u,nz_u,z_u)                               

      implicit none       
        

!    Street parameters
      integer nd_u(nurbm)     ! Number of street direction for each urban class
      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
!------------------------------------------------------------------------

      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 variables
!
      z_u=0.
      nz_u=0
      ss_u=0.
      pb_u=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)

            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,emg_u,emw_u,&
                emr_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b)

! 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 emg_u(nurbm)       ! Emissivity of ground
      real emw_u(nurbm)       ! Emissivity of wall
      real emr_u(nurbm)       ! Emissivity of roof

!    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

!
!Initialize some variables
!
       
       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
       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 urban_map_zrd 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 urban_map_zrd 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
         
       return
      END SUBROUTINE init_para
!==============================================================

!==============================================================
      subroutine angle(along,alat,day,realt,zr,deltar,ah)
!     ----------------     
!      
!         Computation of the solar angles
!         schayes (1982,atm.  env. , p1407)
! Inputs
!========================
! along=longitud
! alat=latitude
! day=julian day (from the beginning of the year)
! realt= time GMT in hours
! Outputs
!============================
! zr=solar zenith angle
! deltar=declination angle
! ah=hour angle
!===============================

      implicit none
      real along,alat, realt, zr, deltar, ah, arg
      real rad,om,radh,initt, pii, drad, alongt, cphi, sphi
      real c1, c2, c3, s1, s2, s3, delta, rmsr2, cd, sid 
      real et, ahor, chor, coznt 
      integer day 


      data rad,om,radh,initt/0.0174533,0.0172142,0.26179939,0/
 
       zr=0.
       deltar=0.
       ah=0.

       pii = 3.14159265358979312
       drad = pii/180.
      
       alongt=along/15.
       cphi=cos(alat*drad)
       sphi=sin(alat*drad)
!
!     declination
!
       arg=om*day
       c1=cos(arg)
       c2=cos(2.*arg)
       c3=cos(3.*arg)
       s1=sin(arg)
       s2=sin(2.*arg)
       s3=sin(3.*arg)
       delta=0.33281-22.984*c1-0.3499*c2-0.1398*c3+3.7872*s1+0.03205*s2+0.07187*s3
       rmsr2=(1./(1.-0.01673*c1))**2
       deltar=delta*rad
       cd=cos(deltar)
       sid=sin(deltar)
!
!     time equation in hours
!
       et=0.0072*c1-0.0528*c2-0.0012*c3-0.1229*s1-0.1565*s2-0.0041*s3
!
!
!   hour angle  
!
      
!      ifh=0
      
 !     ahor=realt-12.+ifh+et+alongt
      ahor=realt-12.+et+alongt
      ah=ahor*radh
      chor=cos(ah)
!
!    zenith angle
!
      coznt=sphi*sid+cphi*cd*chor
      
       zr=acos(coznt)

      return

      END SUBROUTINE angle
!
!====6=8===============================================================72         
!====6=8===============================================================72 

      subroutine upward_rad(ndu,nzu,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,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/m2]  
      real sfw(2*ndm,nz_um)      ! Sensible heat flux from walls [W/m2]
      real sfg(ndm)              ! Sensible heat flux from ground (road) [W/m2]
      real sfr(ndm,nz_um)      ! Sensible heat flux from roofs [W/m2]                      
      real rld                        ! Long wave radiation from the sky [W/m2]
      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 nzu
      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,nwr_u)  ! 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 ndu ! number of street directions
!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,nzu+1
      do id=1,ndu
      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
        if(tw(2*id-1,iz,iw).lt.100.) then
           write(*,*)'in upward_rad ',iz,id,iw,tw(2*id-1,iz,iw) 
           iwrong=0
        endif
        if(tw(2*id,iz,iw).lt.100.) then
           write(*,*)'in upward_rad ',iz,id,iw,tw(2*id,iz,iw) 
           iwrong=0
        endif
      enddo
      enddo
      enddo
      do id=1,ndu
      do iw=1,ng_u
          if(tg(id,iw).lt.100.) then
            write(*,*)'in upward_rad ',id,iw,tg(id,iw) 
            iwrong=0
          endif
      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,ndu          
       rl_emit=rl_emit-( emg_u*sigma*(tg(id,ng_u)**4.)+(1-emg_u)*rlg(id))*ws(id)/(ws(id)+bs(id))/ndu
       rl_inc=rl_inc+rlg(id)*ws(id)/(ws(id)+bs(id))/ndu       
       rs_abs=rs_abs+(1.-albg_u)*rsg(id)*ws(id)/(ws(id)+bs(id))/ndu
         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))/ndu  
 
          do iz=2,nzu
            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))/ndu
            rl_inc=rl_inc+rld*ss(iz)*bs(id)/(ws(id)+bs(id))/ndu            
            rs_abs=rs_abs+(1.-albr_u)*rs*ss(iz)*bs(id)/(ws(id)+bs(id))/ndu
            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))/ndu
         enddo
           
         do iz=1,nzu            
            rl_emit=rl_emit-(emw_u*sigma*( tw(2*id-1,iz,nwr_u)**4.+tw(2*id,iz,nwr_u)**4. )+          &
               (1-emw_u)*( rlw(2*id-1,iz)+rlw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu
            rl_inc=rl_inc+(( rlw(2*id-1,iz)+rlw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu
            rs_abs=rs_abs+((1.-albw_u)*( rsw(2*id-1,iz)+rsw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu 
            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,nwr_u)**4.+tw(2*id,iz,nwr_u)**4. )+(sfw(2*id-1,iz)+sfw(2*id,iz))            
            grdflx_urb=grdflx_urb-gfl*dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu
         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=8===============================================================72
! ===6=8===============================================================72

      subroutine icBEP_XY(iurb,fww_u,fwg_u,fgw_u,fsw_u,             &
                          fws_u,fsg_u,ndu,strd,ws,nzu,z_u)                               

      implicit none       
        
!    Street parameters
      integer ndu     ! Number of street direction for each urban class
      integer iurb

      real strd(ndm)        ! Street length (fix to greater value to the horizontal length of the cells)
      real ws(ndm)          ! Street width [m]

!    Grid parameters
      integer nzu          ! Number of layer in the urban grid
      real z_u(nz_um)       ! Height of the urban grid levels
! -----------------------------------------------------------------------
!     Output
!------------------------------------------------------------------------

!   fww_u,fwg_u,fgw_u,fsw_u,fsg_u 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_u(nz_um,nz_um,ndm,nurbm)         !  from wall to wall
      real fwg_u(nz_um,ndm,nurbm)               !  from wall to ground
      real fgw_u(nz_um,ndm,nurbm)               !  from ground to wall
      real fsw_u(nz_um,ndm,nurbm)               !  from sky to wall
      real fws_u(nz_um,ndm,nurbm)               !  from sky to wall
      real fsg_u(ndm,nurbm)                     !  from sky to ground

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

      integer id

! -----------------------------------------------------------------------
!     This routine compute the view factors
!------------------------------------------------------------------------
!
!Initialize
!
      fww_u=0.
      fwg_u=0.
      fgw_u=0.
      fsw_u=0.
      fws_u=0.
      fsg_u=0.
      
      do id=1,ndu

            call view_factors(iurb,nzu,id,strd(id),z_u,ws(id),  &    
                              fww_u,fwg_u,fgw_u,fsg_u,fsw_u,fws_u) 
      
      enddo               
      return       
      end subroutine icBEP_XY
! ===6=8===============================================================72
! ===6=8===============================================================72

      subroutine icBEPHI_XY(hb_u,hi_urb1D,ss_u,pb_u,nzu,z_u)

      implicit none   
!-----------------------------------------------------------------------
!    Inputs
!-----------------------------------------------------------------------
!    Street parameters
!
      real hi_urb1D(nz_um)    ! The probability that a building has an height h_b
!
!     Grid parameters
!
      real z_u(nz_um)         ! Height of the urban grid levels
! -----------------------------------------------------------------------
!     Output
!------------------------------------------------------------------------

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

! -----------------------------------------------------------------------
!     Local
!------------------------------------------------------------------------
      real hb_u(nz_um)        ! Bulding's heights [m]
      integer iz_u,id,ilu

      real dtot
      real hbmax

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

!Initialize variables
!
      
      nzu=0
      ss_u=0.
      pb_u=0.
      
! Normalisation of the building density

         dtot=0.
         hb_u=0.

         do ilu=1,nz_um
            dtot=dtot+hi_urb1D(ilu)
         enddo
         
         do ilu=1,nz_um
            if (hi_urb1D(ilu)<0.) then
!              write(*,*) 'WARNING, HI_URB1D(ilu) < 0 IN BEP'
               go to 20
            endif
         enddo

         if (dtot.gt.0.) then
            continue
         else
!           write(*,*) 'WARNING, HI_URB1D <= 0 IN BEP'
            go to 20
         endif

         do ilu=1,nz_um
            hi_urb1D(ilu)=hi_urb1D(ilu)/dtot
         enddo
         
         hb_u(1)=dz_u   
         do ilu=2,nz_um
            hb_u(ilu)=dz_u+hb_u(ilu-1)
         enddo
           

! Compute pb and ss 
      
            
         hbmax=0.
       
         do ilu=1,nz_um
            if (hi_urb1D(ilu)>0.and.hi_urb1D(ilu)<=1.) then
                hbmax=hb_u(ilu)
            endif
         enddo
         
         do iz_u=1,nz_um-1
            if(z_u(iz_u+1).gt.hbmax)go to 10
         enddo

10       continue 
        
         nzu=iz_u+1

         if ((nzu+1).gt.nz_um) then 
             write(*,*) 'error, nz_um has to be increased to at least',nzu+1
             stop
         endif

            do iz_u=1,nzu
               ss_u(iz_u)=0.
               do ilu=1,nz_um
                  if(z_u(iz_u).le.hb_u(ilu)                      &    
                    .and.z_u(iz_u+1).gt.hb_u(ilu))then            
                        ss_u(iz_u)=ss_u(iz_u)+hi_urb1D(ilu)
                  endif 
               enddo
            enddo

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

20    continue    
      return
      end subroutine icBEPHI_XY
! ===6=8===============================================================72
! ===6=8===============================================================72
END MODULE module_sf_bep

      FUNCTION bep_nurbm () RESULT (bep_val_nurbm)
         USE module_sf_bep
         IMPLICIT NONE
         INTEGER :: bep_val_nurbm
         bep_val_nurbm = nurbm
      END FUNCTION bep_nurbm
      
      FUNCTION bep_ndm () RESULT (bep_val_ndm)
         USE module_sf_bep
         IMPLICIT NONE
         INTEGER :: bep_val_ndm
         bep_val_ndm = ndm
      END FUNCTION bep_ndm
      
      FUNCTION bep_nz_um () RESULT (bep_val_nz_um)
         USE module_sf_bep
         IMPLICIT NONE
         INTEGER :: bep_val_nz_um
         bep_val_nz_um = nz_um
      END FUNCTION bep_nz_um
      
      FUNCTION bep_ng_u () RESULT (bep_val_ng_u)
         USE module_sf_bep
         IMPLICIT NONE
         INTEGER :: bep_val_ng_u
         bep_val_ng_u = ng_u
      END FUNCTION bep_ng_u
      
      FUNCTION bep_nwr_u () RESULT (bep_val_nwr_u)
         USE module_sf_bep
         IMPLICIT NONE
         INTEGER :: bep_val_nwr_u
         bep_val_nwr_u = nwr_u
      END FUNCTION bep_nwr_u
      
