roof_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 roof_layer_e_budget(PT_ROOF, PQSAT_ROOF, PTI_BLD, PAC_BLD, PTSTEP, &
                                   hbld, phc_roof, ptc_roof, pd_roof, pdn_roof,   &
                                   prhoa, pac_roof, pac_roof_wat, plw_rad, pps,   &
                                   pdelt_roof, pta, pqa, pexna, pexns,            &
                                   pabs_sw_roof, pgsnow_roof, pemis_roof,         &
                                   pflx_bld_roof, pdqs_roof, pabs_lw_roof,        &
                                   phfree_roof, plefree_roof, pimb_roof,          &
                                   pfrac_gr, pg_greenroof_roof,                   &
                                   pf_floor_mass, pf_floor_wall, pf_floor_win,    &
                                   pf_floor_roof, pradht_in,                      &
                                   pts_mass, pt_win2, pts_floor, pti_wall,        &
                                   prad_roof_wall, prad_roof_win, prad_roof_floor,&
                                   prad_roof_mass, pconv_roof_bld, prr,           &
                                   pload_in_roof                                  )
!   ##########################################################################
!
!!****  *ROOF_LAYER_E_BUDGET*  
!!
!!    PURPOSE
!!    -------
!
!     Computes the evoultion of surface temperature of roofs
!         
!     
!!**  METHOD
!     ------
!
!
!
!
!    5 : equation for evolution of Ts_roof
!        *********************************
!
!     dTt_1(t) / dt = 1/(dt_1*Ct_1) * (  Rn - H - LE
!                                      - 2*Kt_1*(Tt_1-Tt_2)/(dt_1 +dt_2)       )
!
!     dTt_k(t) / dt = 1/(dt_k*Ct_k) * (- 2*Kt_k-1*(Tt_k-Tt_k-1)/(dt_k-1 +dt_k) 
!                                      - 2*Kt_k  *(Tt_k-Tt_k+1)/(dt_k+1 +dt_k) )
!
!       with
!
!       K*_k  = (d*_k+ d*_k+1)/(d*_k/k*_k+ d*_k+1/k*_k+1)
!
!       Rn = (dir_Rg + sca_Rg) (1-a) + emis * ( Rat - sigma Ts**4 (t+dt) )
!
!       H  = rho Cp CH V ( Ts (t+dt) - Tas )
!
!       LE = rho Lv CH V ( qs (t+dt) - qas )
!
!      where the as subscript denotes atmospheric values at ground level
!      (and not at first half level)
!
!      The tridiagonal systel is linearized with
!
!       using      Ts**4 (t+dt) = Ts**4 (t) + 4*Ts**3 (t) * ( Ts(t+dt) - Ts(t) )
!
!       and  qs (t+dt) = Hu(t) * qsat(t) + Hu(t) dqsat/dT * ( Ts(t+dt) - Ts(t) )
!
!
!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    MODD_CST
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      V. Masson           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    23/01/98 
!!                  17/10/05 (G. Pigeon) computation of storage inside the roofs
!!                  26/04/12 (G. Pigeon) add term of heating of rain (new arg PRR+XCL)
!!                     09/12 (G. Pigeon) modif of indoor conv. coef and implicit calculation
!!                     10/12 (G. Pigeon) add indoor solar heat load
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_surf_par,  ONLY : xundef
USE modd_csts,ONLY : xcpd, xlvtt, xstefan, xcl
!
USE mode_thermos
!
USE modi_layer_e_budget
USE modi_layer_e_budget_get_coef
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_roof      ! roof layers temperatures
REAL, DIMENSION(:), INTENT(INOUT) :: pqsat_roof     ! q_sat(Ts)
REAL, DIMENSION(:), INTENT(IN)    :: pti_bld        ! inside building temp.
REAL, DIMENSION(:), INTENT(IN)    :: pac_bld        ! aerodynamical resistance
                                                    ! inside building itself
REAL,               INTENT(IN)    :: ptstep         ! time step
 CHARACTER(LEN=3), INTENT(IN)      :: hbld           ! Building Energy model 'DEF' or 'BEM'  
REAL, DIMENSION(:,:), INTENT(IN)  :: phc_roof       ! heat capacity for roof layers
REAL, DIMENSION(:,:), INTENT(IN)  :: ptc_roof       ! thermal conductivity for roof layers
REAL, DIMENSION(:,:), INTENT(IN)  :: pd_roof        ! depth of roof layers
REAL, DIMENSION(:), INTENT(IN)    :: pdn_roof       ! roof snow fraction
REAL, DIMENSION(:), INTENT(IN)    :: prhoa          ! air density
REAL, DIMENSION(:), INTENT(IN)    :: pac_roof       ! aerodynamical conductance
REAL, DIMENSION(:), INTENT(IN)    :: pac_roof_wat   ! aerodynamical conductance (for water)
REAL, DIMENSION(:), INTENT(IN)    :: plw_rad        ! atmospheric infrared radiation
REAL, DIMENSION(:), INTENT(IN)    :: pps            ! pressure at the surface
REAL, DIMENSION(:), INTENT(IN)    :: pdelt_roof     ! fraction of water
REAL, DIMENSION(:), INTENT(IN)    :: pta            ! air temperature at roof level
REAL, DIMENSION(:), INTENT(IN)    :: pqa            ! air specific humidity
                                                    ! at roof level
REAL, DIMENSION(:), INTENT(IN)    :: pexna          ! exner function
REAL, DIMENSION(:), INTENT(IN)    :: pexns          ! surface exner function
REAL, DIMENSION(:), INTENT(IN)    :: pabs_sw_roof   ! absorbed solar radiation
REAL, DIMENSION(:), INTENT(IN)    :: pgsnow_roof    ! roof snow conduction
!                                                   ! heat fluxes at mantel
!                                                   ! base
REAL, DIMENSION(:), INTENT(IN)    :: pemis_roof     ! roof emissivity
REAL, DIMENSION(:), INTENT(IN)    :: pfrac_gr       ! fraction of green roofs
REAL, DIMENSION(:), INTENT(IN)    :: pg_greenroof_roof ! heat conduction flux
!                                                        between greenroof and
!                                                        structural roof
REAL, DIMENSION(:), INTENT(OUT)   :: pflx_bld_roof  ! flux from bld to roof
REAL, DIMENSION(:), INTENT(OUT)   :: pdqs_roof      ! heat storage inside the roofs
REAL, DIMENSION(:), INTENT(OUT)   :: pabs_lw_roof   ! absorbed infra-red rad.
REAL, DIMENSION(:), INTENT(OUT)   :: phfree_roof    ! sensible heat flux of the
                                                    ! snow free part of the roof
REAL, DIMENSION(:), INTENT(OUT)   :: plefree_roof   ! latent heat flux of the
                                                    ! snow free part of the roof
REAL, DIMENSION(:), INTENT(OUT)   :: pimb_roof      ! residual energy imbalance
                                                    ! of the roof for
REAL, DIMENSION(:), INTENT(IN)    :: pf_floor_mass  ! View factor floor-mass
REAL, DIMENSION(:), INTENT(IN)    :: pf_floor_wall  ! View factor floor-wall
REAL, DIMENSION(:), INTENT(IN)    :: pf_floor_win   ! View factor floor-window
REAL, DIMENSION(:), INTENT(IN)    :: pf_floor_roof  ! View factor floor-roof
REAL, DIMENSION(:), INTENT(IN)    :: pradht_in      ! Indoor radiant heat transfer coefficient
                                                    ! [W K-1 m-2]
REAL, DIMENSION(:), INTENT(IN)    :: pts_mass       ! surf. mass temp. (contact with bld air)
REAL, DIMENSION(:), INTENT(IN)    :: pt_win2        ! indoor wind. temp.
REAL, DIMENSION(:), INTENT(IN)    :: pts_floor      ! surf. floor temp. (contact with bld air)
REAL, DIMENSION(:), INTENT(IN)    :: pti_wall       ! indoor wall temp.
REAL, DIMENSION(:), INTENT(OUT)   :: prad_roof_wall ! rad. fluxes from roof to wall [W m-2(roof)]
REAL, DIMENSION(:), INTENT(OUT)   :: prad_roof_win  ! rad. fluxes from roof to win [W m-2(roof)]
REAL, DIMENSION(:), INTENT(OUT)   :: prad_roof_floor! rad. fluxes from roof to floor [W m-2(roof)]
REAL, DIMENSION(:), INTENT(OUT)   :: prad_roof_mass ! rad. fluxes from roof to mass [W m-2(roof)]
REAL, DIMENSION(:), INTENT(OUT)   :: pconv_roof_bld ! conv. fluxes from roof to bld [W m-2(roof)]
REAL, DIMENSION(:), INTENT(IN)    :: prr ! rain rate [kg m-2 s-1]
REAL, DIMENSION(:), INTENT(IN)    :: pload_in_roof ! solar + int heat gain on roof W/m2 [roof]
!
!*      0.2    declarations of local variables
!
REAL :: zimpl = 1.0        ! implicit coefficient
REAL :: zexpl = 0.0        ! explicit coefficient
!
REAL, DIMENSION(SIZE(PTA)) :: zdf_roof ! snow-free fraction
REAL, DIMENSION(SIZE(PTA),SIZE(PT_ROOF,2)) :: za,& ! lower diag.
                                              zb,& ! main  diag.
                                              zc,& ! upper diag.
                                              zy   ! r.h.s.
!
REAL, DIMENSION(SIZE(PTA)) :: zdqsat_roof      ! dq_sat/dTs
REAL, DIMENSION(SIZE(PTA)) :: zrho_acf_roof    ! conductance * rho
REAL, DIMENSION(SIZE(PTA)) :: zrho_acf_roof_wat! conductance * rho (for water)
REAL, DIMENSION(SIZE(PTA)) :: zmtc_o_d_roof_in ! thermal capacity times layer depth
REAL, DIMENSION(SIZE(PTA)) :: zts_roof         ! roof surface temperature at previous time step
REAL, DIMENSION(SIZE(PTA)) :: ztrad_roof       ! roof radiative surface temperature at intermediate time step
REAL, DIMENSION(SIZE(PTA)) :: ztaer_roof       ! roof aerodyn. surface temperature at intermediate time step
REAL, DIMENSION(SIZE(PTA)) :: zheat_rr         ! heat used too cool/heat the rain from the roof
REAL, DIMENSION(SIZE(PTA)) :: zti_roof         ! temperature of internal roof layer used for radiative exchanges
REAL, DIMENSION(SIZE(PTA)) :: zti_roof_conv    ! temperature of internal roof layer used for convective exchanges
REAL, DIMENSION(SIZE(PTA)) :: zchtc_in_roof      ! Indoor roof convec heat transfer coefficient
                                               ! [W K-1 m-2(bld)]
!
INTEGER :: jj
INTEGER :: iroof_layer           ! number of roof layers
INTEGER :: jlayer                ! loop counter
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('ROOF_LAYER_E_BUDGET',0,zhook_handle)
!
prad_roof_wall(:) = xundef
prad_roof_win(:)  = xundef
prad_roof_floor(:)= xundef
prad_roof_mass(:) = xundef
pconv_roof_bld(:) = xundef
!
! *Convection heat transfer coefficients [W m-2 K-1] from EP Engineering Reference
!
iroof_layer = SIZE(pt_roof,2)
!
zchtc_in_roof(:) = chtc_down_doe(pt_roof(:,iroof_layer), pti_bld(:))
DO jj=1,SIZE(zchtc_in_roof)
   zchtc_in_roof(jj) = max(1., zchtc_in_roof(jj))
ENDDO
!
 CALL layer_e_budget_get_coef( pt_roof, ptstep, zimpl, phc_roof, ptc_roof, pd_roof, &
                              za, zb, zc, zy )
!
!
DO jj=1,SIZE(pdn_roof)
  !
  zdf_roof(jj) = 1. - pdn_roof(jj)
  !
  zts_roof(jj) = pt_roof(jj,1)
  zti_roof(jj) = pt_roof(jj,iroof_layer)
  !
  !*      2.     Roof Ts coefficients
  !              --------------------
  !
  zrho_acf_roof(jj) = prhoa(jj) * pac_roof(jj)
  zrho_acf_roof_wat(jj) = prhoa(jj) * pac_roof_wat(jj)
  !
  IF (hbld .EQ. 'DEF') THEN
    zmtc_o_d_roof_in(jj) = 2. * ptc_roof(jj,iroof_layer) / pd_roof(jj,iroof_layer)
    zmtc_o_d_roof_in(jj) = 1./(  1./zmtc_o_d_roof_in(jj) + 1./(xcpd*prhoa(jj)*pac_bld(jj)) ) 
  ENDIF
  !
ENDDO
!
!*      2.1    dqsat/dTs, and humidity for roofs
!              ---------------------------------
!
zdqsat_roof(:) = dqsat(zts_roof(:),pps(:),pqsat_roof(:))
!
!*      2.2    coefficients
!              ------------
! 
DO jj=1,SIZE(pt_roof,1)
  !
  zb(jj,1) = zb(jj,1) + zdf_roof(jj) * (1.-pfrac_gr(jj)) * (                                       &
                        zimpl * ( xcpd/pexns(jj) * zrho_acf_roof(jj)                               &
                                + xlvtt * zrho_acf_roof_wat(jj) * pdelt_roof(jj) * zdqsat_roof(jj) &
                                + xstefan * pemis_roof(jj) * 4.*zts_roof(jj)**3                    &
                                + prr(jj) * xcl )) !! heating/cooling of rain 
  !
  zy(jj,1) = zy(jj,1) + (1.-pfrac_gr(jj))                                                                 &
                      * (pdn_roof(jj)*pgsnow_roof(jj) + zdf_roof(jj) * ( pabs_sw_roof(jj)                 &
                         + xcpd * zrho_acf_roof(jj) * ( pta(jj)/pexna(jj) - zexpl*zts_roof(jj)/pexns(jj)) &
                         + pemis_roof(jj)*plw_rad(jj)                                                     &                 
                         + xlvtt * zrho_acf_roof_wat(jj) * pdelt_roof(jj)                                 &
                           * ( pqa(jj) - pqsat_roof(jj) + zimpl * zdqsat_roof(jj) * zts_roof(jj) )        &
                         + xstefan * pemis_roof(jj) * zts_roof(jj)**4 * ( 3.*zimpl-zexpl )  &
                         + prr(jj) * xcl * (pta(jj) - zexpl * zts_roof(jj)) ) ) & !! heating/cooling of rain
                      +     pfrac_gr(jj)*pg_greenroof_roof(jj)
  !
  IF (hbld=="DEF") THEN
    !
    zb(jj,iroof_layer) = zb(jj,iroof_layer) + zimpl * zmtc_o_d_roof_in(jj)
    !
    zy(jj,iroof_layer) = zy(jj,iroof_layer) &
                         + zmtc_o_d_roof_in(jj) * pti_bld(jj) &
                         - zexpl * zmtc_o_d_roof_in(jj) * pt_roof(jj,iroof_layer)
    !
  ELSEIF (hbld=="BEM") THEN
    !
    zb(jj, iroof_layer) = zb(jj,iroof_layer) + zimpl * &
                       (zchtc_in_roof(jj) * 4./3. + pradht_in(jj) * &
                        (pf_floor_mass(jj) + pf_floor_win(jj) + &
                         pf_floor_wall(jj) + pf_floor_roof(jj) ))

    zy(jj,iroof_layer) = zy(jj,iroof_layer) + &
       zchtc_in_roof(jj) * (pti_bld(jj) - 1./3. *  pt_roof(jj, iroof_layer)*(4*zexpl - 1.)) + &
       pradht_in(jj) * ( &
          pf_floor_mass(jj) * (pts_mass(jj) - zexpl * pt_roof(jj,iroof_layer)) + &
          pf_floor_win(jj) * (pt_win2(jj) - zexpl * pt_roof(jj,iroof_layer)) + &
          pf_floor_wall(jj) * (pti_wall(jj) - zexpl * pt_roof(jj,iroof_layer)) + &
          pf_floor_roof(jj) * (pts_floor(jj)- zexpl * pt_roof(jj,iroof_layer)) ) + &
          pload_in_roof(jj)
    !
  ENDIF
  !
ENDDO
!
!
 CALL layer_e_budget( pt_roof, ptstep, zimpl, phc_roof, ptc_roof, pd_roof, &
                     za, zb, zc, zy, pdqs_roof )
!
!-------------------------------------------------------------------------------
!
!*     diagnostic: computation of flux between bld and internal roof layer
DO jj=1,SIZE(pt_roof,1)
  !
  zti_roof_conv(jj) = 4./3. * zimpl * pt_roof(jj, iroof_layer) + 1./3. * zti_roof(jj) * (4*zexpl -1.)
  zti_roof(jj) = zexpl * zti_roof(jj) + zimpl * pt_roof(jj, iroof_layer) 
  SELECT CASE(hbld)
  CASE("DEF")
     pflx_bld_roof(jj) = zmtc_o_d_roof_in(jj) * (pti_bld(jj) - zti_roof(jj))
  CASE("BEM")
     prad_roof_wall(jj) = pradht_in(jj)     * (zti_roof(jj) - pti_wall(jj))
     prad_roof_win(jj)  = pradht_in(jj)     * (zti_roof(jj) - pt_win2(jj))
     prad_roof_floor(jj)= pradht_in(jj)     * (zti_roof(jj) - pts_floor(jj))
     prad_roof_mass(jj) = pradht_in(jj)     * (zti_roof(jj) - pts_mass(jj))
     pconv_roof_bld(jj) = zchtc_in_roof(jj) * (zti_roof_conv(jj) - pti_bld(jj))
     pflx_bld_roof(jj) = -(prad_roof_wall(jj) + prad_roof_win(jj) + prad_roof_floor(jj) + prad_roof_mass(jj) + pconv_roof_bld(jj))
  ENDSELECT
  
  !
  !*      8.     Infra-red radiation absorbed by roofs
  !              -------------------------------------
  !
  !* radiative surface temperature at intermediate time step
  ztrad_roof(jj) = ( zts_roof(jj)**4 + &
                   4.*zimpl*zts_roof(jj)**3 * (pt_roof(jj,1) - zts_roof(jj)) )**0.25
  !
  !* absorbed LW
  pabs_lw_roof(jj) = pemis_roof(jj) * (plw_rad(jj) - xstefan * ztrad_roof(jj)** 4)
  !
  !*      9.     Sensible heat flux between snow free roof and air
  !              -------------------------------------------------
  !
  !* aerodynamic surface temperature at the intermediate time step
  ztaer_roof(jj) = zexpl * zts_roof(jj) + zimpl * pt_roof(jj,1)
  phfree_roof(jj) = zrho_acf_roof(jj) * xcpd * &
                   ( ztaer_roof(jj)/pexns(jj) - pta(jj)/pexna(jj) )
  !
  zheat_rr(jj) = prr(jj) * xcl * (ztaer_roof(jj) - pta(jj))
  !
  !*      10.     Latent heat flux between snow free roof and air
  !              -------------------------------------------------
  !
  plefree_roof(jj) = zrho_acf_roof_wat(jj) * xlvtt * pdelt_roof(jj) * &
                     ( pqsat_roof(jj) - pqa(jj) +                     &
                       zimpl * zdqsat_roof(jj) * (pt_roof(jj,1) - zts_roof(jj)) ) 
  !
  !      13.     Energy imbalance for verification
  !              ---------------------------------
  pimb_roof(jj) = pabs_sw_roof(jj) + pabs_lw_roof(jj) - pdqs_roof(jj) &
               - zdf_roof(jj) * ( phfree_roof(jj) + plefree_roof(jj)) &
               - pdn_roof(jj) * pgsnow_roof(jj) + pflx_bld_roof(jj)
  !
ENDDO
!
!*      11.     New saturated specified humidity near the roof surface
!              ------------------------------------------------------
!
pqsat_roof(:) =  qsat(pt_roof(:,1),pps(:))
!
!-------------------------------------------------------------------------
IF (lhook) CALL dr_hook('ROOF_LAYER_E_BUDGET',1,zhook_handle)
!-------------------------------------------------------------------------
!
END SUBROUTINE roof_layer_e_budget