SURFEX v8.1
General documentation of Surfex
roof_layer_e_budget.F90
Go to the documentation of this file.
1 !SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
3 !SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SFX_LIC for details. version 1.
5 ! #########
6  SUBROUTINE roof_layer_e_budget(TOP, T, B, PQSAT_ROOF, PAC_BLD, PTSTEP, PDN_ROOF, &
7  PRHOA, PAC_ROOF, PAC_ROOF_WAT, PLW_RAD, PPS, &
8  PDELT_ROOF, PTA, PQA, PEXNA, PEXNS, PABS_SW_ROOF, &
9  PGSNOW_ROOF, PFLX_BLD_ROOF, PDQS_ROOF, PABS_LW_ROOF,&
10  PHFREE_ROOF, PLEFREE_ROOF, PIMB_ROOF, &
11  PG_GREENROOF_ROOF, PRADHT_IN, PTS_FLOOR, PTI_WALL, &
12  PRAD_ROOF_WALL, PRAD_ROOF_WIN, PRAD_ROOF_FLOOR, &
13  PRAD_ROOF_MASS, PCONV_ROOF_BLD, PRR, PLOAD_IN_ROOF )
14 ! ##########################################################################
15 !
16 !!**** *ROOF_LAYER_E_BUDGET*
17 !!
18 !! PURPOSE
19 !! -------
20 !
21 ! Computes the evoultion of surface temperature of roofs
22 !
23 !
24 !!** METHOD
25 ! ------
26 !
27 !
28 !
29 !
30 ! 5 : equation for evolution of Ts_roof
31 ! *********************************
32 !
33 ! dTt_1(t) / dt = 1/(dt_1*Ct_1) * ( Rn - H - LE
34 ! - 2*Kt_1*(Tt_1-Tt_2)/(dt_1 +dt_2) )
35 !
36 ! dTt_k(t) / dt = 1/(dt_k*Ct_k) * (- 2*Kt_k-1*(Tt_k-Tt_k-1)/(dt_k-1 +dt_k)
37 ! - 2*Kt_k *(Tt_k-Tt_k+1)/(dt_k+1 +dt_k) )
38 !
39 ! with
40 !
41 ! K*_k = (d*_k+ d*_k+1)/(d*_k/k*_k+ d*_k+1/k*_k+1)
42 !
43 ! Rn = (dir_Rg + sca_Rg) (1-a) + emis * ( Rat - sigma Ts**4 (t+dt) )
44 !
45 ! H = rho Cp CH V ( Ts (t+dt) - Tas )
46 !
47 ! LE = rho Lv CH V ( qs (t+dt) - qas )
48 !
49 ! where the as subscript denotes atmospheric values at ground level
50 ! (and not at first half level)
51 !
52 ! The tridiagonal systel is linearized with
53 !
54 ! using Ts**4 (t+dt) = Ts**4 (t) + 4*Ts**3 (t) * ( Ts(t+dt) - Ts(t) )
55 !
56 ! and qs (t+dt) = Hu(t) * qsat(t) + Hu(t) dqsat/dT * ( Ts(t+dt) - Ts(t) )
57 !
58 !
59 !
60 !! EXTERNAL
61 !! --------
62 !!
63 !!
64 !! IMPLICIT ARGUMENTS
65 !! ------------------
66 !!
67 !! MODD_CST
68 !!
69 !!
70 !! REFERENCE
71 !! ---------
72 !!
73 !!
74 !! AUTHOR
75 !! ------
76 !!
77 !! V. Masson * Meteo-France *
78 !!
79 !! MODIFICATIONS
80 !! -------------
81 !! Original 23/01/98
82 !! 17/10/05 (G. Pigeon) computation of storage inside the roofs
83 !! 26/04/12 (G. Pigeon) add term of heating of rain (new arg PRR+XCL)
84 !! 09/12 (G. Pigeon) modif of indoor conv. coef and implicit calculation
85 !! 10/12 (G. Pigeon) add indoor solar heat load
86 !-------------------------------------------------------------------------------
87 !
88 !* 0. DECLARATIONS
89 ! ------------
90 !
92 USE modd_teb_n, ONLY : teb_t
93 USE modd_bem_n, ONLY : bem_t
94 !
95 USE modd_surf_par, ONLY : xundef
96 USE modd_csts,ONLY : xcpd, xlvtt, xstefan, xcl
97 !
98 USE mode_thermos
99 !
100 USE modi_layer_e_budget
101 USE modi_layer_e_budget_get_coef
102 USE mode_conv_doe
103 !
104 USE yomhook ,ONLY : lhook, dr_hook
105 USE parkind1 ,ONLY : jprb
106 !
107 IMPLICIT NONE
108 !
109 !* 0.1 declarations of arguments
110 !
111 TYPE(teb_options_t), INTENT(INOUT) :: TOP
112 TYPE(teb_t), INTENT(INOUT) :: T
113 TYPE(bem_t), INTENT(INOUT) :: B
114 !
115 REAL, DIMENSION(:), INTENT(INOUT) :: PQSAT_ROOF ! q_sat(Ts)
116 REAL, DIMENSION(:), INTENT(IN) :: PAC_BLD ! aerodynamical resistance
117  ! inside building itself
118 REAL, INTENT(IN) :: PTSTEP ! time step
119 REAL, DIMENSION(:), INTENT(IN) :: PDN_ROOF ! roof snow fraction
120 REAL, DIMENSION(:), INTENT(IN) :: PRHOA ! air density
121 REAL, DIMENSION(:), INTENT(IN) :: PAC_ROOF ! aerodynamical conductance
122 REAL, DIMENSION(:), INTENT(IN) :: PAC_ROOF_WAT ! aerodynamical conductance (for water)
123 REAL, DIMENSION(:), INTENT(IN) :: PLW_RAD ! atmospheric infrared radiation
124 REAL, DIMENSION(:), INTENT(IN) :: PPS ! pressure at the surface
125 REAL, DIMENSION(:), INTENT(IN) :: PDELT_ROOF ! fraction of water
126 REAL, DIMENSION(:), INTENT(IN) :: PTA ! air temperature at roof level
127 REAL, DIMENSION(:), INTENT(IN) :: PQA ! air specific humidity
128  ! at roof level
129 REAL, DIMENSION(:), INTENT(IN) :: PEXNA ! exner function
130 REAL, DIMENSION(:), INTENT(IN) :: PEXNS ! surface exner function
131 REAL, DIMENSION(:), INTENT(IN) :: PABS_SW_ROOF ! absorbed solar radiation
132 REAL, DIMENSION(:), INTENT(IN) :: PGSNOW_ROOF ! roof snow conduction
133 ! ! heat fluxes at mantel
134 ! ! base
135 REAL, DIMENSION(:), INTENT(IN) :: PG_GREENROOF_ROOF ! heat conduction flux
136 ! between greenroof and
137 ! structural roof
138 REAL, DIMENSION(:), INTENT(OUT) :: PFLX_BLD_ROOF ! flux from bld to roof
139 REAL, DIMENSION(:), INTENT(OUT) :: PDQS_ROOF ! heat storage inside the roofs
140 REAL, DIMENSION(:), INTENT(OUT) :: PABS_LW_ROOF ! absorbed infra-red rad.
141 REAL, DIMENSION(:), INTENT(OUT) :: PHFREE_ROOF ! sensible heat flux of the
142  ! snow free part of the roof
143 REAL, DIMENSION(:), INTENT(OUT) :: PLEFREE_ROOF ! latent heat flux of the
144  ! snow free part of the roof
145 REAL, DIMENSION(:), INTENT(OUT) :: PIMB_ROOF ! residual energy imbalance
146  ! of the roof for
147 REAL, DIMENSION(:), INTENT(IN) :: PRADHT_IN ! Indoor radiant heat transfer coefficient
148  ! [W K-1 m-2]
149 REAL, DIMENSION(:), INTENT(IN) :: PTS_FLOOR ! surf. floor temp. (contact with bld air)
150 REAL, DIMENSION(:), INTENT(IN) :: PTI_WALL ! indoor wall temp.
151 REAL, DIMENSION(:), INTENT(OUT) :: PRAD_ROOF_WALL ! rad. fluxes from roof to wall [W m-2(roof)]
152 REAL, DIMENSION(:), INTENT(OUT) :: PRAD_ROOF_WIN ! rad. fluxes from roof to win [W m-2(roof)]
153 REAL, DIMENSION(:), INTENT(OUT) :: PRAD_ROOF_FLOOR! rad. fluxes from roof to floor [W m-2(roof)]
154 REAL, DIMENSION(:), INTENT(OUT) :: PRAD_ROOF_MASS ! rad. fluxes from roof to mass [W m-2(roof)]
155 REAL, DIMENSION(:), INTENT(OUT) :: PCONV_ROOF_BLD ! conv. fluxes from roof to bld [W m-2(roof)]
156 REAL, DIMENSION(:), INTENT(IN) :: PRR ! rain rate [kg m-2 s-1]
157 REAL, DIMENSION(:), INTENT(IN) :: PLOAD_IN_ROOF ! solar + int heat gain on roof W/m2 [roof]
158 !
159 !* 0.2 declarations of local variables
160 !
161 REAL :: ZIMPL = 1.0 ! implicit coefficient
162 REAL :: ZEXPL = 0.0 ! explicit coefficient
163 !
164 REAL, DIMENSION(SIZE(PTA)) :: ZDF_ROOF ! snow-free fraction
165 REAL, DIMENSION(SIZE(PTA),SIZE(T%XT_ROOF,2)) :: ZA,& ! lower diag.
166  ZB,& ! main diag.
167  ZC,& ! upper diag.
168  ZY ! r.h.s.
169 !
170 REAL, DIMENSION(SIZE(PTA)) :: ZDQSAT_ROOF ! dq_sat/dTs
171 REAL, DIMENSION(SIZE(PTA)) :: ZRHO_ACF_ROOF ! conductance * rho
172 REAL, DIMENSION(SIZE(PTA)) :: ZRHO_ACF_ROOF_WAT! conductance * rho (for water)
173 REAL, DIMENSION(SIZE(PTA)) :: ZMTC_O_D_ROOF_IN ! thermal capacity times layer depth
174 REAL, DIMENSION(SIZE(PTA)) :: ZTS_ROOF ! roof surface temperature at previous time step
175 REAL, DIMENSION(SIZE(PTA)) :: ZTRAD_ROOF ! roof radiative surface temperature at intermediate time step
176 REAL, DIMENSION(SIZE(PTA)) :: ZTAER_ROOF ! roof aerodyn. surface temperature at intermediate time step
177 REAL, DIMENSION(SIZE(PTA)) :: ZHEAT_RR ! heat used too cool/heat the rain from the roof
178 REAL, DIMENSION(SIZE(PTA)) :: ZTI_ROOF ! temperature of internal roof layer used for radiative exchanges
179 REAL, DIMENSION(SIZE(PTA)) :: ZTI_ROOF_CONV ! temperature of internal roof layer used for convective exchanges
180 REAL, DIMENSION(SIZE(PTA)) :: ZCHTC_IN_ROOF ! Indoor roof convec heat transfer coefficient
181  ! [W K-1 m-2(bld)]
182 !
183 INTEGER :: JJ
184 INTEGER :: IROOF_LAYER ! number of roof layers
185 INTEGER :: JLAYER ! loop counter
186 REAL(KIND=JPRB) :: ZHOOK_HANDLE
187 !-------------------------------------------------------------------------------
188 !
189 IF (lhook) CALL dr_hook('ROOF_LAYER_E_BUDGET',0,zhook_handle)
190 !
191 prad_roof_wall(:) = xundef
192 prad_roof_win(:) = xundef
193 prad_roof_floor(:)= xundef
194 prad_roof_mass(:) = xundef
195 pconv_roof_bld(:) = xundef
196 !
197 ! *Convection heat transfer coefficients [W m-2 K-1] from EP Engineering Reference
198 !
199 iroof_layer = SIZE(t%XT_ROOF,2)
200 !
201 zchtc_in_roof(:) = chtc_down_doe(t%XT_ROOF(:,iroof_layer), b%XTI_BLD(:))
202 DO jj=1,SIZE(zchtc_in_roof)
203  zchtc_in_roof(jj) = max(1., zchtc_in_roof(jj))
204 ENDDO
205 !
206  CALL layer_e_budget_get_coef(t%XT_ROOF, ptstep, zimpl, t%XHC_ROOF, t%XTC_ROOF, t%XD_ROOF, &
207  za, zb, zc, zy )
208 !
209 !
210 DO jj=1,SIZE(pdn_roof)
211  !
212  zdf_roof(jj) = 1. - pdn_roof(jj)
213  !
214  zts_roof(jj) = t%XT_ROOF(jj,1)
215  zti_roof(jj) = t%XT_ROOF(jj,iroof_layer)
216  !
217  !* 2. Roof Ts coefficients
218  ! --------------------
219  !
220  zrho_acf_roof(jj) = prhoa(jj) * pac_roof(jj)
221  zrho_acf_roof_wat(jj) = prhoa(jj) * pac_roof_wat(jj)
222  !
223  IF (top%CBEM .EQ. 'DEF') THEN
224  zmtc_o_d_roof_in(jj) = 2. * t%XTC_ROOF(jj,iroof_layer) / t%XD_ROOF (jj,iroof_layer)
225  zmtc_o_d_roof_in(jj) = 1./( 1./zmtc_o_d_roof_in(jj) + 1./(xcpd*prhoa(jj)*pac_bld(jj)) )
226  ENDIF
227  !
228 ENDDO
229 !
230 !* 2.1 dqsat/dTs, and humidity for roofs
231 ! ---------------------------------
232 !
233 zdqsat_roof(:) = dqsat(zts_roof(:),pps(:),pqsat_roof(:))
234 !
235 !* 2.2 coefficients
236 ! ------------
237 !
238 DO jj=1,SIZE(t%XT_ROOF,1)
239  !
240  zb(jj,1) = zb(jj,1) + zdf_roof(jj) * (1.-t%XGREENROOF(jj)) * ( &
241  zimpl * ( xcpd/pexns(jj) * zrho_acf_roof(jj) &
242  + xlvtt * zrho_acf_roof_wat(jj) * pdelt_roof(jj) * zdqsat_roof(jj) &
243  + xstefan * t%XEMIS_ROOF(jj) * 4.*zts_roof(jj)**3 &
244  + prr(jj) * xcl )) !! heating/cooling of rain
245  !
246  zy(jj,1) = zy(jj,1) + (1.-t%XGREENROOF(jj)) &
247  * (pdn_roof(jj)*pgsnow_roof(jj) + zdf_roof(jj) * ( pabs_sw_roof(jj) &
248  + xcpd * zrho_acf_roof(jj) * ( pta(jj)/pexna(jj) - zexpl*zts_roof(jj)/pexns(jj)) &
249  + t%XEMIS_ROOF(jj)*plw_rad(jj) &
250  + xlvtt * zrho_acf_roof_wat(jj) * pdelt_roof(jj) &
251  * ( pqa(jj) - pqsat_roof(jj) + zimpl * zdqsat_roof(jj) * zts_roof(jj) ) &
252  + xstefan * t%XEMIS_ROOF(jj) * zts_roof(jj)**4 * ( 3.*zimpl-zexpl ) &
253  + prr(jj) * xcl * (pta(jj) - zexpl * zts_roof(jj)) ) ) & !! heating/cooling of rain
254  + t%XGREENROOF(jj)*pg_greenroof_roof(jj)
255  !
256  IF (top%CBEM=="DEF") THEN
257  !
258  zb(jj,iroof_layer) = zb(jj,iroof_layer) + zimpl * zmtc_o_d_roof_in(jj)
259  !
260  zy(jj,iroof_layer) = zy(jj,iroof_layer) &
261  + zmtc_o_d_roof_in(jj) * b%XTI_BLD(jj) &
262  - zexpl * zmtc_o_d_roof_in(jj) * t%XT_ROOF(jj,iroof_layer)
263  !
264  ELSEIF (top%CBEM=="BEM") THEN
265  !
266  zb(jj, iroof_layer) = zb(jj,iroof_layer) + zimpl * &
267  (zchtc_in_roof(jj) * 4./3. + pradht_in(jj) * &
268  (b%XF_FLOOR_MASS(jj) + b%XF_FLOOR_WIN(jj) + &
269  b%XF_FLOOR_WALL(jj) + b%XF_FLOOR_ROOF(jj) ))
270 
271  zy(jj,iroof_layer) = zy(jj,iroof_layer) + &
272  zchtc_in_roof(jj) * (b%XTI_BLD(jj) - 1./3. * t%XT_ROOF(jj, iroof_layer)*(4*zexpl - 1.)) + &
273  pradht_in(jj) * ( &
274  b%XF_FLOOR_MASS (jj) * (b%XT_MASS(jj,1) - zexpl * t%XT_ROOF(jj,iroof_layer)) + &
275  b%XF_FLOOR_WIN (jj) * (b%XT_WIN2 (jj) - zexpl * t%XT_ROOF(jj,iroof_layer)) + &
276  b%XF_FLOOR_WALL (jj) * (pti_wall(jj) - zexpl * t%XT_ROOF(jj,iroof_layer)) + &
277  b%XF_FLOOR_ROOF (jj) * (pts_floor(jj) - zexpl * t%XT_ROOF(jj,iroof_layer))) + &
278  pload_in_roof(jj)
279  !
280  ENDIF
281  !
282 ENDDO
283 !
284 !
285  CALL layer_e_budget( t%XT_ROOF, ptstep, zimpl, t%XHC_ROOF, t%XTC_ROOF, t%XD_ROOF, &
286  za, zb, zc, zy, pdqs_roof )
287 !
288 !-------------------------------------------------------------------------------
289 !
290 !* diagnostic: computation of flux between bld and internal roof layer
291 DO jj=1,SIZE(t%XT_ROOF,1)
292  !
293  zti_roof_conv(jj) = 4./3. * zimpl * t%XT_ROOF(jj, iroof_layer) + 1./3. * zti_roof(jj) * (4*zexpl -1.)
294  zti_roof(jj) = zexpl * zti_roof(jj) + zimpl * t%XT_ROOF(jj, iroof_layer)
295  SELECT CASE(top%CBEM)
296  CASE("DEF")
297  pflx_bld_roof(jj) = zmtc_o_d_roof_in(jj) * (b%XTI_BLD(jj) - zti_roof(jj))
298  CASE("BEM")
299  prad_roof_wall(jj) = pradht_in(jj) * (zti_roof(jj) - pti_wall(jj))
300  prad_roof_win(jj) = pradht_in(jj) * (zti_roof(jj) - b%XT_WIN2(jj))
301  prad_roof_floor(jj)= pradht_in(jj) * (zti_roof(jj) - pts_floor(jj))
302  prad_roof_mass(jj) = pradht_in(jj) * (zti_roof(jj) - b%XT_MASS(jj,1))
303  pconv_roof_bld(jj) = zchtc_in_roof(jj) * (zti_roof_conv(jj) - b%XTI_BLD (jj))
304  pflx_bld_roof(jj) = -(prad_roof_wall(jj) + prad_roof_win(jj) + prad_roof_floor(jj) + &
305  prad_roof_mass(jj) + pconv_roof_bld(jj))
306  ENDSELECT
307 
308  !
309  !* 8. Infra-red radiation absorbed by roofs
310  ! -------------------------------------
311  !
312  !* radiative surface temperature at intermediate time step
313  ztrad_roof(jj) = ( zts_roof(jj)**4 + &
314  4.*zimpl*zts_roof(jj)**3 * (t%XT_ROOF(jj,1) - zts_roof(jj)) )**0.25
315  !
316  !* absorbed LW
317  pabs_lw_roof(jj) = t%XEMIS_ROOF(jj) * (plw_rad(jj) - xstefan * ztrad_roof(jj)** 4)
318  !
319  !* 9. Sensible heat flux between snow free roof and air
320  ! -------------------------------------------------
321  !
322  !* aerodynamic surface temperature at the intermediate time step
323  ztaer_roof(jj) = zexpl * zts_roof(jj) + zimpl * t%XT_ROOF(jj,1)
324  phfree_roof(jj) = zrho_acf_roof(jj) * xcpd * &
325  ( ztaer_roof(jj)/pexns(jj) - pta(jj)/pexna(jj) )
326  !
327  zheat_rr(jj) = prr(jj) * xcl * (ztaer_roof(jj) - pta(jj))
328  !
329  !* 10. Latent heat flux between snow free roof and air
330  ! -------------------------------------------------
331  !
332  plefree_roof(jj) = zrho_acf_roof_wat(jj) * xlvtt * pdelt_roof(jj) * &
333  ( pqsat_roof(jj) - pqa(jj) + &
334  zimpl * zdqsat_roof(jj) * (t%XT_ROOF(jj,1) - zts_roof(jj)) )
335  !
336  ! 13. Energy imbalance for verification
337  ! ---------------------------------
338  pimb_roof(jj) = pabs_sw_roof(jj) + pabs_lw_roof(jj) - pdqs_roof(jj) &
339  - zdf_roof(jj) * ( phfree_roof(jj) + plefree_roof(jj)) &
340  - pdn_roof(jj) * pgsnow_roof(jj) + pflx_bld_roof(jj)
341  !
342 ENDDO
343 !
344 !* 11. New saturated specified humidity near the roof surface
345 ! ------------------------------------------------------
346 !
347 pqsat_roof(:) = qsat(t%XT_ROOF(:,1),pps(:))
348 !
349 !-------------------------------------------------------------------------
350 IF (lhook) CALL dr_hook('ROOF_LAYER_E_BUDGET',1,zhook_handle)
351 !-------------------------------------------------------------------------
352 !
353 END SUBROUTINE roof_layer_e_budget
real function, dimension(size(pts)) chtc_down_doe(PTS, PTA)
subroutine roof_layer_e_budget(TOP, T, B, PQSAT_ROOF, PAC_BLD, PTSTEP, PDN_ROOF, PRHOA, PAC_ROOF, PAC_ROOF_WAT, PLW_RAD, PPS, PDELT_ROOF, PTA, PQA, PEXNA, PEXNS, PABS_SW_ROOF, PGSNOW_ROOF, PFLX_BLD_ROOF, PDQS_ROOF, PABS_LW_ROOF, PHFREE_ROOF, PLEFREE_ROOF, PIMB_ROOF, PG_GREENROOF_ROOF, PRADHT_IN, PTS_FLOOR, PTI_WALL, PRAD_ROOF_WALL, PRAD_ROOF_WIN, PRAD_ROOF_FLOOR, PRAD_ROOF_MASS, PCONV_ROOF_BLD, PRR, PLOAD_IN_ROOF)
real, save xcpd
Definition: modd_csts.F90:63
real, save xstefan
Definition: modd_csts.F90:59
real, save xlvtt
Definition: modd_csts.F90:70
real, parameter xundef
integer, parameter jprb
Definition: parkind1.F90:32
subroutine layer_e_budget(PT, PTSTEP, PIMPL, PHC, PTC, PD, PA, PB, PC, PY, PDQS)
real, save xcl
Definition: modd_csts.F90:65
logical lhook
Definition: yomhook.F90:15
subroutine layer_e_budget_get_coef(PT, PTSTEP, PIMPL, PHC, PTC, PD, PA, PB, PC, PY)