wall_layer_e_budget.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 wall_layer_e_budget(TOP, T, B, PT_WL, PTS_WL_B, PTI_WL_B, PTSTEP, PDN_RD,     &
                                   PRHOA, PAC_WL, PAC_BLD, PLW_RAD, PPS, PEXNS, PABS_SW_WL,  &
                                   PT_CANYON, PTS_RD, PTSNOW_RD, PTS_GD, PTS_FL,             &
                                   PLW_WA_TO_WB, PLW_R_TO_W, PLW_G_TO_W, PLW_NR_TO_W,        &
                                   PLW_WIN_TO_W, PLW_S_TO_W, PFLX_BLD_WL, PDQS_WL,           &
                                   PABS_LW_WL, PEMIT_LW_WL, PH_WL, PIMB_WL, PRADHT_IN,       &
                                   PRAD_RF_WL, PRAD_WL_WIN, PRAD_WL_FL, PRAD_WL_MA,          &
                                   PCONV_WL_BLD, PLOAD_IN_WL ) 
!   ##########################################################################
!
!!****  *ROAD_WALL_LAYER_E_BUDGET*  
!!
!!    PURPOSE
!!    -------
!
!     Computes the evoultion of roads and walls surface temperatures
!         
!     
!!**  METHOD
!     ------
!
!    6 : equations for evolution of Ts_road and Ts_wall simultaneously
!        *************************************************************
!
!     dTw_k(t) / dt = 1/(dw_k*Cw_k) * (- 2*Kw_k-1*(Tw_k-Tw_k-1)/(dw_k-1 +dw_k) 
!                                      - 2*Kw_k  *(Tw_k-Tw_k+1)/(dw_k+1 +dw_k) )
!
!     dTw_1(t) / dt = 1/(dw_1*Cw_1) * (  Rn_w - H_w - LE_w 
!                                      - 2*Kw_1*(Tw_1-Tw_2)/(dw_1 +dw_2)       )
!
!
!       with
!
!   K*_k  = (d*_k+ d*_k+1)/(d*_k/k*_k+ d*_k+1/k*_k+1)
!
!   Rn_w = abs_Rg_w 
!  - sigma * emis_w                                                   * Ts_w**4 (t+dt)
!  +         emis_w                       *      SVF_w                * LWR
!  + sigma * emis_w * emis_r              *      SVF_w                * Ts_r**4 (t+dt)
!  + sigma * emis_w * emis_w              * (1-2*SVF_w)               * Ts_w**4 (t+dt)
!  +         emis_w            (1-emis_r) *      SVF_r  *      SVF_w  * LWR
!  +         emis_w            (1-emis_w) *      SVF_w  * (1-2*SVF_w) * LWR
!  + sigma * emis_w * emis_w * (1-emis_r) *      SVF_w  * (1-  SVF_r) * Ts_w**4 (t+dt)
!  + sigma * emis_w * emis_w * (1-emis_w) * (1-2*SVF_w) * (1-2*SVF_w) * Ts_w**4 (t+dt)
!  + sigma * emis_w * emis_r * (1-emis_w) *      SVF_w  * (1-2*SVF_w) * Ts_r**4 (t+dt)
!
!  H_w  = rho Cp CH V ( Ts_w (t+dt) - Ta_canyon )
!
!  LE_w = rho Lv CH V ( qs_w (t+dt) - qa_canyon )
!
!
! The system is implicited (or semi-implicited).
!
! ZIMPL=1    ---> implicit system
! ZIMPL=0.5  ---> semi-implicit system
! ZIMPL=0    ---> explicit system
!
!
!
!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    MODD_CST
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      V. Masson           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    23/01/98 
!!                  21/11/01 (V. Masson and A. Lemonsu) bug of latent flux
!!                           for very strong evaporation (all reservoir emptied
!!                           in one time-step)
!!                     02/11 (V. Masson) splits the routine for road and walls separately
!!                     01/12 (V. Masson) separates the 2 walls
!!                     09/12 (G. Pigeon) modif internal convective coef convection
!!                     10/12 (G. Pigeon) add solar heat gain of indoor wall
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_teb_option_n, ONLY : teb_options_t
USE modd_teb_n, ONLY : teb_t
USE modd_bem_n, ONLY : bem_t
!
USE modd_surf_par, ONLY : xundef
USE modd_csts,ONLY : xcpd, xstefan
!
USE modi_layer_e_budget_get_coef
USE modi_layer_e_budget
USE mode_conv_doe
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
TYPE(teb_options_t), INTENT(INOUT) :: TOP
TYPE(teb_t), INTENT(INOUT) :: T
TYPE(bem_t), INTENT(INOUT) :: B
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: PT_WL    ! wall layers temperatures
REAL, DIMENSION(:),   INTENT(IN)  :: PTS_WL_B ! opposite wall surface temperature
REAL, DIMENSION(:),   INTENT(IN)  :: PTI_WL_B ! opposite wall internal temperature
REAL,               INTENT(IN)    :: PTSTEP       ! time step
REAL, DIMENSION(:), INTENT(IN)    :: PDN_RD     ! snow-covered fraction on roads
REAL, DIMENSION(:), INTENT(IN)    :: PRHOA        ! rho
REAL, DIMENSION(:), INTENT(IN)    :: PAC_WL     ! aerodynamical conductance [m/s]
!                                                 ! between wall and canyon
REAL, DIMENSION(:), INTENT(IN)    :: PAC_BLD      ! aerodynamical conductance
                                                  ! inside the building itself
REAL, DIMENSION(:), INTENT(IN)    :: PLW_RAD      ! atmospheric infrared radiation
REAL, DIMENSION(:), INTENT(IN)    :: PPS          ! pressure at the surface
REAL, DIMENSION(:), INTENT(IN)    :: PEXNS        ! surface Exner function
REAL, DIMENSION(:), INTENT(IN)    :: PABS_SW_WL ! absorbed solar radiation
REAL, DIMENSION(:), INTENT(IN)    :: PT_CANYON    ! air canyon temperature
REAL, DIMENSION(:), INTENT(IN)    :: PTS_RD     ! road surface temperature
REAL, DIMENSION(:), INTENT(IN)    :: PTSNOW_RD  ! road snow temperature
REAL, DIMENSION(:), INTENT(IN)    :: PTS_GD   ! green area surface temperature
REAL, DIMENSION(:), INTENT(IN)    :: PTS_FL    ! floor layers temperatures [K]
!
REAL, DIMENSION(:), INTENT(IN)    :: PLW_WA_TO_WB ! LW interactions wall  -> opposite wall
REAL, DIMENSION(:), INTENT(IN)    :: PLW_R_TO_W   ! LW interactions road -> wall ; DEF formulation 
REAL, DIMENSION(:), INTENT(IN)    :: PLW_G_TO_W   ! LW interactions garden -> wall ; DEF formulation
REAL, DIMENSION(:), INTENT(IN)    :: PLW_S_TO_W   ! LW interactions sky   -> wall 
REAL, DIMENSION(:), INTENT(IN)    :: PLW_NR_TO_W  ! LW interactions road(snow) -> wall 
REAL, DIMENSION(:), INTENT(IN)    :: PLW_WIN_TO_W ! Radiative heat trasfer coeff wall-window 
                                                  ! [W K-1 m-2]
!
REAL, DIMENSION(:), INTENT(OUT)   :: PFLX_BLD_WL! flux from bld to wall
REAL, DIMENSION(:), INTENT(INOUT) :: PDQS_WL    ! heat storage inside the wall 
REAL, DIMENSION(:), INTENT(OUT)   :: PABS_LW_WL ! absorbed infrared rad. [W m-2(wall)]
REAL, DIMENSION(:), INTENT(OUT)   :: PEMIT_LW_WL  ! LW flux emitted by the wall [W m-2(wall)]
REAL, DIMENSION(:), INTENT(OUT)   :: PH_WL      ! Sensible heat flux from wall to air [W/m2(wall)]
                                                  ! wall = facade - glazing
REAL, DIMENSION(:), INTENT(OUT)   :: PIMB_WL    ! wall residual energy imbalance 
                                                  ! for verification [W m-2]
REAL, DIMENSION(:),   INTENT(IN)  :: PRADHT_IN     ! Indoor radiant heat transfer coefficient
                                                    ! [W K-1 m-2]
REAL, DIMENSION(:), INTENT(IN)    :: PRAD_RF_WL ! rad. fluxes from roof to wall [W m-2(roof)]
REAL, DIMENSION(:), INTENT(OUT)   :: PRAD_WL_WIN  ! rad. fluxes from wall to win  [W m-2(wall)]
REAL, DIMENSION(:), INTENT(OUT)   :: PRAD_WL_FL! rad. fluxes from wall to floor [W m-2(wall)]
REAL, DIMENSION(:), INTENT(OUT)   :: PRAD_WL_MA ! rad. fluxes from wall to mass [W m-2(wall)]
REAL, DIMENSION(:), INTENT(OUT)   :: PCONV_WL_BLD ! conv. fluxes from wall to bld [W m-2(wall)]

REAL, DIMENSION(:), INTENT(IN)    :: PLOAD_IN_WL  ! LOAD from solar heat gain + rad int. gains  W/m2 [Wall]


!
!*      0.2    declarations of local variables
!
!
REAL :: ZIMPL=1.0      ! implicit coefficient
REAL :: ZEXPL=0.0      ! explicit coefficient
!
REAL, DIMENSION(SIZE(PPS),SIZE(PT_WL,2)) ::  ZA,& ! lower diag.
                                               ZB,& ! main  diag.
                                               ZC,& ! upper diag.
                                               ZY   ! r.h.s.                       
!
REAL, DIMENSION(SIZE(PPS)) :: ZMTC_O_D_WL_IN
REAL, DIMENSION(SIZE(PPS)) :: ZDF_RD    ! Road snow free fraction
REAL, DIMENSION(SIZE(PPS)) :: ZRHO_ACF_W  ! rho * conductance
!                                         !     * snow-free f.
!
! thermal capacity times layer depth
REAL, DIMENSION(SIZE(PPS)) :: ZTS_WL       ! wall surface temperature
REAL, DIMENSION(SIZE(PPS)) :: ZTI_WL       ! wall indoor surface temperature
REAL, DIMENSION(SIZE(PPS)) :: ZTI_WL_CONV  ! wall indoor surface temperature for conv. flux
REAL, DIMENSION(SIZE(PPS)) :: ZT_SKY         ! sky temperature [K]
!
REAL, DIMENSION(SIZE(PPS)) :: ZTI_RF       ! Indoor roof temperature [K]
REAL, DIMENSION(SIZE(PPS)) :: ZDIF_RAD_WL_RF !diff between the rad flux that should receive the wall from the roof
                                                 ! and what it really receives [W m-2(bld)]
REAL, DIMENSION(SIZE(PPS)) :: ZRAD_WL_RF     ! rad flux between the wall and the roof computed for the wall balance
REAL, DIMENSION(SIZE(PPS)) :: ZF_WL_WL       ! View factor wall-wall inside the building
REAL, DIMENSION(SIZE(PPS)) :: ZCHTC_IN_WL  ! indoor convective heat transfer coeff wall [W m-2 K-1]
INTEGER :: IWL_LAYER           ! number of wall layers
INTEGER :: JJ                    ! loop counter
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('WALL_LAYER_E_BUDGET',0,zhook_handle)
!
pflx_bld_wl  = xundef
pimb_wl      = xundef
prad_wl_win  = xundef
prad_wl_fl= xundef
prad_wl_ma = xundef
pconv_wl_bld = xundef

!
 CALL layer_e_budget_get_coef( pt_wl, ptstep, zimpl, t%XHC_WALL, t%XTC_WALL, t%XD_WALL, &
                              za, zb, zc, zy )
!
iwl_layer = SIZE(pt_wl,2)
!
DO jj=1,SIZE(pdn_rd)
  !
  zdf_rd(jj) = 1. - pdn_rd(jj)
  !
  zts_wl(jj) = pt_wl(jj,1)
  zti_wl(jj) = pt_wl(jj, iwl_layer)
  !
  !*      2.1    outdoor convective flux properties 
  !              ----------------------------------
  !
  zrho_acf_w(jj) = prhoa(jj) * pac_wl(jj)
  !
  !*      2.2    Sky temperature
  !              ---------------
  !
  zt_sky(jj) = (plw_rad(jj)/xstefan)**0.25
  !  
  !*      2.3    indoor average thermal conductivity
  !              -----------------------------------
  !
  IF (top%CBEM .EQ. "DEF") THEN
    zmtc_o_d_wl_in(jj) = 2. * t%XTC_WALL(jj,iwl_layer) / t%XD_WALL (jj,iwl_layer)
    zmtc_o_d_wl_in(jj) = 1./(  1./zmtc_o_d_wl_in(jj)  + 1./(xcpd*prhoa(jj)*pac_bld(jj)))
  ENDIF
ENDDO  
  !  
  !*      2.4    indoor convective coefficient
  !              -----------------------------
  !
  zchtc_in_wl(:) = chtc_vert_doe(pt_wl(:,iwl_layer), b%XTI_BLD(:))
  DO jj=1,SIZE(zchtc_in_wl)
     zchtc_in_wl(jj) = max(1., zchtc_in_wl(jj))
  ENDDO

!
!-------------------------------------------------------------------------------
!
!*      3.    Outer wall layer coefficients
!             ------------------------------
!
DO jj=1,SIZE(pt_wl,1)
  !
  zb(jj,1) = zb(jj,1) + zimpl * xcpd/pexns(jj) * zrho_acf_w(jj)
  !
  zy(jj,1) = zy(jj,1) + pabs_sw_wl(jj)  &
                      + xcpd/pexns(jj) * zrho_acf_w(jj) * ( pt_canyon(jj) - zexpl * zts_wl(jj) )
  !
  !
  zb(jj,1) = zb(jj,1) &
             + zimpl * ( plw_s_to_w(jj) + plw_wa_to_wb(jj)                &
                       + zdf_rd(jj)*plw_r_to_w(jj) +  plw_g_to_w(jj)    &
                       + pdn_rd(jj) *   plw_nr_to_w(jj)                 &
                       + plw_win_to_w(jj)  )    
  !
  zy(jj,1) = zy(jj,1) + &
                      plw_s_to_w(jj) * (zt_sky(jj)     - zexpl * zts_wl(jj))  &
                    + plw_wa_to_wb(jj) * (pts_wl_b(jj) - zexpl * zts_wl(jj))  &
       + zdf_rd(jj) *  plw_r_to_w(jj)  * (pts_rd(jj)   - zexpl * zts_wl(jj))  &
       + pdn_rd(jj) *  plw_nr_to_w(jj) * (ptsnow_rd(jj)- zexpl * zts_wl(jj))  &
                + plw_win_to_w(jj) * (b%XT_WIN1(jj)    - zexpl * zts_wl(jj))
  !
  IF (SIZE(pts_gd)>0) THEN
    zy(jj,1) = zy(jj,1) + plw_g_to_w(jj)  * (pts_gd(jj)   - zexpl * zts_wl(jj))
  ENDIF 
  !
ENDDO
!
!-------------------------------------------------------------------------------
!
!*      4.    Inside wall layer coefficients
!             -----------------------------
!
DO jj=1,SIZE(pt_wl,1)
  !                
  IF (top%CBEM=="DEF") THEN
    !
    zb(jj,iwl_layer) = zb(jj,iwl_layer) + zimpl * zmtc_o_d_wl_in(jj)
    !
    zy(jj,iwl_layer) = zy(jj,iwl_layer) &
                        + zmtc_o_d_wl_in(jj) * b%XTI_BLD(jj) &
                        - zexpl * zmtc_o_d_wl_in(jj) * pt_wl(jj,iwl_layer)
    !
  ELSEIF (top%CBEM=="BEM") THEN
    !
    zf_wl_wl(jj) = 1. - b%XF_WALL_MASS(jj) - b%XF_WALL_WIN(jj) - 2.*b%XF_WALL_FLOOR(jj) 
    !
     zb(jj,iwl_layer) = zb(jj,iwl_layer) + zimpl *               &
                        (zchtc_in_wl(jj) * 4./3. + pradht_in(jj) *           &
                      (  b%XF_WALL_MASS(jj) +     b%XF_WALL_WIN  (jj) &
                          + zf_wl_wl(jj) + 2 * b%XF_WALL_FLOOR(jj)))
    !
    zti_rf(jj) = prad_rf_wl(jj) / pradht_in(jj) + pt_wl(jj,iwl_layer)
    !
    zy(jj,iwl_layer) = zy(jj,iwl_layer) +  &
        zchtc_in_wl(jj) * (b%XTI_BLD(jj) - 1./3. * pt_wl(jj, iwl_layer) * (4 * zexpl -1)) + &
        pradht_in(jj) * ( &
           b%XF_WALL_MASS (jj) * (b%XT_MASS(jj,1) - zexpl * pt_wl(jj,iwl_layer)) + &
           b%XF_WALL_WIN  (jj) * (b%XT_WIN2 (jj)  - zexpl * pt_wl(jj,iwl_layer)) + &
           b%XF_WALL_FLOOR(jj) * (pts_fl(jj)      - zexpl * pt_wl(jj,iwl_layer)) + &
           zf_wl_wl(jj) * (pti_wl_b(jj)          - zexpl * pt_wl(jj,iwl_layer)) + &
           b%XF_WALL_FLOOR(jj) * (zti_rf(jj)      - zexpl * pt_wl(jj,iwl_layer)) )+ &
           pload_in_wl(jj)
    !
  ENDIF
  !
END DO
!
!-------------------------------------------------------------------------------
!
!*      5.    heat conduction calculation
!             ---------------------------
!
 CALL layer_e_budget( pt_wl, ptstep, zimpl, t%XHC_WALL, t%XTC_WALL, t%XD_WALL, &
                     za, zb, zc, zy, pdqs_wl )
!
!-------------------------------------------------------------------------------
!
!*   6.   diagnostics of flux echanged with the wall
!         ------------------------------------------
!
!
!* radiative surface temperature used during the energy balance
zts_wl(:) = zimpl * pt_wl(:,1) + zexpl * zts_wl(:)
!
pabs_lw_wl(:) = plw_s_to_w(:) * (zt_sky(:) - zts_wl(:)) + &
     zdf_rd(:) *plw_r_to_w(:) * (pts_rd(:) - zts_wl(:)) + &
                plw_wa_to_wb(:) * (pts_wl_b(:) - zts_wl(:)) + &
                plw_win_to_w(:) * (b%XT_WIN1(:) - zts_wl(:)) + &
     pdn_rd(:) *plw_nr_to_w(:) * (ptsnow_rd(:) - zts_wl(:))
!
IF (SIZE(pts_gd)>0) THEN
  pabs_lw_wl(:) = pabs_lw_wl(:) + plw_g_to_w(:) * (pts_gd(:) - zts_wl(:))
ENDIF
!
!* emitted lw flux
pemit_lw_wl(:) = xstefan * pt_wl(:,1)**4 + &
                   (1 - t%XEMIS_WALL(:))/t%XEMIS_WALL(:) * pabs_lw_wl(:)
!
!* sensible heat flux to outdoor
ph_wl(:) = zrho_acf_w(:) * xcpd/pexns(:) *  &
             ( zimpl*pt_wl(:,1) + zexpl*zts_wl(:) - pt_canyon(:) )
!
IF (top%CBEM=='BEM') THEN
    !
    !compute ZTI_WALL used in flux calculation
    zti_wl_conv(:) = 4./3. * zimpl * pt_wl(:,iwl_layer) + 1./3. * zti_wl(:) * (4 * zexpl -1.)
    zti_wl(:) = zexpl * zti_wl(:) + zimpl * pt_wl(:,iwl_layer) 
    !
    !compute IR exchanged fluxes with the roof in the wall balance
    zrad_wl_rf(:)  = pradht_in(:)     * (zti_wl(:) - zti_rf(:))
    !compute the difference with the IR flux echanged between roof and wall for
    !the roof balance, this difference will be added to the floor 
    zdif_rad_wl_rf(:)  = zrad_wl_rf(:) + prad_rf_wl(:)

    !compute exchanged fluxes with other surfaces for which the balance is done after
    prad_wl_fl(:) = pradht_in(:)     * (zti_wl(:) - pts_fl(:))
    ! contribution of the difference between the flux from the wall/flux from
    ! the roof
    prad_wl_fl(:)   = prad_wl_fl(:) + zdif_rad_wl_rf(:) 
    prad_wl_ma(:)   = pradht_in(:)     * (zti_wl(:) - b%XT_MASS(:,1))
    prad_wl_win(:)  = pradht_in(:)     * (zti_wl(:) - b%XT_WIN2(:))
    pconv_wl_bld(:) = zchtc_in_wl(:) * (zti_wl_conv(:) - b%XTI_BLD (:))
    !
    !
    pflx_bld_wl(:) = -pconv_wl_bld(:) &
       + b%XF_WALL_FLOOR(:) * prad_rf_wl(:)                             &
       + pradht_in(:) *(b%XF_WALL_MASS (:) * (b%XT_MASS(:,1) - zti_wl(:)) &
                      + b%XF_WALL_WIN  (:) * (b%XT_WIN2 (:)  - zti_wl(:)) &
                      + b%XF_WALL_FLOOR(:) * (pts_fl(:)      - zti_wl(:)) )
!
!* wall energy residual imbalance for verification
!  -----------------------------------------------
  pimb_wl(:) = pabs_sw_wl(:) + pabs_lw_wl(:)  &
               - pdqs_wl(:) - ph_wl(:)         &
               + pflx_bld_wl(:)
!
!
END IF
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('WALL_LAYER_E_BUDGET',1,zhook_handle)
!
END SUBROUTINE wall_layer_e_budget