hvac_autosize.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #############################################################
SUBROUTINE hvac_autosize (B, BOP, G, T, TOP, KI,KLUOUT)
!     #############################################################
!!
!!    PURPOSE
!!    -------
!!
!!    Calibrates HVAC systems for TEB-BEM
!!      
!!    AUTHOR
!!    ------
!!
!!      G. Pigeon           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!    Original    05/2011
!!    modified    08/2013 add solar panels (V. Masson)
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_bem_n, ONLY : bem_t
USE modd_bem_option_n, ONLY : bem_options_t
USE modd_sfx_grid_n, ONLY : grid_t
USE modd_teb_irrig_n, ONLY : teb_irrig_t
USE modd_teb_n, ONLY : teb_t
USE modd_teb_option_n, ONLY : teb_options_t
USE modd_teb_irrig_n, ONLY : teb_irrig_t, teb_irrig_init
USE modd_diag_misc_teb_n, ONLY : diag_misc_teb_t, diag_misc_teb_init
USE modd_diag_misc_teb_options_n, ONLY : diag_misc_teb_options_t, diag_misc_teb_options_init
!
USE modd_csts,  ONLY : xcpd, xpi, xp00, xrd, xstefan
!
USE modi_teb
USE modi_sunpos
USE modi_sw_daycycle
USE modi_urban_lw_coef
USE modi_urban_solar_abs
USE modi_get_sizes_parallel
USE modi_diag_misc_teb_init_n
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
#ifdef AIX64 
USE omp_lib
#endif
!
IMPLICIT NONE
!
#ifndef AIX64
  include 'omp_lib.h'
#endif
!
!*    0.1    Declaration of arguments
!            ------------------------
!
TYPE(bem_t), INTENT(INOUT) :: B
TYPE(bem_options_t), INTENT(INOUT) :: BOP
TYPE(grid_t), INTENT(INOUT) :: G
TYPE(teb_t), INTENT(INOUT) :: T
TYPE(teb_options_t), INTENT(INOUT) :: TOP
!
INTEGER,       INTENT(IN)    :: KI     ! number of points
INTEGER,       INTENT(IN)    :: KLUOUT ! output listing logical unit
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
!-------------------------------------------------------------------------
!
TYPE(diag_misc_teb_t) :: YDMTC
TYPE(diag_misc_teb_t) :: YDMT
TYPE(diag_misc_teb_options_t) :: YDMTO
TYPE(teb_irrig_t)     :: YTIR
!
! Options TOP
!
LOGICAL :: GCANOPY
 CHARACTER(LEN=3) :: YBEM
 CHARACTER(LEN=6) :: YZ0H 
 CHARACTER(LEN=5) :: YCH_BEM  
!
! Teb Fields T
!
REAL, DIMENSION(KI) :: ZGD
!
! Bem Options BOP
!
 CHARACTER(LEN=6)    :: YCOOL_COIL
 CHARACTER(LEN=6)    :: YHEAT_COIL
!
! Bem!
 CHARACTER(LEN=4),DIMENSION(KI) :: YNATVENT
REAL, DIMENSION(KI) :: ZF_WATER_COND
REAL, DIMENSION(KI) :: ZAUX_MAX
REAL, DIMENSION(KI) :: ZT_WIN2
LOGICAL, DIMENSION(KI) :: GNATVENT_NIGHT
!
! Road irrigation (not used here)
!
TYPE(teb_irrig_t) :: YIR
!
!-------------------------------------------------------------------------
!
REAL, PARAMETER     :: ZTSTEP = 300.
INTEGER, PARAMETER  :: JPYEAR = 2004   ! Current year (UTC)
!
!local variable
!
INTEGER             :: IMONTH = 7  ! Current month (UTC)
INTEGER             :: IDAY   = 12 ! Current day (UTC)
REAL                :: ZTIME  = 0.  ! Time at start of the run (s)
!
 CHARACTER(LEN=3)    :: YIMPLICIT_WIND = 'NEW'
!
! Loop control indexes
!
INTEGER             :: JJ
INTEGER             :: JFORC_STEP
INTEGER             :: INB_STEP_ATM
!
! Intermediate variables
!
REAL, DIMENSION(KI) :: ZU_RF
REAL, DIMENSION(KI) :: ZU_WL
!
REAL,DIMENSION(KI) :: ZT1    ! intermediate variable
REAL,DIMENSION(KI) :: ZTN    ! intermediate variable
!
REAL, DIMENSION(KI) :: ZT_SKY
!
REAL, DIMENSION(KI) :: ZTOT_SW
REAL, DIMENSION(KI) :: ZTOUT_EQ
!
! Arguments to TEB
!
REAL, DIMENSION(KI) :: ZT_CAN
REAL, DIMENSION(KI) :: ZQ_CAN
REAL, DIMENSION(KI) :: ZU_CAN
REAL, DIMENSION(KI) :: ZZ_LOWCAN
!
REAL, DIMENSION(KI) :: ZPS
REAL, DIMENSION(KI) :: ZPA
REAL, DIMENSION(KI) :: ZEXNS
REAL, DIMENSION(KI) :: ZEXNA
REAL, DIMENSION(KI) :: ZTA
REAL, DIMENSION(KI) :: ZQA
REAL, DIMENSION(KI) :: ZRHOA
REAL, DIMENSION(KI) :: ZLW_RAD
REAL, DIMENSION(KI) :: ZRR
REAL, DIMENSION(KI) :: ZSR
REAL, DIMENSION(KI) :: ZZREF
!
REAL, DIMENSION(KI) :: ZPEW_A_COEF
REAL, DIMENSION(KI) :: ZPEW_B_COEF
!
REAL, DIMENSION(KI) :: ZTS_GD
!
REAL, DIMENSION(KI) :: ZDF_RF
REAL, DIMENSION(KI) :: ZDN_RF
REAL, DIMENSION(KI) :: ZDN_RD
REAL, DIMENSION(KI) :: ZDF_RD
!
REAL, DIMENSION(KI) :: ZQSAT_RF
REAL, DIMENSION(KI) :: ZQSAT_RD
REAL, DIMENSION(KI) :: ZDELT_RF
REAL, DIMENSION(KI) :: ZDELT_RD
!
 CHARACTER(LEN=4)    :: YSNOW_RF
 CHARACTER(LEN=4)    :: YSNOW_RD
!
REAL, DIMENSION(KI) :: ZLW_WA_TO_WB   ! longwave exchange coefficients
REAL, DIMENSION(KI) :: ZLW_WA_TO_R
REAL, DIMENSION(KI) :: ZLW_WB_TO_R
REAL, DIMENSION(KI) :: ZLW_WA_TO_NR
REAL, DIMENSION(KI) :: ZLW_WB_TO_NR
REAL, DIMENSION(KI) :: ZLW_R_TO_WA
REAL, DIMENSION(KI) :: ZLW_R_TO_WB
REAL, DIMENSION(KI) :: ZLW_G_TO_WA
REAL, DIMENSION(KI) :: ZLW_G_TO_WB
REAL, DIMENSION(KI) :: ZLW_S_TO_WA
REAL, DIMENSION(KI) :: ZLW_S_TO_WB
REAL, DIMENSION(KI) :: ZLW_S_TO_R
REAL, DIMENSION(KI) :: ZLW_S_TO_NR
REAL, DIMENSION(KI) :: ZLW_NR_TO_WA
REAL, DIMENSION(KI) :: ZLW_NR_TO_WB
REAL, DIMENSION(KI) :: ZLW_NR_TO_WIN
REAL, DIMENSION(KI) :: ZLW_WA_TO_WIN
REAL, DIMENSION(KI) :: ZLW_WB_TO_WIN
REAL, DIMENSION(KI) :: ZLW_G_TO_WIN
REAL, DIMENSION(KI) :: ZLW_R_TO_WIN
REAL, DIMENSION(KI) :: ZLW_S_TO_WIN
REAL, DIMENSION(KI) :: ZLW_WIN_TO_WA
REAL, DIMENSION(KI) :: ZLW_WIN_TO_WB
REAL, DIMENSION(KI) :: ZLW_WIN_TO_R
REAL, DIMENSION(KI) :: ZLW_WIN_TO_NR
!
REAL, DIMENSION(KI) :: ZLEW_RF     ! latent heat flux over roof (snow)
!
REAL, DIMENSION(KI) :: ZLEW_RD     ! latent heat flux over road (snow)
!
REAL, DIMENSION(KI) :: ZLE_WL_A    ! latent heat flux over wall
REAL, DIMENSION(KI) :: ZLE_WL_B    ! latent heat flux over wall
 
!
REAL, DIMENSION(KI) :: ZRNSNOW_RF  ! net radiation over snow
REAL, DIMENSION(KI) :: ZHSNOW_RF   ! sensible heat flux over snow
REAL, DIMENSION(KI) :: ZLESNOW_RF  ! latent heat flux over snow
REAL, DIMENSION(KI) :: ZGSNOW_RF   ! flux under the snow
REAL, DIMENSION(KI) :: ZMELT_RF    ! snow melt
!
REAL, DIMENSION(KI) :: ZRNSNOW_RD  ! net radiation over snow
REAL, DIMENSION(KI) :: ZHSNOW_RD   ! sensible heat flux over snow
REAL, DIMENSION(KI) :: ZLESNOW_RD  ! latent heat flux over snow
REAL, DIMENSION(KI) :: ZGSNOW_RD   ! flux under the snow
REAL, DIMENSION(KI) :: ZMELT_RD    ! snow melt
!
REAL, DIMENSION(KI) :: ZUW_RD      ! Momentum flux for roads
REAL, DIMENSION(KI) :: ZUW_RF      ! Momentum flux for roofs
REAL, DIMENSION(KI) :: ZDUWDU_RD   !
REAL, DIMENSION(KI) :: ZDUWDU_RF   !
REAL, DIMENSION(KI) :: ZUSTAR_TWN   ! friction velocity over town
!
REAL, DIMENSION(KI) :: ZCD           ! town averaged drag coefficient
REAL, DIMENSION(KI) :: ZCDN          ! town averaged neutral drag coefficient
REAL, DIMENSION(KI) :: ZCH_TWN      ! town averaged heat transfer coefficient
REAL, DIMENSION(KI) :: ZRI_TWN      ! town averaged Richardson number
REAL, DIMENSION(KI) :: ZRESA_TWN    ! town aerodynamical resistance
REAL, DIMENSION(KI) :: ZAC_RF      ! roof conductance
REAL, DIMENSION(KI) :: ZAC_RD       ! road conductance
REAL, DIMENSION(KI) :: ZAC_WL      ! wall conductance
REAL, DIMENSION(KI) :: ZAC_TOP       ! top conductance
REAL, DIMENSION(KI) :: ZAC_GD     ! garden conductance
REAL, DIMENSION(KI) :: ZAC_RF_WAT  ! roof water conductance
REAL, DIMENSION(KI) :: ZAC_RD_WAT  ! roof water conductance 
REAL, DIMENSION(KI) :: ZEMIT_LW_FAC
REAL, DIMENSION(KI) :: ZEMIT_LW_RD
REAL, DIMENSION(KI) :: ZT_RAD_IND
REAL, DIMENSION(KI) :: ZHU_BLD
REAL, DIMENSION(KI) :: ZTSUN
!
! Arguments to urban_solar_abs
!
REAL, DIMENSION(KI) :: ZZENITH
REAL, DIMENSION(KI) :: ZAZIM
!
REAL, DIMENSION(KI) :: ZDIR_SW
REAL, DIMENSION(KI) :: ZSCA_SW
!
REAL, DIMENSION(KI) :: ZFRAC_PANEL
REAL, DIMENSION(KI) :: ZALB_PANEL
REAL, DIMENSION(KI) :: ZSVF_GD
!
REAL, DIMENSION(KI) :: ZREC_SW_RD
REAL, DIMENSION(KI) :: ZREC_SW_SNOW_RD
REAL, DIMENSION(KI) :: ZREC_SW_WL_A
REAL, DIMENSION(KI) :: ZREC_SW_WL_B
REAL, DIMENSION(KI) :: ZREC_SW_GD
REAL, DIMENSION(KI) :: ZREC_SW_RF
!
REAL, DIMENSION(KI) :: ZDIR_ALB_TWN
REAL, DIMENSION(KI) :: ZSCA_ALB_TWN
REAL, DIMENSION(KI) :: ZSW_RAD_GD
REAL, DIMENSION(KI) :: ZREC_SW_WIN
REAL, DIMENSION(KI) :: ZREF_SW_GRND
REAL, DIMENSION(KI) :: ZREF_SW_FAC
!
REAL, DIMENSION(KI) :: ZALB_GR
REAL, DIMENSION(KI) :: ZALB_GD
!
! Arguments to urban_lw_coef
!
REAL, DIMENSION(KI) :: ZEMIS_GD
!
REAL, DIMENSION(KI) :: ZLW_WA_TO_G
REAL, DIMENSION(KI) :: ZLW_WB_TO_G
REAL, DIMENSION(KI) :: ZLW_S_TO_G
REAL, DIMENSION(KI) :: ZLW_WIN_TO_G
!
!! GREGOIRE 13/03
!
!
!new for shading
REAL, DIMENSION(KI) :: ZE_SHADING
LOGICAL, DIMENSION(KI) :: GSHAD_DAY
LOGICAL, DIMENSION(KI) :: GSHADE
!
!
! Case greenroof
REAL, DIMENSION(KI) :: ZUW_GR
REAL, DIMENSION(KI) :: ZRN_GR
REAL, DIMENSION(KI) :: ZH_GR
REAL, DIMENSION(KI) :: ZLE_GR
REAL, DIMENSION(KI) :: ZGFLUX_GR
REAL, DIMENSION(KI) :: ZRUNOFF_GR 
REAL, DIMENSION(KI) :: ZDRAIN_GR 
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!*      0.     Initialization
!              --------------
!
IF (lhook) CALL dr_hook('HVAC_AUTOSIZE',0,zhook_handle)
!
!    Design parameters
!
zps   = 101325.
zpa = zps
zexns = (zps/xp00)**(xrd/xcpd)
zexna = (zpa/xp00)**(xrd/xcpd)
zqa   = 0.011
zlw_rad= 300.
zrr = 0.0 
zsr = 0.0 
zzref =  50.
!
zpew_a_coef  = 0.5  
zpew_b_coef  = 0.5 
!
! initial value for air temperature and outdoor wall/roof/window/road temperature
zt_can(:) = 10.7/2 * sin(2*xpi/(24*3600) * (ztime+16*3600)) + (b%XT_SIZE_MAX(:)-10.7/2)
!
zq_can = 0.011
zu_can = 2.5 
zz_lowcan = zzref
!
!
zts_gd = 300.
!
zdn_rf = 0.0
zdf_rf = 1.0
zdn_rd = 0.0
zdf_rd = 1.0
!
zqsat_rf = 0.015
zqsat_rd = 0.015
zdelt_rf = 0.0 
zdelt_rd = 0.0 
!
!
!
!*      A.     Autosize of the heating system
!              ---------------------------------
!
zu_rf(:) = 0.0
DO jj=1,top%NROOF_LAYER
  zu_rf(:) = zu_rf(:) + t%XD_ROOF(:,jj)/t%XTC_ROOF(:,jj) 
END DO
zu_rf(:) = zu_rf(:) + 1./10. + 1./25.         
zu_rf(:) = 1. / zu_rf(:)
!
zu_wl(:) = 0.0
DO jj=1,top%NWALL_LAYER
  zu_wl(:) = zu_wl(:) + t%XD_WALL(:,jj)/t%XTC_WALL(:,jj)
END DO
zu_wl(:) = zu_wl(:) + 1./10. + 1./25.         
zu_wl(:) = 1. / zu_wl(:)
!
zrhoa = 1.30
zt_sky = 253.15
ztout_eq(:) = (b%XT_SIZE_MIN(:) + zt_sky(:))/2.
!
!   Heating Coil Capacity [W m-2(bld)]
b%XCAP_SYS_HEAT(:) = zu_wl(:) * t%XWALL_O_BLD(:) * (b%XTHEAT_TARGET(:) - ztout_eq(:)) &
               + b%XU_WIN(:)  * b%XGLAZ_O_BLD(:) * (b%XTHEAT_TARGET(:) - ztout_eq(:)) &
               + zu_rf(:)                        * (b%XTHEAT_TARGET(:) - ztout_eq(:)) &
               - b%XQIN(:) * t%XBLD_HEIGHT(:) / b%XFLOOR_HEIGHT(:)*          &
                 (1 - b%XQIN_FLAT(:))                                        &
               + b%XINF(:)    * t%XBLD_HEIGHT(:) / 3600* zrhoa(:) * xcpd *   &
                 (b%XTHEAT_TARGET(:) - b%XT_SIZE_MIN(:))                     &
               + b%XV_VENT(:) * t%XBLD_HEIGHT(:) / 3600* zrhoa(:) * xcpd *   &
                 (b%XTHEAT_TARGET(:) - b%XT_SIZE_MIN(:))
!
!   Rated air flow rate [kg s-1 m-2(bld)]
b%XM_SYS_RAT(:) = b%XCAP_SYS_HEAT(:)/xcpd/(323.15 - b%XTHEAT_TARGET(:))
!
!    Initialization
b%XCAP_SYS_RAT  = 0.
!
!*      B.     Autosize of the cooling system
!              -----------------------------------
!
zrhoa = 1.15
!
!-------------------------------------------------------
!
!    Initial values
!
! Options TOP
!
ybem = top%CBEM
top%CBEM = "BEM" 
!
yz0h = top%CZ0H
top%CZ0H = 'KAND07'
!
gcanopy = top%LCANOPY
top%LCANOPY = .false.
!
ych_bem = top%CCH_BEM
top%CCH_BEM = 'DOE-2'
!
!
! Teb Fields T
!
ALLOCATE(t%XWS_ROOF(ki))
ALLOCATE(t%XWS_ROAD(ki))
t%XWS_ROOF = 0.0
t%XWS_ROAD = 0.0
!
ALLOCATE(t%TSNOW_ROOF%WSNOW(ki,1))
ALLOCATE(t%TSNOW_ROOF%T    (ki,1))
ALLOCATE(t%TSNOW_ROOF%ALB  (ki))
ALLOCATE(t%TSNOW_ROOF%RHO  (ki,1))
ALLOCATE(t%TSNOW_ROOF%TS   (ki))
ALLOCATE(t%TSNOW_ROOF%EMIS (ki))
!
ysnow_rf = t%TSNOW_ROOF%SCHEME 
t%TSNOW_ROOF%SCHEME = 'NONE'
t%TSNOW_ROOF%WSNOW  = 0.0
t%TSNOW_ROOF%T      = 273.0
t%TSNOW_ROOF%RHO    = 0.0
t%TSNOW_ROOF%ALB    = 0.8
t%TSNOW_ROOF%TS     = 273.0
t%TSNOW_ROOF%EMIS   = 0.0
!
ALLOCATE(t%TSNOW_ROAD%WSNOW(ki,1))
ALLOCATE(t%TSNOW_ROAD%T    (ki,1))
ALLOCATE(t%TSNOW_ROAD%ALB  (ki))
ALLOCATE(t%TSNOW_ROAD%RHO  (ki,1))
ALLOCATE(t%TSNOW_ROAD%TS   (ki))
ALLOCATE(t%TSNOW_ROAD%EMIS (ki))
!
ysnow_rd = t%TSNOW_ROAD%SCHEME 
t%TSNOW_ROAD%SCHEME = 'NONE'
t%TSNOW_ROAD%WSNOW  = 0.0
t%TSNOW_ROAD%T      = 273.0
t%TSNOW_ROAD%RHO    = 0.0
t%TSNOW_ROAD%ALB    = 0.8
t%TSNOW_ROAD%TS     = 273.0
t%TSNOW_ROAD%EMIS   = 1.0
!
zgd = t%XGARDEN
t%XGARDEN = 0.
!
!------------------------------------------------
!BEM 
!
ycool_coil = bop%CCOOL_COIL
bop%CCOOL_COIL = 'IDEAL '
yheat_coil = bop%CHEAT_COIL
bop%CHEAT_COIL ='IDEAL '
!
!!TI_BLD
ALLOCATE(b%XTI_BLD(ki))
!! !ZTI_BLD   = 297.16 ! indoor air temperature
DO jj=1,ki
    b%XTI_BLD(jj) = max(b%XTHEAT_TARGET(jj),zt_can(jj)) ! indoor air temperature
ENDDO
!
!T_FLOOR, T_MASS
ALLOCATE(b%XT_FLOOR(ki,bop%NFLOOR_LAYER))
ALLOCATE(b%XT_MASS (ki,bop%NFLOOR_LAYER))
!
DO jj=1,bop%NFLOOR_LAYER
   b%XT_FLOOR(:,jj)  = b%XTI_BLD(:) ! building floor temperature
   b%XT_MASS(:,jj)   = b%XTI_BLD(:) ! building mass temperature
ENDDO
!
!OUTDOOR WINDOW TEMPERATURE
ALLOCATE(b%XT_WIN1(ki))
b%XT_WIN1(:) = zt_can(:)
!! 
ALLOCATE(b%XT_WIN2(ki))
b%XT_WIN2(:) = b%XTI_BLD(:)
!
ynatvent(:) = b%CNATVENT(:)
b%CNATVENT(:) = 'NONE'
zf_water_cond(:) = b%XF_WATER_COND(:)
b%XF_WATER_COND(:) = 0.
zaux_max = b%XAUX_MAX
b%XAUX_MAX = 0.
!
ALLOCATE(b%XQI_BLD(ki))
b%XQI_BLD = 0.011
!
gnatvent_night(:) = b%LNATVENT_NIGHT(:)
b%LNATVENT_NIGHT(:) = .false.
!
!
! TEB
!
ALLOCATE(t%XT_ROOF  (ki,top%NROOF_LAYER))
ALLOCATE(t%XT_ROAD  (ki,top%NROAD_LAYER))
ALLOCATE(t%XT_WALL_A(ki,top%NWALL_LAYER))
ALLOCATE(t%XT_WALL_B(ki,top%NWALL_LAYER))
!
!RF
t%XT_ROOF  (:,top%NROOF_LAYER)   = b%XTI_BLD(:)   ! roof layers temperatures
t%XT_ROOF(:,1) = zt_can(:)
zt1(:) = t%XT_ROOF(:,1)
ztn(:) = t%XT_ROOF(:,top%NROOF_LAYER)
IF (top%NROOF_LAYER .GT. 2) CALL interp_proftwl(zt1, ztn, t%XD_ROOF, t%XT_ROOF)
!
!RD
DO jj=1,top%NROAD_LAYER
   t%XT_ROAD(:,jj) = zt_can(:)
ENDDO
!
!WL_A
t%XT_WALL_A(:,top%NWALL_LAYER)   = b%XTI_BLD(:)   ! wall layers temperatures
t%XT_WALL_A(:,1) = zt_can(:)
zt1(:)=t%XT_WALL_A(:,1)
ztn(:)=t%XT_WALL_A(:,top%NWALL_LAYER)
IF (top%NWALL_LAYER .GT. 2) CALL interp_proftwl(zt1, ztn, t%XD_WALL, t%XT_WALL_A)
!
!WL_B
t%XT_WALL_B = t%XT_WALL_A
!
!
!* road watering (not used)
yir%LPAR_RD_IRRIG   = .false.
!
ALLOCATE(yir%XRD_START_MONTH(ki))
ALLOCATE(yir%XRD_END_MONTH  (ki))
ALLOCATE(yir%XRD_START_HOUR (ki))
ALLOCATE(yir%XRD_END_HOUR   (ki))
ALLOCATE(yir%XRD_24H_IRRIG  (ki))
!
yir%XRD_START_MONTH = 1.
yir%XRD_END_MONTH   = 1.
yir%XRD_START_HOUR  = 0.
yir%XRD_END_HOUR    = 24.
yir%XRD_24H_IRRIG   = 0.
!
!----------------------------------------------
!
! greenroofs   are not taken into account in the building's HVAC equipment sizing process
!
zrn_gr(:) = 0.
zh_gr(:) = 0.
zle_gr(:) = 0.
zgflux_gr(:) = 0.
zrunoff_gr(:) = 0.
zdrain_gr(:) = 0.
!
 CALL diag_misc_teb_init(ydmt)
!
zuw_gr(:) = 0.
!
zac_gd(:) = 0.
!
zfrac_panel  = 0.
zalb_panel   = 0.1
zsvf_gd(:) = 0.
zalb_gd(:) = 0.
zalb_gr(:) = 0.
gshade(:) = .false.
gshad_day(:) = .false.
ze_shading(:) = 0.
!
 CALL diag_misc_teb_options_init(ydmto)
 CALL diag_misc_teb_init(ydmtc)
!
 ydmto%LSURF_MISC_BUDGET = .true.
 CALL diag_misc_teb_init_n(ydmtc, ydmt, ydmto, top, ki, 0)
!
!* one supposes zero conduction heat flux between the greenroof and the nnroof.
ydmt%XG_GREENROOF_ROOF(:) = 0.
!
ydmt%XQIN(:) = b%XQIN(:)
!
ydmt%XT_SYS(:) = b%XTI_BLD(:)
ydmt%XQ_SYS(:) = b%XQI_BLD(:)
!
ydmt%XTHEAT_TARGET(:) = b%XTHEAT_TARGET(:)
ydmt%XTCOOL_TARGET(:) = b%XTCOOL_TARGET(:)
!
inb_step_atm = 3600*24*4/ztstep
DO jforc_step= 1,inb_step_atm
!
! Daily outdoor air temperature cycle
  zt_can(:) = 10.7/2 * sin(2*xpi/(24*3600) * (ztime+16*3600))  &
                 + (b%XT_SIZE_MAX(:)-10.7/2)
  zta(:) = zt_can(:)
!
!
!*      B.1     Solar radiation
!               ---------------
!
  CALL sunpos(jpyear, imonth, iday, ztime, g%XLON, g%XLAT, ztsun, zzenith, zazim)
!
  CALL sw_daycycle(ki, zzenith, ztot_sw)
!
  zdir_sw(:) = 0.88 * ztot_sw(:) * 0.85 ! manual adjustment
  zsca_sw(:) = 0.12 * ztot_sw(:) * 0.85 ! manual adjustment
  WHERE (zdir_sw < 0.0) 
    zdir_sw = 0.0
  END WHERE
  WHERE (zsca_sw < 0.0) 
    zsca_sw = 0.0
  END WHERE
!
!
! solar panels are not taken into account in the building's HVAC equipment sizing process
  CALL urban_solar_abs(top, t, b, ydmt, zdir_sw, zsca_sw, zzenith, zazim,&
                       zfrac_panel, zalb_panel, zalb_gd, zsvf_gd,&
                       zalb_gr, zdn_rf, zdf_rf, zdn_rd,          &
                       zdf_rd, zrec_sw_rd, zrec_sw_snow_rd,      &
                       zrec_sw_wl_a, zrec_sw_wl_b, zrec_sw_gd,   &
                       zrec_sw_rf, zdir_alb_twn, zsca_alb_twn,   &
                       zsw_rad_gd, zrec_sw_win, zref_sw_grnd,    &
                       zref_sw_fac, ze_shading, gshad_day, gshade  )
!
!
!*      B.2     LW properties
!               -------------
!
  zemis_gd = 1.0
  CALL urban_lw_coef(b, t, zlw_rad, zemis_gd,                       &
                     t%TSNOW_ROAD%TS, zts_gd,                       &  
                     zlw_wa_to_wb, zlw_wa_to_r, zlw_wb_to_r,            &
                     zlw_wa_to_nr,zlw_wb_to_nr,                         &
                     zlw_wa_to_g, zlw_wb_to_g,                          &
                     zlw_wa_to_win, zlw_wb_to_win,                      &
                     zlw_r_to_wa, zlw_r_to_wb, zlw_r_to_win,            &
                     zlw_g_to_wa, zlw_g_to_wb, zlw_g_to_win,            &
                     zlw_s_to_wa, zlw_s_to_wb, zlw_s_to_r,              &
                     zlw_s_to_nr, zlw_s_to_g, zlw_s_to_win,             &
                     zlw_win_to_wa, zlw_win_to_wb, zlw_win_to_r,        &
                     zlw_win_to_nr, zlw_win_to_g,                       &
                     zlw_nr_to_wa, zlw_nr_to_wb, zlw_nr_to_win          )
!
!*      B.3     TEB simulation
!               -------------
!
  CALL teb  (top, t, bop, b, yir, ydmt, yimplicit_wind, ztsun,               &
             zt_can, zq_can, zu_can, zt_can, zq_can, zu_can, zz_lowcan,      &
             zpew_a_coef, zpew_b_coef, zpew_a_coef, zpew_b_coef, zps, zpa,   &
             zexns, zexna, zta, zqa, zrhoa, zlw_rad, zrr, zsr, zzref, zzref, &
             zu_can, t%XH_TRAFFIC, t%XLE_TRAFFIC, ztstep, zdf_rf, zdn_rf,    &
             zdf_rd, zdn_rd, zqsat_rf, zqsat_rd, zdelt_rf, zdelt_rd, zts_gd, &
             zlew_rf, zuw_gr, zlew_rd, zle_wl_a, zle_wl_b,zrnsnow_rf,        &
             zhsnow_rf, zlesnow_rf, zgsnow_rf, zmelt_rf, zrn_gr, zh_gr,      &
             zle_gr, zgflux_gr, zdrain_gr, zrunoff_gr, zrnsnow_rd,           &
             zhsnow_rd, zlesnow_rd, zgsnow_rd, zmelt_rd, zuw_rd, zuw_rf,     &
             zduwdu_rd,  zduwdu_rf, zustar_twn, zcd, zcdn, zch_twn, zri_twn, &
             zresa_twn, zac_rf, zac_rd, zac_wl, zac_top, zac_gd, zac_rf_wat, &
             zac_rd_wat, zlw_wa_to_wb, zlw_wa_to_r, zlw_wb_to_r,             &
             zlw_wa_to_nr, zlw_wb_to_nr, zlw_r_to_wa, zlw_r_to_wb,           &
             zlw_g_to_wa, zlw_g_to_wb, zlw_s_to_wa, zlw_s_to_wb, zlw_s_to_r, &
             zlw_s_to_nr, zlw_nr_to_wa, zlw_nr_to_wb, zlw_nr_to_win,         &
             zlw_wa_to_win, zlw_wb_to_win, zlw_g_to_win, zlw_r_to_win,       &
             zlw_s_to_win, zlw_win_to_wa, zlw_win_to_wb, zlw_win_to_r,       &
             zlw_win_to_nr, iday, zemit_lw_fac, zemit_lw_rd, zt_rad_ind,     &
             zhu_bld, ztime, ze_shading  )
! 
!   Time update
  ztime = ztime + ztstep
  IF (ztime >= 86400) THEN
    ztime = 0.0
    iday = iday + 1
  END IF
!
ENDDO
!
b%XQIN(:) = ydmt%XQIN(:)
b%XTHEAT_TARGET(:) = ydmt%XTHEAT_TARGET(:)
b%XTCOOL_TARGET(:) = ydmt%XTCOOL_TARGET(:)
!
 CALL teb_irrig_init(yir)
 CALL diag_misc_teb_init(ydmt)
 CALL diag_misc_teb_init(ydmtc)
!----------------------------------------------------
!
! Options
!
top%CZ0H = yz0h
top%CBEM = ybem
top%LCANOPY = gcanopy
top%CCH_BEM = ych_bem
!
! Teb Fields T
!
DEALLOCATE(t%XWS_ROOF,t%XWS_ROAD)
!
t%TSNOW_ROOF%SCHEME = ysnow_rf
DEALLOCATE(t%TSNOW_ROOF%WSNOW)
DEALLOCATE(t%TSNOW_ROOF%T)
DEALLOCATE(t%TSNOW_ROOF%ALB)
DEALLOCATE(t%TSNOW_ROOF%RHO)
DEALLOCATE(t%TSNOW_ROOF%TS)
DEALLOCATE(t%TSNOW_ROOF%EMIS)
!
t%TSNOW_ROAD%SCHEME = ysnow_rd
DEALLOCATE(t%TSNOW_ROAD%WSNOW)
DEALLOCATE(t%TSNOW_ROAD%T)
DEALLOCATE(t%TSNOW_ROAD%ALB)
DEALLOCATE(t%TSNOW_ROAD%RHO)
DEALLOCATE(t%TSNOW_ROAD%TS)
DEALLOCATE(t%TSNOW_ROAD%EMIS)
!
t%XGARDEN = zgd 
!
!------------------------------------------------
!BEM 
!
bop%CCOOL_COIL = ycool_coil
bop%CHEAT_COIL = yheat_coil
!
!!TI_BLD
DEALLOCATE(b%XTI_BLD)
!
!T_FLOOR, T_MASS
DEALLOCATE(b%XT_FLOOR,b%XT_MASS)
!
!OUTDOOR WINDOW TEMPERATURE
DEALLOCATE(b%XT_WIN1)
DEALLOCATE(b%XT_WIN2)
!
b%CNATVENT(:) = ynatvent(:) 
b%XF_WATER_COND(:) = zf_water_cond(:)
b%XAUX_MAX = zaux_max
!
DEALLOCATE(b%XQI_BLD)
!
b%LNATVENT_NIGHT(:) = gnatvent_night(:)
!
! TEB
!
DEALLOCATE(t%XT_ROOF,t%XT_ROAD,t%XT_WALL_A,t%XT_WALL_B)
!
!--------------------------------------------------------
!
! -----------------------------------------------------------
! Print autosize results
! -----------------------------------------------------------
!
WRITE(kluout,*) ' '
WRITE(kluout,*) '    --------------------------------'
WRITE(kluout,*) '      HVAC AUTOSIZE CALCULATIONS '
WRITE(kluout,*) ' '
WRITE(kluout,*) '    Rated mass flow rate:'
WRITE(kluout,*) '    ',maxval(b%XM_SYS_RAT), 'kg s-1 m-2(bld)'
WRITE(kluout,*) '    ',minval(b%XM_SYS_RAT), 'kg s-1 m-2(bld)'
WRITE(kluout,*) '    Rated cooling system capacity:'
WRITE(kluout,*) '    ',maxval(b%XCAP_SYS_RAT), 'W m-2(bld)'
WRITE(kluout,*) '    ',minval(b%XCAP_SYS_RAT), 'W m-2(bld)'
WRITE(kluout,*) '    Rated heating sysem capacity:'
WRITE(kluout,*) '    ',maxval(b%XCAP_SYS_HEAT), 'W m-2(bld)'
WRITE(kluout,*) '    ',minval(b%XCAP_SYS_HEAT), 'W m-2(bld)'
WRITE(kluout,*) '    --------------------------------'
WRITE(kluout,*) ' '
IF (lhook) CALL dr_hook('HVAC_AUTOSIZE',1,zhook_handle)
!
CONTAINS 
!
SUBROUTINE interp_proftwl(PT1, PTN, PD, PT)
!interpolation of vertical profile for 'wall' : roof/wall
!arguments
REAL, DIMENSION(:), INTENT(IN)    :: PT1 !temperature layer 1
REAL, DIMENSION(:), INTENT(IN)    :: PTN !temperature layer N
REAL, DIMENSION(:,:), INTENT(IN)  :: PD  !depth of all layers
REAL, DIMENSION(:,:), INTENT(OUT) :: PT  !temperature of all layers
!local variables
INTEGER :: ILAYER ! number of layers
REAL, DIMENSION(SIZE(PT1)) :: ZDN ! total depth from mid layer 1 to mid layer n
REAL, DIMENSION(SIZE(PT1)) :: ZD  ! sequential depth in the calculation
INTEGER :: JJ, JI

ilayer=SIZE(pd,2)
DO ji=1,ki
   zdn(ji) = 0.5 * pd(ji,1)
   DO jj=2,ilayer-1
      zdn(ji) = zdn(ji) + pd(ji,jj)
   ENDDO
   zdn(ji) = zdn(ji) + 0.5 * pd(ji,ilayer)
ENDDO
DO ji=1,ki
   zd(ji) = 0.5*pd(ji,1)
   DO jj=2,ilayer-1
      zd(ji) = zd(ji) + 0.5*pd(ji,jj)
      pt(ji,jj) = pt1(ji) + (ptn(ji)-pt1(ji)) / zdn(ji) * zd(ji) 
      zd(ji) = zd(ji) + 0.5 * pd(ji,jj)
   ENDDO
   pt(ji,1) = pt1(ji)
   pt(ji,ilayer) = ptn(ji)
ENDDO
END SUBROUTINE interp_proftwl

END SUBROUTINE hvac_autosize