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(HBEM, PT_WALL, PTS_WALL_B, PTI_WALL_B, PTSTEP,          &
                                   phc_wall, ptc_wall, pd_wall, pdn_road, prhoa, pac_wall, &
                                   pac_bld, pti_bld, plw_rad, pps, pexns,     &
                                   pabs_sw_wall, pt_canyon, pts_road, ptsnow_road,         &
                                   pts_garden, pts_mass, pts_floor, pemis_wall,            &
                                   plw_wa_to_wb, plw_r_to_w, plw_g_to_w, plw_nr_to_w,      &
                                   plw_win_to_w, plw_s_to_w,                               &
                                   pt_win1,                                                &
                                   pflx_bld_wall, pdqs_wall, pt_win2, pabs_lw_wall,        &
                                   pemit_lw_wall, ph_wall, pimb_wall, pf_wall_mass,        &
                                   pf_wall_floor, pf_wall_win, pradht_in, prad_roof_wall,  &
                                   prad_wall_win, prad_wall_floor,                         &
                                   prad_wall_mass, pconv_wall_bld, pload_in_wall ) 
!   ##########################################################################
!
!!****  *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_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
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: pt_wall    ! wall layers temperatures
REAL, DIMENSION(:),   INTENT(IN)  :: pts_wall_b ! opposite wall surface temperature
REAL, DIMENSION(:),   INTENT(IN)  :: pti_wall_b ! opposite wall internal temperature
REAL,               INTENT(IN)    :: ptstep       ! time step
 CHARACTER(LEN=3), INTENT(IN)      :: hbem         ! Building Energy model 'DEF' or 'BEM'
REAL, DIMENSION(:,:), INTENT(IN)  :: phc_wall     ! heat capacity for wall layers
REAL, DIMENSION(:,:), INTENT(IN)  :: ptc_wall     ! thermal conductivity for wall layers
REAL, DIMENSION(:,:), INTENT(IN)  :: pd_wall      ! depth of wall layers
REAL, DIMENSION(:), INTENT(IN)    :: pdn_road     ! snow-covered fraction on roads
REAL, DIMENSION(:), INTENT(IN)    :: prhoa        ! rho
REAL, DIMENSION(:), INTENT(IN)    :: pac_wall     ! aerodynamical conductance [m/s]
!                                                 ! between wall and canyon
REAL, DIMENSION(:), INTENT(IN)    :: pac_bld      ! aerodynamical conductance
                                                  ! inside the building itself
REAL, DIMENSION(:), INTENT(IN)    :: pti_bld      ! inside building temperature
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_wall ! absorbed solar radiation
REAL, DIMENSION(:), INTENT(IN)    :: pt_canyon    ! air canyon temperature
REAL, DIMENSION(:), INTENT(IN)    :: pts_road     ! road surface temperature
REAL, DIMENSION(:), INTENT(IN)    :: ptsnow_road  ! road snow temperature
REAL, DIMENSION(:), INTENT(IN)    :: pts_garden   ! green area surface temperature
REAL, DIMENSION(:), INTENT(IN)    :: pts_mass     ! surface mass temperature  [K]
REAL, DIMENSION(:), INTENT(IN)    :: pts_floor    ! floor layers temperatures [K]
REAL, DIMENSION(:), INTENT(IN)    :: pemis_wall   ! wall emissivity
!
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(IN)    :: pt_win1      ! outdoor window temperature [K]
REAL, DIMENSION(:), INTENT(OUT)   :: pflx_bld_wall! flux from bld to wall
REAL, DIMENSION(:), INTENT(INOUT) :: pdqs_wall    ! heat storage inside the wall 
REAL, DIMENSION(:), INTENT(IN)    :: pt_win2      ! indoor window temperature [K]
REAL, DIMENSION(:), INTENT(OUT)   :: pabs_lw_wall ! absorbed infrared rad. [W m-2(wall)]
REAL, DIMENSION(:), INTENT(OUT)   :: pemit_lw_wall  ! LW flux emitted by the wall [W m-2(wall)]
REAL, DIMENSION(:), INTENT(OUT)   :: ph_wall      ! Sensible heat flux from wall to air [W/m2(wall)]
                                                  ! wall = facade - glazing
REAL, DIMENSION(:), INTENT(OUT)   :: pimb_wall    ! wall residual energy imbalance 
                                                  ! for verification [W m-2]
REAL, DIMENSION(:),   INTENT(IN)  :: pf_wall_floor ! View factor wall-floor
REAL, DIMENSION(:),   INTENT(IN)  :: pf_wall_mass  ! View factor wall-mass
REAL, DIMENSION(:),   INTENT(IN)  :: pf_wall_win   ! View factor wall-win
REAL, DIMENSION(:),   INTENT(IN)  :: pradht_in     ! Indoor radiant heat transfer coefficient
                                                    ! [W K-1 m-2]
REAL, DIMENSION(:), INTENT(IN)    :: prad_roof_wall ! rad. fluxes from roof to wall [W m-2(roof)]
REAL, DIMENSION(:), INTENT(OUT)   :: prad_wall_win  ! rad. fluxes from wall to win  [W m-2(wall)]
REAL, DIMENSION(:), INTENT(OUT)   :: prad_wall_floor! rad. fluxes from wall to floor [W m-2(wall)]
REAL, DIMENSION(:), INTENT(OUT)   :: prad_wall_mass ! rad. fluxes from wall to mass [W m-2(wall)]
REAL, DIMENSION(:), INTENT(OUT)   :: pconv_wall_bld ! conv. fluxes from wall to bld [W m-2(wall)]

REAL, DIMENSION(:), INTENT(IN)    :: pload_in_wall  ! 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_WALL,2)) ::  za,& ! lower diag.
                                               zb,& ! main  diag.
                                               zc,& ! upper diag.
                                               zy   ! r.h.s.                       
!
REAL, DIMENSION(SIZE(PPS)) :: zmtc_o_d_wall_in
REAL, DIMENSION(SIZE(PPS)) :: zdf_road    ! 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_wall       ! wall surface temperature
REAL, DIMENSION(SIZE(PPS)) :: zti_wall       ! wall indoor surface temperature
REAL, DIMENSION(SIZE(PPS)) :: zti_wall_conv  ! wall indoor surface temperature for conv. flux
REAL, DIMENSION(SIZE(PPS)) :: zt_sky         ! sky temperature [K]
!
REAL, DIMENSION(SIZE(PPS)) :: zti_roof       ! Indoor roof temperature [K]
REAL, DIMENSION(SIZE(PPS)) :: zdif_rad_wall_roof !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_wall_roof     ! rad flux between the wall and the roof computed for the wall balance
REAL, DIMENSION(SIZE(PPS)) :: zf_wall_wall       ! View factor wall-wall inside the building
REAL, DIMENSION(SIZE(PPS)) :: zchtc_in_wall  ! indoor convective heat transfer coeff wall [W m-2 K-1]
INTEGER :: iwall_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_wall  = xundef
pimb_wall      = xundef
prad_wall_win  = xundef
prad_wall_floor= xundef
prad_wall_mass = xundef
pconv_wall_bld = xundef

!
 CALL layer_e_budget_get_coef( pt_wall, ptstep, zimpl, phc_wall, ptc_wall, pd_wall, &
                              za, zb, zc, zy )
!
iwall_layer = SIZE(pt_wall,2)
!
DO jj=1,SIZE(pdn_road)
  !
  zdf_road(jj) = 1. - pdn_road(jj)
  !
  zts_wall(jj) = pt_wall(jj,1)
  zti_wall(jj) = pt_wall(jj, iwall_layer)
  !
  !*      2.1    outdoor convective flux properties 
  !              ----------------------------------
  !
  zrho_acf_w(jj) = prhoa(jj) * pac_wall(jj)
  !
  !*      2.2    Sky temperature
  !              ---------------
  !
  zt_sky(jj) = (plw_rad(jj)/xstefan)**0.25
  !  
  !*      2.3    indoor average thermal conductivity
  !              -----------------------------------
  !
  IF (hbem .EQ. "DEF") THEN
    zmtc_o_d_wall_in(jj) = 2. * ptc_wall(jj,iwall_layer) / pd_wall(jj,iwall_layer)
    zmtc_o_d_wall_in(jj) = 1./(  1./zmtc_o_d_wall_in(jj)  + 1./(xcpd*prhoa(jj)*pac_bld(jj)))
  ENDIF
ENDDO  
  !  
  !*      2.4    indoor convective coefficient
  !              -----------------------------
  !
  zchtc_in_wall(:) = chtc_vert_doe(pt_wall(:,iwall_layer), pti_bld(:))
  DO jj=1,SIZE(zchtc_in_wall)
     zchtc_in_wall(jj) = max(1., zchtc_in_wall(jj))
  ENDDO

!
!-------------------------------------------------------------------------------
!
!*      3.    Outer wall layer coefficients
!             ------------------------------
!
DO jj=1,SIZE(pt_wall,1)
  !
  zb(jj,1) = zb(jj,1) + zimpl * xcpd/pexns(jj) * zrho_acf_w(jj)
  !
  zy(jj,1) = zy(jj,1) + pabs_sw_wall(jj)  &
                      + xcpd/pexns(jj) * zrho_acf_w(jj) * ( pt_canyon(jj) - zexpl * zts_wall(jj) )
  !
  !
  zb(jj,1) = zb(jj,1) &
             + zimpl * ( plw_s_to_w(jj) + plw_wa_to_wb(jj)                &
                       + zdf_road(jj)*plw_r_to_w(jj) +  plw_g_to_w(jj)    &
                       + pdn_road(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_wall(jj))  &
                       + plw_wa_to_wb(jj) * (pts_wall_b(jj) - zexpl * zts_wall(jj))  &
       + zdf_road(jj) *  plw_r_to_w(jj)  * (pts_road(jj)   - zexpl * zts_wall(jj))  &
                       + plw_g_to_w(jj)  * (pts_garden(jj) - zexpl * zts_wall(jj)) &
       + pdn_road(jj) *  plw_nr_to_w(jj) * (ptsnow_road(jj)- zexpl * zts_wall(jj))  &
                       + plw_win_to_w(jj) * (pt_win1(jj)    - zexpl * zts_wall(jj))  

  !
ENDDO
!
!-------------------------------------------------------------------------------
!
!*      4.    Inside wall layer coefficients
!             -----------------------------
!
DO jj=1,SIZE(pt_wall,1)
  !                
  IF (hbem=="DEF") THEN
    !
    zb(jj,iwall_layer) = zb(jj,iwall_layer) + zimpl * zmtc_o_d_wall_in(jj)
    !
    zy(jj,iwall_layer) = zy(jj,iwall_layer) &
                        + zmtc_o_d_wall_in(jj) * pti_bld(jj) &
                        - zexpl * zmtc_o_d_wall_in(jj) * pt_wall(jj,iwall_layer)
    !
  ELSEIF (hbem=="BEM") THEN
    !
    zf_wall_wall(jj) = 1. - pf_wall_mass(jj) - pf_wall_win(jj) - 2.*pf_wall_floor(jj) 
    !
     zb(jj,iwall_layer) = zb(jj,iwall_layer) + zimpl *               &
                        (zchtc_in_wall(jj) * 4./3. + pradht_in(jj) *           &
                         (  pf_wall_mass(jj) +     pf_wall_win(jj) &
                          + zf_wall_wall(jj) + 2 * pf_wall_floor(jj)))
    !
    zti_roof(jj) = prad_roof_wall(jj) / pradht_in(jj) + pt_wall(jj,iwall_layer)
    zy(jj,iwall_layer) = zy(jj,iwall_layer) +  &
        zchtc_in_wall(jj) * (pti_bld(jj) - 1./3. * pt_wall(jj, iwall_layer) * (4 * zexpl -1)) + &
        pradht_in(jj) * ( &
           pf_wall_mass(jj) * (pts_mass(jj)  - zexpl * pt_wall(jj,iwall_layer)) + &
           pf_wall_win(jj) * (pt_win2(jj)  - zexpl * pt_wall(jj,iwall_layer)) + &
           pf_wall_floor(jj) * (pts_floor(jj) - zexpl * pt_wall(jj,iwall_layer)) + &
           zf_wall_wall(jj) * (pti_wall_b(jj)- zexpl * pt_wall(jj,iwall_layer)) + &
           pf_wall_floor(jj) * (zti_roof(jj) - zexpl * pt_wall(jj,iwall_layer)) )+ &
           pload_in_wall(jj)
    !
  ENDIF
  !
END DO
!
!-------------------------------------------------------------------------------
!
!*      5.    heat conduction calculation
!             ---------------------------
!
 CALL layer_e_budget( pt_wall, ptstep, zimpl, phc_wall, ptc_wall, pd_wall, &
                     za, zb, zc, zy, pdqs_wall )
!
!-------------------------------------------------------------------------------
!
!*   6.   diagnostics of flux echanged with the wall
!         ------------------------------------------
!
!
!* radiative surface temperature used during the energy balance
zts_wall(:) = zimpl * pt_wall(:,1) + zexpl * zts_wall(:)
!
pabs_lw_wall(:) = plw_s_to_w(:) * (zt_sky(:) - zts_wall(:)) + &
     zdf_road(:) *plw_r_to_w(:) * (pts_road(:) - zts_wall(:)) + &
                  plw_g_to_w(:) * (pts_garden(:) - zts_wall(:)) + &
                  plw_wa_to_wb(:) * (pts_wall_b(:) - zts_wall(:)) + &
                  plw_win_to_w(:) * (pt_win1(:) - zts_wall(:)) + &
     pdn_road(:) *plw_nr_to_w(:) * (ptsnow_road(:) - zts_wall(:))
!
!* emitted lw flux
pemit_lw_wall(:) = xstefan * pt_wall(:,1)**4 + &
                   (1 - pemis_wall(:))/pemis_wall(:) * pabs_lw_wall(:)
!
!* sensible heat flux to outdoor
ph_wall(:) = zrho_acf_w(:) * xcpd/pexns(:) *  &
             ( zimpl*pt_wall(:,1) + zexpl*zts_wall(:) - pt_canyon(:) )
!
IF (hbem=='BEM') THEN
    !
    !compute ZTI_WALL used in flux calculation
    zti_wall_conv(:) = 4./3. * zimpl * pt_wall(:,iwall_layer) + 1./3. * zti_wall(:) * (4 * zexpl -1.)
    zti_wall(:) = zexpl * zti_wall(:) + zimpl * pt_wall(:,iwall_layer) 
    !
    !compute IR exchanged fluxes with the roof in the wall balance
    zrad_wall_roof(:)  = pradht_in(:)     * (zti_wall(:) - zti_roof(:))
    !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_wall_roof(:)  = zrad_wall_roof(:) + prad_roof_wall(:)

    !compute exchanged fluxes with other surfaces for which the balance is done after
    prad_wall_floor(:) = pradht_in(:)     * (zti_wall(:) - pts_floor(:))
    ! contribution of the difference between the flux from the wall/flux from
    ! the roof
    prad_wall_floor(:) = prad_wall_floor(:) + zdif_rad_wall_roof(:) 
    prad_wall_mass(:)  = pradht_in(:)     * (zti_wall(:) - pts_mass(:))
    prad_wall_win(:)   = pradht_in(:)     * (zti_wall(:) - pt_win2(:))
    pconv_wall_bld(:)  = zchtc_in_wall(:) * (zti_wall_conv(:) - pti_bld(:))
    !
    !
    pflx_bld_wall(:) = -pconv_wall_bld(:) &
       + pf_wall_floor(:) * prad_roof_wall(:)                             &
       + pradht_in(:) *(pf_wall_mass(:) * (pts_mass(:) - zti_wall(:)) &
                       + pf_wall_win(:) * (pt_win2(:)  - zti_wall(:)) &
                       + pf_wall_floor(:) * (pts_floor(:) - zti_wall(:)) )
!
!* wall energy residual imbalance for verification
!  -----------------------------------------------
  pimb_wall(:) = pabs_sw_wall(:) + pabs_lw_wall(:)  &
               - pdqs_wall(:) - ph_wall(:)         &
               + pflx_bld_wall(:)
!
!
END IF
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('WALL_LAYER_E_BUDGET',1,zhook_handle)
!
END SUBROUTINE wall_layer_e_budget