1 | !> @file biometeorology_mod.f90 |
---|
2 | !--------------------------------------------------------------------------------! |
---|
3 | ! This file is part of PALM-4U. |
---|
4 | ! |
---|
5 | ! PALM-4U is free software: you can redistribute it and/or modify it under the |
---|
6 | ! terms of the GNU General Public License as published by the Free Software |
---|
7 | ! Foundation, either version 3 of the License, or (at your option) any later |
---|
8 | ! version. |
---|
9 | ! |
---|
10 | ! PALM-4U is distributed in the hope that it will be useful, but WITHOUT ANY |
---|
11 | ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR |
---|
12 | ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. |
---|
13 | ! |
---|
14 | ! You should have received a copy of the GNU General Public License along with |
---|
15 | ! PALM. If not, see <http://www.gnu.org/licenses/>. |
---|
16 | ! |
---|
17 | ! Copyright 2018-2018 Deutscher Wetterdienst (DWD) |
---|
18 | ! Copyright 2018-2018 Institute of Computer Science, Academy of Sciences, Prague |
---|
19 | ! Copyright 2018-2018 Leibniz Universitaet Hannover |
---|
20 | !--------------------------------------------------------------------------------! |
---|
21 | ! |
---|
22 | ! Current revisions: |
---|
23 | ! ------------------ |
---|
24 | ! |
---|
25 | ! |
---|
26 | ! Former revisions: |
---|
27 | ! ----------------- |
---|
28 | ! $Id: biometeorology_mod.f90 3614 2018-12-10 07:05:46Z schwenkel $ |
---|
29 | ! unused variables removed |
---|
30 | ! |
---|
31 | ! 3593 2018-12-03 13:51:13Z kanani |
---|
32 | ! Bugfix: additional tmrt_grid allocation in case bio_mrt not selected as ouput, |
---|
33 | ! replace degree symbol by degree_C |
---|
34 | ! |
---|
35 | ! 3582 2018-11-29 19:16:36Z suehring |
---|
36 | ! Consistently use bio_fill_value everywhere, |
---|
37 | ! move allocation and initialization of output variables to bio_check_data_output |
---|
38 | ! and bio_3d_data_averaging, |
---|
39 | ! dom_dwd_user, Schrempf: |
---|
40 | ! - integration of UVEM module part from r3474 (edited) |
---|
41 | ! - split UTCI regression into 6 parts |
---|
42 | ! - all data_output_3d is now explicity casted to sp |
---|
43 | ! |
---|
44 | ! 3525 2018-11-14 16:06:14Z kanani |
---|
45 | ! Clean up, renaming from "biom" to "bio", summary of thermal index calculation |
---|
46 | ! into one module (dom_dwd_user) |
---|
47 | ! |
---|
48 | ! 3479 2018-11-01 16:00:30Z gronemeier |
---|
49 | ! - reworked check for output quantities |
---|
50 | ! - assign new palm-error numbers |
---|
51 | ! - set unit according to data standard. |
---|
52 | ! |
---|
53 | ! 3475 2018-10-30 21:16:31Z kanani |
---|
54 | ! Add option for using MRT from RTM instead of simplified global value |
---|
55 | ! |
---|
56 | ! 3464 2018-10-30 18:08:55Z kanani |
---|
57 | ! From branch resler@3462, pavelkrc: |
---|
58 | ! make use of basic_constants_and_equations_mod |
---|
59 | ! |
---|
60 | ! 3448 2018-10-29 18:14:31Z kanani |
---|
61 | ! Initial revision |
---|
62 | ! |
---|
63 | ! |
---|
64 | ! |
---|
65 | ! Authors: |
---|
66 | ! -------- |
---|
67 | ! @author Dominik Froehlich <dominik.froehlich@dwd.de>, thermal indices |
---|
68 | ! @author Jaroslav Resler <resler@cs.cas.cz>, mean radiant temperature |
---|
69 | ! @author Michael Schrempf <schrempf@muk.uni-hannover.de>, uv exposure |
---|
70 | ! |
---|
71 | ! |
---|
72 | ! Description: |
---|
73 | ! ------------ |
---|
74 | !> Biometeorology module consisting of two parts: |
---|
75 | !> 1.: Human thermal comfort module calculating thermal perception of a sample |
---|
76 | !> human being under the current meteorological conditions. |
---|
77 | !> 2.: Calculation of vitamin-D weighted UV exposure |
---|
78 | !> |
---|
79 | !> @todo Alphabetical sorting of "USE ..." lists, "ONLY" list, variable declarations |
---|
80 | !> (per subroutine: first all CHARACTERs, then INTEGERs, LOGICALs, REALs, ) |
---|
81 | !> @todo Comments start with capital letter --> "!-- Include..." |
---|
82 | !> @todo uv_vitd3dose-->new output type necessary (cumulative) |
---|
83 | !> @todo consider upwelling radiation in UV |
---|
84 | !> |
---|
85 | !> @note nothing now |
---|
86 | !> |
---|
87 | !> @bug no known bugs by now |
---|
88 | !------------------------------------------------------------------------------! |
---|
89 | MODULE biometeorology_mod |
---|
90 | |
---|
91 | USE arrays_3d, & |
---|
92 | ONLY: pt, p, u, v, w, q |
---|
93 | |
---|
94 | USE averaging, & |
---|
95 | ONLY: pt_av, q_av, u_av, v_av, w_av |
---|
96 | |
---|
97 | USE basic_constants_and_equations_mod, & |
---|
98 | ONLY: c_p, degc_to_k, l_v, magnus, sigma_sb, pi |
---|
99 | |
---|
100 | USE control_parameters, & |
---|
101 | ONLY: average_count_3d, biometeorology, dz, dz_stretch_factor, & |
---|
102 | dz_stretch_level, humidity, initializing_actions, nz_do3d, & |
---|
103 | simulated_time, surface_pressure |
---|
104 | |
---|
105 | USE date_and_time_mod, & |
---|
106 | ONLY: calc_date_and_time, day_of_year, time_utc |
---|
107 | |
---|
108 | USE grid_variables, & |
---|
109 | ONLY: ddx, dx, ddy, dy |
---|
110 | |
---|
111 | USE indices, & |
---|
112 | ONLY: nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg, & |
---|
113 | nysg, nyng |
---|
114 | |
---|
115 | USE kinds !< Set precision of INTEGER and REAL arrays according to PALM |
---|
116 | |
---|
117 | USE netcdf_data_input_mod, & |
---|
118 | ONLY: netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f, & |
---|
119 | uvem_irradiance_f, uvem_integration_f, building_obstruction_f |
---|
120 | ! |
---|
121 | !-- Import radiation model to obtain input for mean radiant temperature |
---|
122 | USE radiation_model_mod, & |
---|
123 | ONLY: ix, iy, iz, id, mrt_nlevels, mrt_include_sw, & |
---|
124 | mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation, & |
---|
125 | radiation_interactions, rad_sw_in, & |
---|
126 | rad_sw_out, rad_lw_in, rad_lw_out |
---|
127 | |
---|
128 | USE surface_mod, & |
---|
129 | ONLY: get_topography_top_index_ji |
---|
130 | |
---|
131 | IMPLICIT NONE |
---|
132 | |
---|
133 | PRIVATE |
---|
134 | |
---|
135 | ! |
---|
136 | !-- Declare all global variables within the module (alphabetical order) |
---|
137 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tmrt_grid !< tmrt results (degree_C) |
---|
138 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: perct !< PT results (degree_C) |
---|
139 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: utci !< UTCI results (degree_C) |
---|
140 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pet !< PET results (degree_C) |
---|
141 | ! |
---|
142 | !-- Grids for averaged thermal indices |
---|
143 | REAL(wp), DIMENSION(:), ALLOCATABLE :: mrt_av_grid !< time average mean |
---|
144 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: perct_av !< PT results (aver. input) (degree_C) |
---|
145 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: utci_av !< UTCI results (aver. input) (degree_C) |
---|
146 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pet_av !< PET results (aver. input) (degree_C) |
---|
147 | |
---|
148 | |
---|
149 | INTEGER( iwp ) :: bio_cell_level !< cell level biom calculates for |
---|
150 | REAL ( wp ) :: bio_output_height !< height output is calculated in m |
---|
151 | REAL ( wp ) :: time_bio_results !< the time, the last set of biom results have been calculated for |
---|
152 | REAL ( wp ), PARAMETER :: human_absorb = 0.7_wp !< SW absorbtivity of a human body (Fanger 1972) |
---|
153 | REAL ( wp ), PARAMETER :: human_emiss = 0.97_wp !< LW emissivity of a human body after (Fanger 1972) |
---|
154 | REAL ( wp ), PARAMETER :: bio_fill_value = -9999._wp !< set module fill value, replace by global fill value as soon as available |
---|
155 | ! |
---|
156 | !-- |
---|
157 | LOGICAL :: aver_perct = .FALSE. !< switch: do perct averaging in this module? (if .FALSE. this is done globally) |
---|
158 | LOGICAL :: aver_q = .FALSE. !< switch: do e averaging in this module? |
---|
159 | LOGICAL :: aver_u = .FALSE. !< switch: do u averaging in this module? |
---|
160 | LOGICAL :: aver_v = .FALSE. !< switch: do v averaging in this module? |
---|
161 | LOGICAL :: aver_w = .FALSE. !< switch: do w averaging in this module? |
---|
162 | LOGICAL :: average_trigger_perct = .FALSE. !< update averaged input on call to bio_perct? |
---|
163 | LOGICAL :: average_trigger_utci = .FALSE. !< update averaged input on call to bio_utci? |
---|
164 | LOGICAL :: average_trigger_pet = .FALSE. !< update averaged input on call to bio_pet? |
---|
165 | |
---|
166 | LOGICAL :: thermal_comfort = .TRUE. !< Turn all thermal indices on or off |
---|
167 | LOGICAL :: bio_perct = .TRUE. !< Turn index PT (instant. input) on or off |
---|
168 | LOGICAL :: bio_perct_av = .TRUE. !< Turn index PT (averaged input) on or off |
---|
169 | LOGICAL :: bio_pet = .TRUE. !< Turn index PET (instant. input) on or off |
---|
170 | LOGICAL :: bio_pet_av = .TRUE. !< Turn index PET (averaged input) on or off |
---|
171 | LOGICAL :: bio_utci = .TRUE. !< Turn index UTCI (instant. input) on or off |
---|
172 | LOGICAL :: bio_utci_av = .TRUE. !< Turn index UTCI (averaged input) on or off |
---|
173 | |
---|
174 | ! |
---|
175 | !-- UVEM parameters from here |
---|
176 | ! |
---|
177 | !-- Declare all global variables within the module (alphabetical order) |
---|
178 | INTEGER(iwp) :: ai = 0 !< loop index in azimuth direction |
---|
179 | INTEGER(iwp) :: bi = 0 !< loop index of bit location within an 8bit-integer (one Byte) |
---|
180 | INTEGER(iwp) :: clothing = 1 !< clothing (0=unclothed, 1=Arms,Hands,Face free, 3=Hand,Face free) |
---|
181 | INTEGER(iwp) :: iq = 0 !< loop index of irradiance quantity |
---|
182 | INTEGER(iwp) :: pobi = 0 !< loop index of the position of corresponding byte within ibset byte vektor |
---|
183 | INTEGER(iwp) :: obstruction_direct_beam = 0 !< Obstruction information for direct beam |
---|
184 | INTEGER(iwp) :: zi = 0 !< loop index in zenith direction |
---|
185 | |
---|
186 | INTEGER(KIND=1), DIMENSION(0:44) :: obstruction_temp1 = 0 !< temporary obstruction information stored with ibset |
---|
187 | INTEGER(iwp), DIMENSION(0:359) :: obstruction_temp2 = 0 !< restored temporary obstruction information from ibset file |
---|
188 | |
---|
189 | INTEGER(iwp), DIMENSION(0:35,0:9) :: obstruction = 1 !< final 2D obstruction information array |
---|
190 | |
---|
191 | LOGICAL :: consider_obstructions = .TRUE. !< namelist parameter (see documentation) |
---|
192 | LOGICAL :: sun_in_south = .FALSE. !< namelist parameter (see documentation) |
---|
193 | LOGICAL :: turn_to_sun = .TRUE. !< namelist parameter (see documentation) |
---|
194 | LOGICAL :: uv_exposure = .FALSE. !< namelist parameter (see documentation) |
---|
195 | |
---|
196 | REAL(wp) :: diffuse_exposure = 0.0_wp !< calculated exposure by diffuse radiation |
---|
197 | REAL(wp) :: direct_exposure = 0.0_wp !< calculated exposure by direct solar beam |
---|
198 | REAL(wp) :: orientation_angle = 0.0_wp !< orientation of front/face of the human model |
---|
199 | REAL(wp) :: projection_area_direct_beam = 0.0_wp !< projection area for direct solar beam |
---|
200 | REAL(wp) :: saa = 180.0_wp !< solar azimuth angle |
---|
201 | REAL(wp) :: startpos_human = 0.0_wp !< start value for azimuth interpolation of human geometry array |
---|
202 | REAL(wp) :: startpos_saa_float = 0.0_wp !< start value for azimuth interpolation of radiance array |
---|
203 | REAL(wp) :: sza = 20.0_wp !< solar zenith angle |
---|
204 | REAL(wp) :: xfactor = 0.0_wp !< relative x-position used for interpolation |
---|
205 | REAL(wp) :: yfactor = 0.0_wp !< relative y-position used for interpolation |
---|
206 | |
---|
207 | REAL(wp), DIMENSION(0:2) :: irradiance = 0.0_wp !< iradiance values extracted from irradiance lookup table |
---|
208 | |
---|
209 | REAL(wp), DIMENSION(0:2,0:90) :: irradiance_lookup_table = 0.0_wp !< irradiance lookup table |
---|
210 | REAL(wp), DIMENSION(0:35,0:9) :: integration_array = 0.0_wp !< solid angle factors for hemispherical integration |
---|
211 | REAL(wp), DIMENSION(0:35,0:9) :: projection_area = 0.0_wp !< projection areas of a human (all directions) |
---|
212 | REAL(wp), DIMENSION(0:35,0:9) :: projection_area_lookup_table = 0.0_wp !< human geometry lookup table (projection areas) |
---|
213 | REAL(wp), DIMENSION(0:71,0:9) :: projection_area_direct_temp = 0.0_wp !< temporary projection area for direct solar beam |
---|
214 | REAL(wp), DIMENSION(0:71,0:9) :: projection_area_temp = 0.0_wp !< temporary projection area for all directions |
---|
215 | REAL(wp), DIMENSION(0:35,0:9) :: radiance_array = 0.0_wp !< radiance extracted from radiance_lookup_table |
---|
216 | REAL(wp), DIMENSION(0:71,0:9) :: radiance_array_temp = 0.0_wp !< temporary radiance data |
---|
217 | |
---|
218 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vitd3_exposure !< result variable for instantaneous vitamin-D weighted exposures |
---|
219 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vitd3_exposure_av !< result variable for summation of vitamin-D weighted exposures |
---|
220 | |
---|
221 | REAL(wp), DIMENSION(0:35,0:9,0:90) :: radiance_lookup_table = 0.0_wp !< radiance lookup table |
---|
222 | |
---|
223 | ! |
---|
224 | !-- INTERFACES that must be available to other modules (alphabetical order) |
---|
225 | |
---|
226 | PUBLIC bio_3d_data_averaging, bio_check_data_output, & |
---|
227 | bio_calculate_mrt_grid, bio_calculate_thermal_index_maps, bio_calc_ipt, & |
---|
228 | bio_check_parameters, bio_data_output_3d, bio_data_output_2d, & |
---|
229 | bio_define_netcdf_grid, bio_get_thermal_index_input_ij, bio_header, & |
---|
230 | bio_init, bio_parin, bio_perct, bio_perct_av, bio_pet, & |
---|
231 | bio_pet_av, bio_utci, bio_utci_av, thermal_comfort, time_bio_results, & |
---|
232 | ! |
---|
233 | !-- UVEM PUBLIC variables and methods |
---|
234 | uvem_calc_exposure, uv_exposure |
---|
235 | |
---|
236 | ! |
---|
237 | !-- PALM interfaces: |
---|
238 | ! |
---|
239 | !-- 3D averaging for HTCM _INPUT_ variables |
---|
240 | INTERFACE bio_3d_data_averaging |
---|
241 | MODULE PROCEDURE bio_3d_data_averaging |
---|
242 | END INTERFACE bio_3d_data_averaging |
---|
243 | ! |
---|
244 | !-- Calculate mtr from rtm fluxes and assign into 2D grid |
---|
245 | INTERFACE bio_calculate_mrt_grid |
---|
246 | MODULE PROCEDURE bio_calculate_mrt_grid |
---|
247 | END INTERFACE bio_calculate_mrt_grid |
---|
248 | ! |
---|
249 | !-- Calculate static thermal indices PT, UTCI and/or PET |
---|
250 | INTERFACE bio_calculate_thermal_index_maps |
---|
251 | MODULE PROCEDURE bio_calculate_thermal_index_maps |
---|
252 | END INTERFACE bio_calculate_thermal_index_maps |
---|
253 | ! |
---|
254 | !-- Calculate the dynamic index iPT (to be caled by the agent model) |
---|
255 | INTERFACE bio_calc_ipt |
---|
256 | MODULE PROCEDURE bio_calc_ipt |
---|
257 | END INTERFACE bio_calc_ipt |
---|
258 | ! |
---|
259 | !-- Data output checks for 2D/3D data to be done in check_parameters |
---|
260 | INTERFACE bio_check_data_output |
---|
261 | MODULE PROCEDURE bio_check_data_output |
---|
262 | END INTERFACE bio_check_data_output |
---|
263 | ! |
---|
264 | !-- Input parameter checks to be done in check_parameters |
---|
265 | INTERFACE bio_check_parameters |
---|
266 | MODULE PROCEDURE bio_check_parameters |
---|
267 | END INTERFACE bio_check_parameters |
---|
268 | ! |
---|
269 | !-- Data output of 2D quantities |
---|
270 | INTERFACE bio_data_output_2d |
---|
271 | MODULE PROCEDURE bio_data_output_2d |
---|
272 | END INTERFACE bio_data_output_2d |
---|
273 | ! |
---|
274 | !-- no 3D data, thus, no averaging of 3D data, removed |
---|
275 | INTERFACE bio_data_output_3d |
---|
276 | MODULE PROCEDURE bio_data_output_3d |
---|
277 | END INTERFACE bio_data_output_3d |
---|
278 | ! |
---|
279 | !-- Definition of data output quantities |
---|
280 | INTERFACE bio_define_netcdf_grid |
---|
281 | MODULE PROCEDURE bio_define_netcdf_grid |
---|
282 | END INTERFACE bio_define_netcdf_grid |
---|
283 | ! |
---|
284 | !-- Obtains all relevant input values to estimate local thermal comfort/stress |
---|
285 | INTERFACE bio_get_thermal_index_input_ij |
---|
286 | MODULE PROCEDURE bio_get_thermal_index_input_ij |
---|
287 | END INTERFACE bio_get_thermal_index_input_ij |
---|
288 | ! |
---|
289 | !-- Output of information to the header file |
---|
290 | INTERFACE bio_header |
---|
291 | MODULE PROCEDURE bio_header |
---|
292 | END INTERFACE bio_header |
---|
293 | ! |
---|
294 | !-- Initialization actions |
---|
295 | INTERFACE bio_init |
---|
296 | MODULE PROCEDURE bio_init |
---|
297 | END INTERFACE bio_init |
---|
298 | ! |
---|
299 | !-- Reading of NAMELIST parameters |
---|
300 | INTERFACE bio_parin |
---|
301 | MODULE PROCEDURE bio_parin |
---|
302 | END INTERFACE bio_parin |
---|
303 | |
---|
304 | ! |
---|
305 | !-- Calculate UV exposure grid |
---|
306 | INTERFACE uvem_calc_exposure |
---|
307 | MODULE PROCEDURE uvem_calc_exposure |
---|
308 | END INTERFACE uvem_calc_exposure |
---|
309 | |
---|
310 | CONTAINS |
---|
311 | |
---|
312 | |
---|
313 | !------------------------------------------------------------------------------! |
---|
314 | ! Description: |
---|
315 | ! ------------ |
---|
316 | !> Sum up and time-average biom input quantities as well as allocate |
---|
317 | !> the array necessary for storing the average. |
---|
318 | !> There is a considerable difference to the 3d_data_averaging subroutines |
---|
319 | !> used by other modules: |
---|
320 | !> For the thermal indices, the module needs to average the input conditions |
---|
321 | !> not the result! |
---|
322 | !------------------------------------------------------------------------------! |
---|
323 | SUBROUTINE bio_3d_data_averaging( mode, variable ) |
---|
324 | |
---|
325 | IMPLICIT NONE |
---|
326 | |
---|
327 | CHARACTER (LEN=*) :: mode !< averaging mode: allocate, sum, or average |
---|
328 | CHARACTER (LEN=*) :: variable !< The variable in question |
---|
329 | |
---|
330 | INTEGER(iwp) :: i !< Running index, x-dir |
---|
331 | INTEGER(iwp) :: j !< Running index, y-dir |
---|
332 | INTEGER(iwp) :: k !< Running index, z-dir |
---|
333 | |
---|
334 | |
---|
335 | IF ( mode == 'allocate' ) THEN |
---|
336 | |
---|
337 | SELECT CASE ( TRIM( variable ) ) |
---|
338 | |
---|
339 | CASE ( 'bio_mrt' ) |
---|
340 | IF ( .NOT. ALLOCATED( mrt_av_grid ) ) THEN |
---|
341 | ALLOCATE( mrt_av_grid(nmrtbl) ) |
---|
342 | ENDIF |
---|
343 | mrt_av_grid = 0.0_wp |
---|
344 | |
---|
345 | CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' ) |
---|
346 | |
---|
347 | ! |
---|
348 | !-- Averaging, as well as the allocation of the required grids must be |
---|
349 | ! done only once, independent from for how many thermal indices |
---|
350 | ! averaged output is desired. |
---|
351 | ! Therefore wee need to memorize which index is the one that controls |
---|
352 | ! the averaging (what must be the first thermal index called). |
---|
353 | ! Indices are in unknown order as depending on the input file, |
---|
354 | ! determine first index to average und update only once |
---|
355 | |
---|
356 | !-- Only proceed here if this was not done for any index before. This |
---|
357 | ! is done only once during the whole model run. |
---|
358 | IF ( .NOT. average_trigger_perct .AND. & |
---|
359 | .NOT. average_trigger_utci .AND. & |
---|
360 | .NOT. average_trigger_pet ) THEN |
---|
361 | ! |
---|
362 | !-- Allocate the required grids |
---|
363 | IF ( .NOT. ALLOCATED( perct_av ) ) THEN |
---|
364 | ALLOCATE( perct_av (nys:nyn,nxl:nxr) ) |
---|
365 | ENDIF |
---|
366 | perct_av = REAL( bio_fill_value, KIND = wp ) |
---|
367 | |
---|
368 | IF ( .NOT. ALLOCATED( utci_av ) ) THEN |
---|
369 | ALLOCATE( utci_av (nys:nyn,nxl:nxr) ) |
---|
370 | ENDIF |
---|
371 | utci_av = REAL( bio_fill_value, KIND = wp ) |
---|
372 | |
---|
373 | IF ( .NOT. ALLOCATED( pet_av ) ) THEN |
---|
374 | ALLOCATE( pet_av (nys:nyn,nxl:nxr) ) |
---|
375 | ENDIF |
---|
376 | pet_av = REAL( bio_fill_value, KIND = wp ) |
---|
377 | ! |
---|
378 | !-- Memorize the first index called to control averaging |
---|
379 | IF ( TRIM( variable ) == 'bio_perct*' ) THEN |
---|
380 | average_trigger_perct = .TRUE. |
---|
381 | ENDIF |
---|
382 | IF ( TRIM( variable ) == 'bio_utci*' ) THEN |
---|
383 | average_trigger_utci = .TRUE. |
---|
384 | ENDIF |
---|
385 | IF ( TRIM( variable ) == 'bio_pet*' ) THEN |
---|
386 | average_trigger_pet = .TRUE. |
---|
387 | ENDIF |
---|
388 | ENDIF |
---|
389 | ! |
---|
390 | !-- Only continue if var is the index, that controls averaging. |
---|
391 | ! Break immediatelly (doing nothing) for the other indices. |
---|
392 | IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') RETURN |
---|
393 | IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*') RETURN |
---|
394 | IF ( average_trigger_pet .AND. TRIM( variable ) /= 'bio_pet*') RETURN |
---|
395 | ! |
---|
396 | !-- Now memorize which of the input grids are not averaged by other |
---|
397 | ! modules. Set averaging switch to .TRUE. in that case. |
---|
398 | IF ( .NOT. ALLOCATED( pt_av ) ) THEN |
---|
399 | ALLOCATE( pt_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) |
---|
400 | aver_perct = .TRUE. |
---|
401 | pt_av = 0.0_wp |
---|
402 | ENDIF |
---|
403 | |
---|
404 | IF ( .NOT. ALLOCATED( q_av ) ) THEN |
---|
405 | ALLOCATE( q_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) |
---|
406 | aver_q = .TRUE. |
---|
407 | q_av = 0.0_wp |
---|
408 | ENDIF |
---|
409 | |
---|
410 | IF ( .NOT. ALLOCATED( u_av ) ) THEN |
---|
411 | ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
412 | aver_u = .TRUE. |
---|
413 | u_av = 0.0_wp |
---|
414 | ENDIF |
---|
415 | |
---|
416 | IF ( .NOT. ALLOCATED( v_av ) ) THEN |
---|
417 | ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
418 | aver_v = .TRUE. |
---|
419 | v_av = 0.0_wp |
---|
420 | ENDIF |
---|
421 | |
---|
422 | IF ( .NOT. ALLOCATED( w_av ) ) THEN |
---|
423 | ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
424 | aver_w = .TRUE. |
---|
425 | w_av = 0.0_wp |
---|
426 | ENDIF |
---|
427 | |
---|
428 | CASE ( 'uvem_vitd3dose*' ) |
---|
429 | IF ( .NOT. ALLOCATED( vitd3_exposure_av ) ) THEN |
---|
430 | ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) ) |
---|
431 | ENDIF |
---|
432 | vitd3_exposure_av = 0.0_wp |
---|
433 | |
---|
434 | CASE DEFAULT |
---|
435 | CONTINUE |
---|
436 | |
---|
437 | END SELECT |
---|
438 | |
---|
439 | ELSEIF ( mode == 'sum' ) THEN |
---|
440 | |
---|
441 | SELECT CASE ( TRIM( variable ) ) |
---|
442 | |
---|
443 | CASE ( 'bio_mrt' ) |
---|
444 | IF ( ALLOCATED( mrt_av_grid ) ) THEN |
---|
445 | |
---|
446 | IF ( mrt_include_sw ) THEN |
---|
447 | mrt_av_grid(:) = mrt_av_grid(:) + & |
---|
448 | (( human_absorb * mrtinsw(:) + human_emiss * mrtinlw(:)) & |
---|
449 | / (human_emiss * sigma_sb)) ** .25_wp - degc_to_k |
---|
450 | ELSE |
---|
451 | mrt_av_grid(:) = mrt_av_grid(:) + & |
---|
452 | (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp & |
---|
453 | - degc_to_k |
---|
454 | ENDIF |
---|
455 | ENDIF |
---|
456 | |
---|
457 | CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' ) |
---|
458 | ! |
---|
459 | !-- Only continue if updateindex, see above |
---|
460 | IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') & |
---|
461 | RETURN |
---|
462 | IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*') & |
---|
463 | RETURN |
---|
464 | IF ( average_trigger_pet .AND. TRIM( variable ) /= 'bio_pet*') & |
---|
465 | RETURN |
---|
466 | |
---|
467 | IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN |
---|
468 | DO i = nxl, nxr |
---|
469 | DO j = nys, nyn |
---|
470 | DO k = nzb, nzt+1 |
---|
471 | pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i) |
---|
472 | ENDDO |
---|
473 | ENDDO |
---|
474 | ENDDO |
---|
475 | ENDIF |
---|
476 | |
---|
477 | IF ( ALLOCATED( q_av ) .AND. aver_q ) THEN |
---|
478 | DO i = nxl, nxr |
---|
479 | DO j = nys, nyn |
---|
480 | DO k = nzb, nzt+1 |
---|
481 | q_av(k,j,i) = q_av(k,j,i) + q(k,j,i) |
---|
482 | ENDDO |
---|
483 | ENDDO |
---|
484 | ENDDO |
---|
485 | ENDIF |
---|
486 | |
---|
487 | IF ( ALLOCATED( u_av ) .AND. aver_u ) THEN |
---|
488 | DO i = nxlg, nxrg !< yes, ghost points are required here! |
---|
489 | DO j = nysg, nyng |
---|
490 | DO k = nzb, nzt+1 |
---|
491 | u_av(k,j,i) = u_av(k,j,i) + u(k,j,i) |
---|
492 | ENDDO |
---|
493 | ENDDO |
---|
494 | ENDDO |
---|
495 | ENDIF |
---|
496 | |
---|
497 | IF ( ALLOCATED( v_av ) .AND. aver_v ) THEN |
---|
498 | DO i = nxlg, nxrg !< yes, ghost points are required here! |
---|
499 | DO j = nysg, nyng |
---|
500 | DO k = nzb, nzt+1 |
---|
501 | v_av(k,j,i) = v_av(k,j,i) + v(k,j,i) |
---|
502 | ENDDO |
---|
503 | ENDDO |
---|
504 | ENDDO |
---|
505 | ENDIF |
---|
506 | |
---|
507 | IF ( ALLOCATED( w_av ) .AND. aver_w ) THEN |
---|
508 | DO i = nxlg, nxrg !< yes, ghost points are required here! |
---|
509 | DO j = nysg, nyng |
---|
510 | DO k = nzb, nzt+1 |
---|
511 | w_av(k,j,i) = w_av(k,j,i) + w(k,j,i) |
---|
512 | ENDDO |
---|
513 | ENDDO |
---|
514 | ENDDO |
---|
515 | ENDIF |
---|
516 | ! |
---|
517 | !-- This is a cumulated dose. No mode == 'average' for this quantity. |
---|
518 | CASE ( 'uvem_vitd3dose*' ) |
---|
519 | IF ( ALLOCATED( vitd3_exposure_av ) ) THEN |
---|
520 | DO i = nxlg, nxrg |
---|
521 | DO j = nysg, nyng |
---|
522 | vitd3_exposure_av(j,i) = vitd3_exposure_av(j,i) + vitd3_exposure(j,i) |
---|
523 | ENDDO |
---|
524 | ENDDO |
---|
525 | ENDIF |
---|
526 | |
---|
527 | CASE DEFAULT |
---|
528 | CONTINUE |
---|
529 | |
---|
530 | END SELECT |
---|
531 | |
---|
532 | ELSEIF ( mode == 'average' ) THEN |
---|
533 | |
---|
534 | SELECT CASE ( TRIM( variable ) ) |
---|
535 | |
---|
536 | CASE ( 'bio_mrt' ) |
---|
537 | IF ( ALLOCATED( mrt_av_grid ) ) THEN |
---|
538 | mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp ) |
---|
539 | ENDIF |
---|
540 | |
---|
541 | CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' ) |
---|
542 | ! |
---|
543 | !-- Only continue if update index, see above |
---|
544 | IF ( average_trigger_perct .AND. & |
---|
545 | TRIM( variable ) /= 'bio_perct*' ) RETURN |
---|
546 | IF ( average_trigger_utci .AND. & |
---|
547 | TRIM( variable ) /= 'bio_utci*' ) RETURN |
---|
548 | IF ( average_trigger_pet .AND. & |
---|
549 | TRIM( variable ) /= 'bio_pet*' ) RETURN |
---|
550 | |
---|
551 | IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN |
---|
552 | DO i = nxl, nxr |
---|
553 | DO j = nys, nyn |
---|
554 | DO k = nzb, nzt+1 |
---|
555 | pt_av(k,j,i) = pt_av(k,j,i) / & |
---|
556 | REAL( average_count_3d, KIND=wp ) |
---|
557 | ENDDO |
---|
558 | ENDDO |
---|
559 | ENDDO |
---|
560 | ENDIF |
---|
561 | |
---|
562 | IF ( ALLOCATED( q_av ) .AND. aver_q ) THEN |
---|
563 | DO i = nxl, nxr |
---|
564 | DO j = nys, nyn |
---|
565 | DO k = nzb, nzt+1 |
---|
566 | q_av(k,j,i) = q_av(k,j,i) / & |
---|
567 | REAL( average_count_3d, KIND=wp ) |
---|
568 | ENDDO |
---|
569 | ENDDO |
---|
570 | ENDDO |
---|
571 | ENDIF |
---|
572 | |
---|
573 | IF ( ALLOCATED( u_av ) .AND. aver_u ) THEN |
---|
574 | DO i = nxlg, nxrg !< yes, ghost points are required here! |
---|
575 | DO j = nysg, nyng |
---|
576 | DO k = nzb, nzt+1 |
---|
577 | u_av(k,j,i) = u_av(k,j,i) / & |
---|
578 | REAL( average_count_3d, KIND=wp ) |
---|
579 | ENDDO |
---|
580 | ENDDO |
---|
581 | ENDDO |
---|
582 | ENDIF |
---|
583 | |
---|
584 | IF ( ALLOCATED( v_av ) .AND. aver_v ) THEN |
---|
585 | DO i = nxlg, nxrg |
---|
586 | DO j = nysg, nyng |
---|
587 | DO k = nzb, nzt+1 |
---|
588 | v_av(k,j,i) = v_av(k,j,i) / & |
---|
589 | REAL( average_count_3d, KIND=wp ) |
---|
590 | ENDDO |
---|
591 | ENDDO |
---|
592 | ENDDO |
---|
593 | ENDIF |
---|
594 | |
---|
595 | IF ( ALLOCATED( w_av ) .AND. aver_w ) THEN |
---|
596 | DO i = nxlg, nxrg |
---|
597 | DO j = nysg, nyng |
---|
598 | DO k = nzb, nzt+1 |
---|
599 | w_av(k,j,i) = w_av(k,j,i) / & |
---|
600 | REAL( average_count_3d, KIND=wp ) |
---|
601 | ENDDO |
---|
602 | ENDDO |
---|
603 | ENDDO |
---|
604 | ENDIF |
---|
605 | ! |
---|
606 | !-- Udate all thermal index grids with updated averaged input |
---|
607 | CALL bio_calculate_thermal_index_maps ( .TRUE. ) |
---|
608 | |
---|
609 | ! |
---|
610 | !-- No averaging for UVEM since we are calculating a dose (only sum is |
---|
611 | !-- calculated and saved to av.nc file) |
---|
612 | |
---|
613 | END SELECT |
---|
614 | |
---|
615 | ENDIF |
---|
616 | |
---|
617 | |
---|
618 | END SUBROUTINE bio_3d_data_averaging |
---|
619 | |
---|
620 | |
---|
621 | |
---|
622 | !------------------------------------------------------------------------------! |
---|
623 | ! Description: |
---|
624 | ! ------------ |
---|
625 | !> Check data output for biometeorology model |
---|
626 | !------------------------------------------------------------------------------! |
---|
627 | SUBROUTINE bio_check_data_output( var, unit, i, ilen, k ) |
---|
628 | |
---|
629 | USE control_parameters, & |
---|
630 | ONLY: data_output, message_string |
---|
631 | |
---|
632 | IMPLICIT NONE |
---|
633 | |
---|
634 | CHARACTER (LEN=*) :: unit !< The unit for the variable var |
---|
635 | CHARACTER (LEN=*) :: var !< The variable in question |
---|
636 | |
---|
637 | INTEGER(iwp) :: i !< |
---|
638 | INTEGER(iwp) :: ilen !< |
---|
639 | INTEGER(iwp) :: k !< |
---|
640 | |
---|
641 | SELECT CASE ( TRIM( var ) ) |
---|
642 | ! |
---|
643 | !-- Allocate a temporary array with the desired output dimensions. |
---|
644 | CASE ( 'bio_mrt' ) |
---|
645 | unit = 'degree_C' |
---|
646 | IF ( .NOT. ALLOCATED( tmrt_grid ) ) THEN |
---|
647 | ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) ) |
---|
648 | ENDIF |
---|
649 | tmrt_grid = REAL( bio_fill_value, KIND = wp ) |
---|
650 | |
---|
651 | CASE ( 'bio_perct*' ) |
---|
652 | unit = 'degree_C' |
---|
653 | IF ( .NOT. ALLOCATED( perct ) ) THEN |
---|
654 | ALLOCATE( perct (nys:nyn,nxl:nxr) ) |
---|
655 | ENDIF |
---|
656 | perct = REAL( bio_fill_value, KIND = wp ) |
---|
657 | |
---|
658 | CASE ( 'bio_utci*' ) |
---|
659 | unit = 'degree_C' |
---|
660 | IF ( .NOT. ALLOCATED( utci ) ) THEN |
---|
661 | ALLOCATE( utci (nys:nyn,nxl:nxr) ) |
---|
662 | ENDIF |
---|
663 | utci = REAL( bio_fill_value, KIND = wp ) |
---|
664 | |
---|
665 | CASE ( 'bio_pet*' ) |
---|
666 | unit = 'degree_C' |
---|
667 | IF ( .NOT. ALLOCATED( pet ) ) THEN |
---|
668 | ALLOCATE( pet (nys:nyn,nxl:nxr) ) |
---|
669 | ENDIF |
---|
670 | pet = REAL( bio_fill_value, KIND = wp ) |
---|
671 | |
---|
672 | CASE ( 'uvem_vitd3*' ) |
---|
673 | IF ( .NOT. uv_exposure ) THEN |
---|
674 | message_string = 'output of "' // TRIM( var ) // '" requi' // & |
---|
675 | 'res a namelist &uvexposure_par' |
---|
676 | CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 ) |
---|
677 | ENDIF |
---|
678 | IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN |
---|
679 | message_string = 'illegal value for data_output: "' // & |
---|
680 | TRIM( var ) // '" & only 2d-horizontal ' // & |
---|
681 | 'cross sections are allowed for this value' |
---|
682 | CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 ) |
---|
683 | ENDIF |
---|
684 | unit = 'IU/s' |
---|
685 | IF ( .NOT. ALLOCATED( vitd3_exposure ) ) THEN |
---|
686 | ALLOCATE( vitd3_exposure(nysg:nyng,nxlg:nxrg) ) |
---|
687 | ENDIF |
---|
688 | vitd3_exposure = 0.0_wp |
---|
689 | |
---|
690 | CASE ( 'uvem_vitd3dose*' ) |
---|
691 | IF ( .NOT. uv_exposure ) THEN |
---|
692 | message_string = 'output of "' // TRIM( var ) // '" requi' // & |
---|
693 | 'res a namelist &uvexposure_par' |
---|
694 | CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 ) |
---|
695 | ENDIF |
---|
696 | IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN |
---|
697 | message_string = 'illegal value for data_output: "' // & |
---|
698 | TRIM( var ) // '" & only 2d-horizontal ' // & |
---|
699 | 'cross sections are allowed for this value' |
---|
700 | CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 ) |
---|
701 | ENDIF |
---|
702 | unit = 'IU/av-h' |
---|
703 | IF ( .NOT. ALLOCATED( vitd3_exposure_av ) ) THEN |
---|
704 | ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) ) |
---|
705 | ENDIF |
---|
706 | vitd3_exposure_av = 0.0_wp |
---|
707 | |
---|
708 | CASE DEFAULT |
---|
709 | unit = 'illegal' |
---|
710 | |
---|
711 | END SELECT |
---|
712 | |
---|
713 | IF ( thermal_comfort .AND. unit == 'degree_C' ) THEN |
---|
714 | ! |
---|
715 | !-- Futher checks if output belongs to biometeorology |
---|
716 | ! Break if required modules "radiation" and "humidity" are not running. |
---|
717 | IF ( .NOT. radiation ) THEN |
---|
718 | message_string = 'output of "' // TRIM( var ) // '" require' & |
---|
719 | // 's radiation = .TRUE.' |
---|
720 | CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 ) |
---|
721 | unit = 'illegal' |
---|
722 | ENDIF |
---|
723 | IF ( mrt_nlevels == 0 ) THEN |
---|
724 | message_string = 'output of "' // TRIM( var ) // '" require' & |
---|
725 | // 's mrt_nlevels > 0' |
---|
726 | CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 ) |
---|
727 | unit = 'illegal' |
---|
728 | ENDIF |
---|
729 | |
---|
730 | |
---|
731 | ENDIF |
---|
732 | |
---|
733 | END SUBROUTINE bio_check_data_output |
---|
734 | |
---|
735 | !------------------------------------------------------------------------------! |
---|
736 | ! Description: |
---|
737 | ! ------------ |
---|
738 | !> Check parameters routine for biom module |
---|
739 | !> check_parameters 1302 |
---|
740 | !------------------------------------------------------------------------------! |
---|
741 | SUBROUTINE bio_check_parameters |
---|
742 | |
---|
743 | USE control_parameters, & |
---|
744 | ONLY: message_string |
---|
745 | |
---|
746 | |
---|
747 | IMPLICIT NONE |
---|
748 | |
---|
749 | ! |
---|
750 | !-- Disabled as long as radiation model not available |
---|
751 | |
---|
752 | IF ( thermal_comfort .AND. .NOT. humidity ) THEN |
---|
753 | message_string = 'The estimation of thermal comfort requires ' // & |
---|
754 | 'air humidity information, but humidity module ' // & |
---|
755 | 'is disabled!' |
---|
756 | CALL message( 'check_parameters', 'PA0561', 0, 0, 0, 6, 0 ) |
---|
757 | ENDIF |
---|
758 | |
---|
759 | END SUBROUTINE bio_check_parameters |
---|
760 | |
---|
761 | |
---|
762 | !------------------------------------------------------------------------------! |
---|
763 | ! |
---|
764 | ! Description: |
---|
765 | ! ------------ |
---|
766 | !> Subroutine defining 2D output variables |
---|
767 | !> data_output_2d 1188ff |
---|
768 | !------------------------------------------------------------------------------! |
---|
769 | SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf, & |
---|
770 | two_d, nzb_do, nzt_do ) |
---|
771 | |
---|
772 | |
---|
773 | USE kinds |
---|
774 | |
---|
775 | |
---|
776 | IMPLICIT NONE |
---|
777 | ! |
---|
778 | !-- Input variables |
---|
779 | CHARACTER (LEN=*), INTENT(IN) :: variable !< Char identifier to select var for output |
---|
780 | INTEGER(iwp), INTENT(IN) :: av !< Use averaged data? 0 = no, 1 = yes? |
---|
781 | INTEGER(iwp), INTENT(IN) :: nzb_do !< Unused. 2D. nz bottom to nz top |
---|
782 | INTEGER(iwp), INTENT(IN) :: nzt_do !< Unused. |
---|
783 | ! |
---|
784 | !-- Output variables |
---|
785 | CHARACTER (LEN=*), INTENT(OUT) :: grid !< Grid type (always "zu1" for biom) |
---|
786 | LOGICAL, INTENT(OUT) :: found !< Output found? |
---|
787 | LOGICAL, INTENT(OUT) :: two_d !< Flag parameter that indicates 2D variables, horizontal cross sections, must be .TRUE. |
---|
788 | REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< Temp. result grid to return |
---|
789 | ! |
---|
790 | !-- Internal variables |
---|
791 | INTEGER(iwp) :: i !< Running index, x-dir |
---|
792 | INTEGER(iwp) :: j !< Running index, y-dir |
---|
793 | INTEGER(iwp) :: k !< Running index, z-dir |
---|
794 | INTEGER(iwp) :: l !< Running index, radiation grid |
---|
795 | |
---|
796 | |
---|
797 | found = .TRUE. |
---|
798 | local_pf = bio_fill_value |
---|
799 | |
---|
800 | SELECT CASE ( TRIM( variable ) ) |
---|
801 | |
---|
802 | |
---|
803 | CASE ( 'bio_mrt_xy' ) |
---|
804 | grid = 'zu1' |
---|
805 | two_d = .FALSE. !< can be calculated for several levels |
---|
806 | local_pf = REAL( bio_fill_value, KIND = wp ) |
---|
807 | DO l = 1, nmrtbl |
---|
808 | i = mrtbl(ix,l) |
---|
809 | j = mrtbl(iy,l) |
---|
810 | k = mrtbl(iz,l) |
---|
811 | IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR. & |
---|
812 | i < nxl .OR. i > nxr ) CYCLE |
---|
813 | IF ( av == 0 ) THEN |
---|
814 | IF ( mrt_include_sw ) THEN |
---|
815 | local_pf(i,j,k) = ((human_absorb * mrtinsw(l) + & |
---|
816 | human_emiss * mrtinlw(l)) / & |
---|
817 | (human_emiss * sigma_sb)) ** .25_wp - degc_to_k |
---|
818 | ELSE |
---|
819 | local_pf(i,j,k) = (human_emiss * mrtinlw(l) / & |
---|
820 | sigma_sb) ** .25_wp - degc_to_k |
---|
821 | ENDIF |
---|
822 | ELSE |
---|
823 | local_pf(i,j,k) = mrt_av_grid(l) |
---|
824 | ENDIF |
---|
825 | ENDDO |
---|
826 | |
---|
827 | |
---|
828 | CASE ( 'bio_perct*_xy' ) ! 2d-array |
---|
829 | grid = 'zu1' |
---|
830 | two_d = .TRUE. |
---|
831 | IF ( av == 0 ) THEN |
---|
832 | DO i = nxl, nxr |
---|
833 | DO j = nys, nyn |
---|
834 | local_pf(i,j,nzb+1) = perct(j,i) |
---|
835 | ENDDO |
---|
836 | ENDDO |
---|
837 | ELSE |
---|
838 | DO i = nxl, nxr |
---|
839 | DO j = nys, nyn |
---|
840 | local_pf(i,j,nzb+1) = perct_av(j,i) |
---|
841 | ENDDO |
---|
842 | ENDDO |
---|
843 | END IF |
---|
844 | |
---|
845 | |
---|
846 | CASE ( 'bio_utci*_xy' ) ! 2d-array |
---|
847 | grid = 'zu1' |
---|
848 | two_d = .TRUE. |
---|
849 | IF ( av == 0 ) THEN |
---|
850 | DO i = nxl, nxr |
---|
851 | DO j = nys, nyn |
---|
852 | local_pf(i,j,nzb+1) = utci(j,i) |
---|
853 | ENDDO |
---|
854 | ENDDO |
---|
855 | ELSE |
---|
856 | DO i = nxl, nxr |
---|
857 | DO j = nys, nyn |
---|
858 | local_pf(i,j,nzb+1) = utci_av(j,i) |
---|
859 | ENDDO |
---|
860 | ENDDO |
---|
861 | END IF |
---|
862 | |
---|
863 | |
---|
864 | CASE ( 'bio_pet*_xy' ) ! 2d-array |
---|
865 | grid = 'zu1' |
---|
866 | two_d = .TRUE. |
---|
867 | IF ( av == 0 ) THEN |
---|
868 | DO i = nxl, nxr |
---|
869 | DO j = nys, nyn |
---|
870 | local_pf(i,j,nzb+1) = pet(j,i) |
---|
871 | ENDDO |
---|
872 | ENDDO |
---|
873 | ELSE |
---|
874 | DO i = nxl, nxr |
---|
875 | DO j = nys, nyn |
---|
876 | local_pf(i,j,nzb+1) = pet_av(j,i) |
---|
877 | ENDDO |
---|
878 | ENDDO |
---|
879 | END IF |
---|
880 | |
---|
881 | ! |
---|
882 | !-- Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein. |
---|
883 | !-- However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged. |
---|
884 | CASE ( 'uvem_vitd3*_xy' ) ! 2d-array |
---|
885 | IF ( av == 0 ) THEN |
---|
886 | DO i = nxl, nxr |
---|
887 | DO j = nys, nyn |
---|
888 | local_pf(i,j,nzb+1) = vitd3_exposure(j,i) |
---|
889 | ENDDO |
---|
890 | ENDDO |
---|
891 | ENDIF |
---|
892 | |
---|
893 | two_d = .TRUE. |
---|
894 | grid = 'zu1' |
---|
895 | |
---|
896 | CASE ( 'uvem_vitd3dose*_xy' ) ! 2d-array |
---|
897 | IF ( av == 1 ) THEN |
---|
898 | DO i = nxl, nxr |
---|
899 | DO j = nys, nyn |
---|
900 | local_pf(i,j,nzb+1) = vitd3_exposure_av(j,i) |
---|
901 | ENDDO |
---|
902 | ENDDO |
---|
903 | ENDIF |
---|
904 | |
---|
905 | two_d = .TRUE. |
---|
906 | grid = 'zu1' |
---|
907 | |
---|
908 | |
---|
909 | CASE DEFAULT |
---|
910 | found = .FALSE. |
---|
911 | grid = 'none' |
---|
912 | |
---|
913 | END SELECT |
---|
914 | |
---|
915 | |
---|
916 | END SUBROUTINE bio_data_output_2d |
---|
917 | |
---|
918 | |
---|
919 | !------------------------------------------------------------------------------! |
---|
920 | ! |
---|
921 | ! Description: |
---|
922 | ! ------------ |
---|
923 | !> Subroutine defining 3D output variables (dummy, always 2d!) |
---|
924 | !> data_output_3d 709ff |
---|
925 | !------------------------------------------------------------------------------! |
---|
926 | SUBROUTINE bio_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) |
---|
927 | |
---|
928 | USE indices |
---|
929 | |
---|
930 | USE kinds |
---|
931 | |
---|
932 | |
---|
933 | IMPLICIT NONE |
---|
934 | ! |
---|
935 | !-- Input variables |
---|
936 | CHARACTER (LEN=*), INTENT(IN) :: variable !< Char identifier to select var for output |
---|
937 | INTEGER(iwp), INTENT(IN) :: av !< Use averaged data? 0 = no, 1 = yes? |
---|
938 | INTEGER(iwp), INTENT(IN) :: nzb_do !< Unused. 2D. nz bottom to nz top |
---|
939 | INTEGER(iwp), INTENT(IN) :: nzt_do !< Unused. |
---|
940 | ! |
---|
941 | !-- Output variables |
---|
942 | LOGICAL, INTENT(OUT) :: found !< Output found? |
---|
943 | REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< Temp. result grid to return |
---|
944 | ! |
---|
945 | !-- Internal variables |
---|
946 | INTEGER(iwp) :: l !< Running index, radiation grid |
---|
947 | INTEGER(iwp) :: i !< Running index, x-dir |
---|
948 | INTEGER(iwp) :: j !< Running index, y-dir |
---|
949 | INTEGER(iwp) :: k !< Running index, z-dir |
---|
950 | |
---|
951 | ! REAL(wp) :: mrt !< Buffer for mean radiant temperature |
---|
952 | |
---|
953 | found = .TRUE. |
---|
954 | |
---|
955 | SELECT CASE ( TRIM( variable ) ) |
---|
956 | |
---|
957 | CASE ( 'bio_mrt' ) |
---|
958 | local_pf = REAL( bio_fill_value, KIND = sp ) |
---|
959 | DO l = 1, nmrtbl |
---|
960 | i = mrtbl(ix,l) |
---|
961 | j = mrtbl(iy,l) |
---|
962 | k = mrtbl(iz,l) |
---|
963 | IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR. & |
---|
964 | i < nxl .OR. i > nxr ) CYCLE |
---|
965 | IF ( av == 0 ) THEN |
---|
966 | IF ( mrt_include_sw ) THEN |
---|
967 | local_pf(i,j,k) = REAL(((human_absorb * mrtinsw(l) + & |
---|
968 | human_emiss * mrtinlw(l)) / & |
---|
969 | (human_emiss * sigma_sb)) ** .25_wp - & |
---|
970 | degc_to_k, kind=sp ) |
---|
971 | ELSE |
---|
972 | local_pf(i,j,k) = REAL((human_emiss * mrtinlw(l) / & |
---|
973 | sigma_sb) ** .25_wp - degc_to_k, kind=sp ) !< why not (human_emiss * sigma_sb) as above? |
---|
974 | ENDIF |
---|
975 | ELSE |
---|
976 | local_pf(i,j,k) = REAL(mrt_av_grid(l), kind=sp) |
---|
977 | ENDIF |
---|
978 | ENDDO |
---|
979 | |
---|
980 | CASE DEFAULT |
---|
981 | found = .FALSE. |
---|
982 | |
---|
983 | END SELECT |
---|
984 | |
---|
985 | END SUBROUTINE bio_data_output_3d |
---|
986 | |
---|
987 | !------------------------------------------------------------------------------! |
---|
988 | ! Description: |
---|
989 | ! ------------ |
---|
990 | !> Subroutine defining appropriate grid for netcdf variables. |
---|
991 | !> It is called out from subroutine netcdf_interface_mod. |
---|
992 | !> netcdf_interface_mod 918ff |
---|
993 | !------------------------------------------------------------------------------! |
---|
994 | SUBROUTINE bio_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) |
---|
995 | |
---|
996 | IMPLICIT NONE |
---|
997 | ! |
---|
998 | !-- Input variables |
---|
999 | CHARACTER (LEN=*), INTENT(IN) :: var !< Name of output variable |
---|
1000 | ! |
---|
1001 | !-- Output variables |
---|
1002 | CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< x grid of output variable |
---|
1003 | CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< y grid of output variable |
---|
1004 | CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< z grid of output variable |
---|
1005 | |
---|
1006 | LOGICAL, INTENT(OUT) :: found !< Flag if output var is found |
---|
1007 | ! |
---|
1008 | !-- Local variables |
---|
1009 | LOGICAL :: is2d !< Var is 2d? |
---|
1010 | |
---|
1011 | INTEGER(iwp) :: l !< Length of the var array |
---|
1012 | |
---|
1013 | |
---|
1014 | found = .FALSE. |
---|
1015 | grid_x = 'none' |
---|
1016 | grid_y = 'none' |
---|
1017 | grid_z = 'none' |
---|
1018 | |
---|
1019 | l = MAX( 2, LEN_TRIM( var ) ) |
---|
1020 | is2d = ( var(l-1:l) == 'xy' ) |
---|
1021 | |
---|
1022 | IF ( var(1:4) == 'bio_' ) THEN |
---|
1023 | found = .TRUE. |
---|
1024 | grid_x = 'x' |
---|
1025 | grid_y = 'y' |
---|
1026 | grid_z = 'zu' |
---|
1027 | IF ( is2d .AND. var(1:7) /= 'bio_mrt' ) grid_z = 'zu1' |
---|
1028 | ENDIF |
---|
1029 | |
---|
1030 | IF ( is2d .AND. var(1:4) == 'uvem' ) THEN |
---|
1031 | grid_x = 'x' |
---|
1032 | grid_y = 'y' |
---|
1033 | grid_z = 'zu1' |
---|
1034 | ENDIF |
---|
1035 | |
---|
1036 | END SUBROUTINE bio_define_netcdf_grid |
---|
1037 | |
---|
1038 | !------------------------------------------------------------------------------! |
---|
1039 | ! Description: |
---|
1040 | ! ------------ |
---|
1041 | !> Header output for biom module |
---|
1042 | !> header 982 |
---|
1043 | !------------------------------------------------------------------------------! |
---|
1044 | SUBROUTINE bio_header( io ) |
---|
1045 | |
---|
1046 | IMPLICIT NONE |
---|
1047 | ! |
---|
1048 | !-- Input variables |
---|
1049 | INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file |
---|
1050 | ! |
---|
1051 | !-- Internal variables |
---|
1052 | CHARACTER (LEN=86) :: output_height_chr !< String for output height |
---|
1053 | |
---|
1054 | WRITE( output_height_chr, '(F8.1,7X)' ) bio_output_height |
---|
1055 | ! |
---|
1056 | !-- Write biom header |
---|
1057 | WRITE( io, 1 ) |
---|
1058 | WRITE( io, 2 ) TRIM( output_height_chr ) |
---|
1059 | WRITE( io, 3 ) TRIM( ACHAR( bio_cell_level ) ) |
---|
1060 | |
---|
1061 | 1 FORMAT (//' Human thermal comfort module information:'/ & |
---|
1062 | ' ------------------------------'/) |
---|
1063 | 2 FORMAT (' --> All indices calculated for a height of (m): ', A ) |
---|
1064 | 3 FORMAT (' --> This corresponds to cell level : ', A ) |
---|
1065 | |
---|
1066 | END SUBROUTINE bio_header |
---|
1067 | |
---|
1068 | |
---|
1069 | !------------------------------------------------------------------------------! |
---|
1070 | ! Description: |
---|
1071 | ! ------------ |
---|
1072 | !> Initialization of the HTCM |
---|
1073 | !> init_3d_model 1987ff |
---|
1074 | !------------------------------------------------------------------------------! |
---|
1075 | SUBROUTINE bio_init |
---|
1076 | |
---|
1077 | USE control_parameters, & |
---|
1078 | ONLY: message_string |
---|
1079 | |
---|
1080 | USE netcdf_data_input_mod, & |
---|
1081 | ONLY: netcdf_data_input_uvem |
---|
1082 | |
---|
1083 | IMPLICIT NONE |
---|
1084 | ! |
---|
1085 | !-- Internal vriables |
---|
1086 | REAL ( wp ) :: height !< current height in meters |
---|
1087 | ! |
---|
1088 | !-- Determine cell level corresponding to 1.1 m above ground level |
---|
1089 | ! (gravimetric center of sample human) |
---|
1090 | |
---|
1091 | time_bio_results = 0.0_wp |
---|
1092 | bio_cell_level = 0_iwp |
---|
1093 | bio_output_height = 0.5_wp * dz(1) |
---|
1094 | height = 0.0_wp |
---|
1095 | |
---|
1096 | bio_cell_level = INT ( 1.099_wp / dz(1) ) |
---|
1097 | bio_output_height = bio_output_height + bio_cell_level * dz(1) |
---|
1098 | |
---|
1099 | IF ( .NOT. radiation_interactions ) THEN |
---|
1100 | message_string = 'The mrt calculation requires ' // & |
---|
1101 | 'enabled radiation_interactions but it ' // & |
---|
1102 | 'is disabled!' |
---|
1103 | CALL message( 'check_parameters', 'PAHU03', 0, 0, -1, 6, 0 ) |
---|
1104 | ENDIF |
---|
1105 | |
---|
1106 | ! |
---|
1107 | !-- Init UVEM and load lookup tables |
---|
1108 | CALL netcdf_data_input_uvem |
---|
1109 | |
---|
1110 | END SUBROUTINE bio_init |
---|
1111 | |
---|
1112 | |
---|
1113 | !------------------------------------------------------------------------------! |
---|
1114 | ! Description: |
---|
1115 | ! ------------ |
---|
1116 | !> Parin for &biometeorology_parameters for reading biomet parameters |
---|
1117 | !------------------------------------------------------------------------------! |
---|
1118 | SUBROUTINE bio_parin |
---|
1119 | |
---|
1120 | IMPLICIT NONE |
---|
1121 | |
---|
1122 | ! |
---|
1123 | !-- Internal variables |
---|
1124 | CHARACTER (LEN=80) :: line !< Dummy string for current line in parameter file |
---|
1125 | |
---|
1126 | NAMELIST /biometeorology_parameters/ bio_pet, & |
---|
1127 | bio_pet_av, & |
---|
1128 | bio_perct, & |
---|
1129 | bio_perct_av, & |
---|
1130 | bio_utci, & |
---|
1131 | bio_utci_av, & |
---|
1132 | thermal_comfort, & |
---|
1133 | ! |
---|
1134 | !-- UVEM namelist parameters |
---|
1135 | clothing, & |
---|
1136 | consider_obstructions, & |
---|
1137 | orientation_angle, & |
---|
1138 | sun_in_south, & |
---|
1139 | turn_to_sun, & |
---|
1140 | uv_exposure |
---|
1141 | |
---|
1142 | |
---|
1143 | !-- Try to find biometeorology_parameters namelist |
---|
1144 | REWIND ( 11 ) |
---|
1145 | line = ' ' |
---|
1146 | DO WHILE ( INDEX( line, '&biometeorology_parameters' ) == 0 ) |
---|
1147 | READ ( 11, '(A)', END = 20 ) line |
---|
1148 | ENDDO |
---|
1149 | BACKSPACE ( 11 ) |
---|
1150 | |
---|
1151 | ! |
---|
1152 | !-- Read biometeorology_parameters namelist |
---|
1153 | READ ( 11, biometeorology_parameters, ERR = 10, END = 20 ) |
---|
1154 | |
---|
1155 | ! |
---|
1156 | !-- Set flag that indicates that the biomet_module is switched on |
---|
1157 | biometeorology = .TRUE. |
---|
1158 | |
---|
1159 | GOTO 20 |
---|
1160 | |
---|
1161 | ! |
---|
1162 | !-- In case of error |
---|
1163 | 10 BACKSPACE( 11 ) |
---|
1164 | READ( 11 , '(A)') line |
---|
1165 | CALL parin_fail_message( 'biometeorology_parameters', line ) |
---|
1166 | |
---|
1167 | ! |
---|
1168 | !-- Complete |
---|
1169 | 20 CONTINUE |
---|
1170 | |
---|
1171 | |
---|
1172 | END SUBROUTINE bio_parin |
---|
1173 | |
---|
1174 | !------------------------------------------------------------------------------! |
---|
1175 | ! Description: |
---|
1176 | ! ------------ |
---|
1177 | !> Calculate biometeorology MRT for all 2D grid |
---|
1178 | !------------------------------------------------------------------------------! |
---|
1179 | SUBROUTINE bio_calculate_mrt_grid ( av ) |
---|
1180 | |
---|
1181 | IMPLICIT NONE |
---|
1182 | |
---|
1183 | LOGICAL, INTENT(IN) :: av !< use averaged input? |
---|
1184 | ! |
---|
1185 | !-- Internal variables |
---|
1186 | INTEGER(iwp) :: i !< Running index, x-dir, radiation coordinates |
---|
1187 | INTEGER(iwp) :: j !< Running index, y-dir, radiation coordinates |
---|
1188 | INTEGER(iwp) :: k !< Running index, y-dir, radiation coordinates |
---|
1189 | INTEGER(iwp) :: l !< Running index, radiation coordinates |
---|
1190 | |
---|
1191 | ! |
---|
1192 | !-- Calculate biometeorology MRT from local radiation fluxes calculated by RTM and assign |
---|
1193 | !-- into 2D grid. Depending on selected output quantities, tmrt_grid might not have been |
---|
1194 | !-- allocated in bio_check_data_output yet. |
---|
1195 | IF ( .NOT. ALLOCATED( tmrt_grid ) ) THEN |
---|
1196 | ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) ) |
---|
1197 | ENDIF |
---|
1198 | tmrt_grid = REAL( bio_fill_value, KIND = wp ) |
---|
1199 | |
---|
1200 | DO l = 1, nmrtbl |
---|
1201 | i = mrtbl(ix,l) |
---|
1202 | j = mrtbl(iy,l) |
---|
1203 | k = mrtbl(iz,l) |
---|
1204 | IF ( k - get_topography_top_index_ji( j, i, 's' ) == bio_cell_level + & |
---|
1205 | 1_iwp) THEN |
---|
1206 | IF ( mrt_include_sw ) THEN |
---|
1207 | tmrt_grid(j,i) = ((human_absorb*mrtinsw(l) + & |
---|
1208 | human_emiss*mrtinlw(l)) / & |
---|
1209 | (human_emiss*sigma_sb)) ** .25_wp - degc_to_k |
---|
1210 | ELSE |
---|
1211 | tmrt_grid(j,i) = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp & |
---|
1212 | - degc_to_k |
---|
1213 | ENDIF |
---|
1214 | ENDIF |
---|
1215 | ENDDO |
---|
1216 | |
---|
1217 | END SUBROUTINE bio_calculate_mrt_grid |
---|
1218 | |
---|
1219 | |
---|
1220 | !------------------------------------------------------------------------------! |
---|
1221 | ! Description: |
---|
1222 | ! ------------ |
---|
1223 | !> Calculate static thermal indices for 2D grid point i, j |
---|
1224 | !------------------------------------------------------------------------------! |
---|
1225 | SUBROUTINE bio_get_thermal_index_input_ij( average_input, i, j, ta, vp, ws, & |
---|
1226 | pair, tmrt ) |
---|
1227 | |
---|
1228 | IMPLICIT NONE |
---|
1229 | ! |
---|
1230 | !-- Input variables |
---|
1231 | LOGICAL, INTENT ( IN ) :: average_input !< Determine averaged input conditions? |
---|
1232 | INTEGER(iwp), INTENT ( IN ) :: i !< Running index, x-dir |
---|
1233 | INTEGER(iwp), INTENT ( IN ) :: j !< Running index, y-dir |
---|
1234 | ! |
---|
1235 | !-- Output parameters |
---|
1236 | REAL(wp), INTENT ( OUT ) :: tmrt !< Mean radiant temperature (degree_C) |
---|
1237 | REAL(wp), INTENT ( OUT ) :: ta !< Air temperature (degree_C) |
---|
1238 | REAL(wp), INTENT ( OUT ) :: vp !< Vapour pressure (hPa) |
---|
1239 | REAL(wp), INTENT ( OUT ) :: ws !< Wind speed (local level) (m/s) |
---|
1240 | REAL(wp), INTENT ( OUT ) :: pair !< Air pressure (hPa) |
---|
1241 | ! |
---|
1242 | !-- Internal variables |
---|
1243 | INTEGER(iwp) :: k !< Running index, z-dir |
---|
1244 | ! INTEGER(iwp) :: ir !< Running index, x-dir, radiation coordinates |
---|
1245 | ! INTEGER(iwp) :: jr !< Running index, y-dir, radiation coordinates |
---|
1246 | ! INTEGER(iwp) :: kr !< Running index, y-dir, radiation coordinates |
---|
1247 | INTEGER(iwp) :: k_wind !< Running index, z-dir, wind speed only |
---|
1248 | ! INTEGER(iwp) :: l !< Running index, radiation coordinates |
---|
1249 | |
---|
1250 | REAL(wp) :: vp_sat !< Saturation vapor pressure (hPa) |
---|
1251 | |
---|
1252 | ! |
---|
1253 | !-- Determine cell level closest to 1.1m above ground |
---|
1254 | ! by making use of truncation due to int cast |
---|
1255 | k = get_topography_top_index_ji(j, i, 's') + bio_cell_level !< Vertical cell center closest to 1.1m |
---|
1256 | |
---|
1257 | ! |
---|
1258 | !-- Avoid non-representative horizontal u and v of 0.0 m/s too close to ground |
---|
1259 | k_wind = k |
---|
1260 | IF ( bio_cell_level < 1_iwp ) THEN |
---|
1261 | k_wind = k + 1_iwp |
---|
1262 | ENDIF |
---|
1263 | ! |
---|
1264 | !-- Determine local values: |
---|
1265 | IF ( average_input ) THEN |
---|
1266 | ! |
---|
1267 | !-- Calculate ta from Tp assuming dry adiabatic laps rate |
---|
1268 | ta = pt_av(k, j, i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k |
---|
1269 | |
---|
1270 | vp = bio_fill_value |
---|
1271 | IF ( humidity .AND. ALLOCATED( q_av ) ) THEN |
---|
1272 | vp = q_av(k, j, i) |
---|
1273 | ENDIF |
---|
1274 | |
---|
1275 | ws = ( 0.5_wp * ABS( u_av(k_wind, j, i) + u_av(k_wind, j, i+1) ) + & |
---|
1276 | 0.5_wp * ABS( v_av(k_wind, j, i) + v_av(k_wind, j+1, i) ) + & |
---|
1277 | 0.5_wp * ABS( w_av(k_wind, j, i) + w_av(k_wind+1, j, i) ) ) |
---|
1278 | ELSE |
---|
1279 | ! |
---|
1280 | !-- Calculate ta from Tp assuming dry adiabatic laps rate |
---|
1281 | ta = pt(k, j, i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k |
---|
1282 | |
---|
1283 | vp = bio_fill_value |
---|
1284 | IF ( humidity ) THEN |
---|
1285 | vp = q(k, j, i) |
---|
1286 | ENDIF |
---|
1287 | |
---|
1288 | ws = ( 0.5_wp * ABS( u(k_wind, j, i) + u(k_wind, j, i+1) ) + & |
---|
1289 | 0.5_wp * ABS( v(k_wind, j, i) + v(k_wind, j+1, i) ) + & |
---|
1290 | 0.5_wp * ABS( w(k_wind, j, i) + w(k_wind+1, j, i) ) ) |
---|
1291 | |
---|
1292 | ENDIF |
---|
1293 | ! |
---|
1294 | !-- Local air pressure |
---|
1295 | pair = surface_pressure |
---|
1296 | ! |
---|
1297 | !-- Calculate water vapour pressure at saturation and convert to hPa |
---|
1298 | !-- The magnus formula is limited to temperatures up to 333.15 K to |
---|
1299 | ! avoid negative values of vp_sat |
---|
1300 | IF ( vp > -998._wp ) THEN |
---|
1301 | vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) ) |
---|
1302 | vp = vp * pair / ( vp + 0.622_wp ) |
---|
1303 | IF ( vp > vp_sat ) vp = vp_sat |
---|
1304 | ENDIF |
---|
1305 | ! |
---|
1306 | !-- local mtr value at [i,j] |
---|
1307 | tmrt = bio_fill_value !< this can be a valid result (e.g. for inside some ostacle) |
---|
1308 | IF ( radiation ) THEN |
---|
1309 | ! |
---|
1310 | !-- Use MRT from RTM precalculated in tmrt_grid |
---|
1311 | tmrt = tmrt_grid(j,i) |
---|
1312 | ENDIF |
---|
1313 | |
---|
1314 | END SUBROUTINE bio_get_thermal_index_input_ij |
---|
1315 | |
---|
1316 | |
---|
1317 | !------------------------------------------------------------------------------! |
---|
1318 | ! Description: |
---|
1319 | ! ------------ |
---|
1320 | !> Calculate static thermal indices for any point within a 2D grid |
---|
1321 | !> time_integration.f90: 1065ff |
---|
1322 | !------------------------------------------------------------------------------! |
---|
1323 | SUBROUTINE bio_calculate_thermal_index_maps ( av ) |
---|
1324 | |
---|
1325 | IMPLICIT NONE |
---|
1326 | ! |
---|
1327 | !-- Input attributes |
---|
1328 | LOGICAL, INTENT ( IN ) :: av !< Calculate based on averaged input conditions? |
---|
1329 | ! |
---|
1330 | !-- Internal variables |
---|
1331 | INTEGER(iwp) :: i, j !< Running index |
---|
1332 | |
---|
1333 | REAL(wp) :: clo !< Clothing index (no dimension) |
---|
1334 | REAL(wp) :: ta !< Air temperature (degree_C) |
---|
1335 | REAL(wp) :: vp !< Vapour pressure (hPa) |
---|
1336 | REAL(wp) :: ws !< Wind speed (local level) (m/s) |
---|
1337 | REAL(wp) :: pair !< Air pressure (hPa) |
---|
1338 | REAL(wp) :: perct_ij !< Perceived temperature (degree_C) |
---|
1339 | REAL(wp) :: utci_ij !< Universal thermal climate index (degree_C) |
---|
1340 | REAL(wp) :: pet_ij !< Physiologically equivalent temperature (degree_C) |
---|
1341 | REAL(wp) :: tmrt_ij !< Mean radiant temperature (degree_C) |
---|
1342 | |
---|
1343 | ! |
---|
1344 | !-- fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...) |
---|
1345 | IF ( simulated_time > 0.0_wp ) THEN |
---|
1346 | CALL bio_calculate_mrt_grid ( av ) |
---|
1347 | ENDIF |
---|
1348 | |
---|
1349 | DO i = nxl, nxr |
---|
1350 | DO j = nys, nyn |
---|
1351 | ! |
---|
1352 | !-- Determine local input conditions |
---|
1353 | tmrt_ij = bio_fill_value |
---|
1354 | vp = bio_fill_value |
---|
1355 | ! |
---|
1356 | !-- Determine input only if |
---|
1357 | IF ( simulated_time > 0.0_wp ) THEN |
---|
1358 | CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp, & |
---|
1359 | ws, pair, tmrt_ij ) |
---|
1360 | END IF |
---|
1361 | ! |
---|
1362 | !-- Only proceed if input is available |
---|
1363 | pet_ij = bio_fill_value !< set fail value, e.g. valid for within |
---|
1364 | perct_ij = bio_fill_value !< some obstacle |
---|
1365 | utci_ij = bio_fill_value |
---|
1366 | IF ( .NOT. ( tmrt_ij <= -998._wp .OR. vp <= -998._wp ) ) THEN |
---|
1367 | ! |
---|
1368 | !-- Calculate static thermal indices based on local tmrt |
---|
1369 | clo = bio_fill_value |
---|
1370 | |
---|
1371 | IF ( bio_perct ) THEN |
---|
1372 | ! |
---|
1373 | !-- Estimate local perceived temperature |
---|
1374 | CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair, clo, & |
---|
1375 | perct_ij ) |
---|
1376 | ENDIF |
---|
1377 | |
---|
1378 | IF ( bio_utci ) THEN |
---|
1379 | ! |
---|
1380 | !-- Estimate local universal thermal climate index |
---|
1381 | CALL calculate_utci_static( ta, vp, ws, tmrt_ij, & |
---|
1382 | bio_output_height, utci_ij ) |
---|
1383 | ENDIF |
---|
1384 | |
---|
1385 | IF ( bio_pet ) THEN |
---|
1386 | ! |
---|
1387 | !-- Estimate local physiologically equivalent temperature |
---|
1388 | CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair, pet_ij ) |
---|
1389 | ENDIF |
---|
1390 | END IF |
---|
1391 | |
---|
1392 | |
---|
1393 | IF ( av ) THEN |
---|
1394 | ! |
---|
1395 | !-- Write results for selected averaged indices |
---|
1396 | IF ( bio_perct_av ) THEN |
---|
1397 | perct_av(j, i) = perct_ij |
---|
1398 | END IF |
---|
1399 | IF ( bio_utci_av ) THEN |
---|
1400 | utci_av(j, i) = utci_ij |
---|
1401 | END IF |
---|
1402 | IF ( bio_pet_av ) THEN |
---|
1403 | pet_av(j, i) = pet_ij |
---|
1404 | END IF |
---|
1405 | ELSE |
---|
1406 | ! |
---|
1407 | !-- Write result for selected indices |
---|
1408 | IF ( bio_perct ) THEN |
---|
1409 | perct(j, i) = perct_ij |
---|
1410 | END IF |
---|
1411 | IF ( bio_utci ) THEN |
---|
1412 | utci(j, i) = utci_ij |
---|
1413 | END IF |
---|
1414 | IF ( bio_pet ) THEN |
---|
1415 | pet(j, i) = pet_ij |
---|
1416 | END IF |
---|
1417 | END IF |
---|
1418 | |
---|
1419 | END DO |
---|
1420 | END DO |
---|
1421 | |
---|
1422 | END SUBROUTINE bio_calculate_thermal_index_maps |
---|
1423 | |
---|
1424 | !------------------------------------------------------------------------------! |
---|
1425 | ! Description: |
---|
1426 | ! ------------ |
---|
1427 | !> Calculate dynamic thermal indices (currently only iPT, but expandable) |
---|
1428 | !------------------------------------------------------------------------------! |
---|
1429 | SUBROUTINE bio_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage, & |
---|
1430 | t_clo, clo, actlev, age, weight, height, work, sex, ipt ) |
---|
1431 | |
---|
1432 | IMPLICIT NONE |
---|
1433 | ! |
---|
1434 | !-- Input parameters |
---|
1435 | REAL(wp), INTENT ( IN ) :: ta !< Air temperature (degree_C) |
---|
1436 | REAL(wp), INTENT ( IN ) :: vp !< Vapour pressure (hPa) |
---|
1437 | REAL(wp), INTENT ( IN ) :: ws !< Wind speed (local level) (m/s) |
---|
1438 | REAL(wp), INTENT ( IN ) :: pair !< Air pressure (hPa) |
---|
1439 | REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (degree_C) |
---|
1440 | REAL(wp), INTENT ( IN ) :: dt !< Time past since last calculation (s) |
---|
1441 | REAL(wp), INTENT ( IN ) :: age !< Age of agent (y) |
---|
1442 | REAL(wp), INTENT ( IN ) :: weight !< Weight of agent (Kg) |
---|
1443 | REAL(wp), INTENT ( IN ) :: height !< Height of agent (m) |
---|
1444 | REAL(wp), INTENT ( IN ) :: work !< Mechanical workload of agent |
---|
1445 | ! (without metabolism!) (W) |
---|
1446 | INTEGER(iwp), INTENT ( IN ) :: sex !< Sex of agent (1 = male, 2 = female) |
---|
1447 | ! |
---|
1448 | !-- Both, input and output parameters |
---|
1449 | Real(wp), INTENT ( INOUT ) :: energy_storage !< Energy storage (W/m²) |
---|
1450 | Real(wp), INTENT ( INOUT ) :: t_clo !< Clothing temperature (degree_C) |
---|
1451 | Real(wp), INTENT ( INOUT ) :: clo !< Current clothing in sulation |
---|
1452 | Real(wp), INTENT ( INOUT ) :: actlev !< Individuals activity level |
---|
1453 | ! per unit surface area (W/m²) |
---|
1454 | ! |
---|
1455 | !-- Output parameters |
---|
1456 | REAL(wp), INTENT ( OUT ) :: ipt !< Instationary perceived temp. (degree_C) |
---|
1457 | ! |
---|
1458 | !-- If clo equals the initial value, this is the initial call |
---|
1459 | IF ( clo <= -998._wp ) THEN |
---|
1460 | ! |
---|
1461 | !-- Initialize instationary perceived temperature with personalized |
---|
1462 | ! PT as an initial guess, set actlev and clo |
---|
1463 | CALL ipt_init ( age, weight, height, sex, work, actlev, clo, & |
---|
1464 | ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo, & |
---|
1465 | ipt ) |
---|
1466 | ELSE |
---|
1467 | ! |
---|
1468 | !-- Estimate local instatinoary perceived temperature |
---|
1469 | CALL ipt_cycle ( ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo, & |
---|
1470 | clo, actlev, work, ipt ) |
---|
1471 | ENDIF |
---|
1472 | |
---|
1473 | END SUBROUTINE bio_calc_ipt |
---|
1474 | |
---|
1475 | |
---|
1476 | !------------------------------------------------------------------------------! |
---|
1477 | ! Description: |
---|
1478 | ! ------------ |
---|
1479 | !> SUBROUTINE for calculating UTCI Temperature (UTCI) |
---|
1480 | !> computed by a 6th order approximation |
---|
1481 | ! |
---|
1482 | !> UTCI regression equation after |
---|
1483 | !> Bröde P, Fiala D, Blazejczyk K, Holmér I, Jendritzky G, Kampmann B, Tinz B, |
---|
1484 | !> Havenith G (2012) Deriving the operational procedure for the Universal Thermal |
---|
1485 | !> Climate Index (UTCI). International Journal of Biometeorology 56 (3):481-494. |
---|
1486 | !> doi:10.1007/s00484-011-0454-1 |
---|
1487 | ! |
---|
1488 | !> original source available at: |
---|
1489 | !> www.utci.org |
---|
1490 | !------------------------------------------------------------------------------! |
---|
1491 | SUBROUTINE calculate_utci_static( ta_in, vp, ws_hag, tmrt, hag, utci_ij ) |
---|
1492 | |
---|
1493 | IMPLICIT NONE |
---|
1494 | ! |
---|
1495 | !-- Type of input of the argument list |
---|
1496 | REAL(WP), INTENT ( IN ) :: ta_in !< Local air temperature (degree_C) |
---|
1497 | REAL(WP), INTENT ( IN ) :: vp !< Loacl vapour pressure (hPa) |
---|
1498 | REAL(WP), INTENT ( IN ) :: ws_hag !< Incident wind speed (m/s) |
---|
1499 | REAL(WP), INTENT ( IN ) :: tmrt !< Local mean radiant temperature (degree_C) |
---|
1500 | REAL(WP), INTENT ( IN ) :: hag !< Height of wind speed input (m) |
---|
1501 | ! |
---|
1502 | !-- Type of output of the argument list |
---|
1503 | REAL(wp), INTENT ( OUT ) :: utci_ij !< Universal Thermal Climate Index (degree_C) |
---|
1504 | |
---|
1505 | REAL(WP) :: ta !< air temperature modified by offset (degree_C) |
---|
1506 | REAL(WP) :: pa !< air pressure in kPa (kPa) |
---|
1507 | REAL(WP) :: d_tmrt !< delta-tmrt (degree_C) |
---|
1508 | REAL(WP) :: va !< wind speed at 10 m above ground level (m/s) |
---|
1509 | REAL(WP) :: offset !< utci deviation by ta cond. exceeded (degree_C) |
---|
1510 | REAL(WP) :: part_ta !< Air temperature related part of the regression |
---|
1511 | REAL(WP) :: ta2 !< 2 times ta |
---|
1512 | REAL(WP) :: ta3 !< 3 times ta |
---|
1513 | REAL(WP) :: ta4 !< 4 times ta |
---|
1514 | REAL(WP) :: ta5 !< 5 times ta |
---|
1515 | REAL(WP) :: ta6 !< 6 times ta |
---|
1516 | REAL(WP) :: part_va !< Vapour pressure related part of the regression |
---|
1517 | REAL(WP) :: va2 !< 2 times va |
---|
1518 | REAL(WP) :: va3 !< 3 times va |
---|
1519 | REAL(WP) :: va4 !< 4 times va |
---|
1520 | REAL(WP) :: va5 !< 5 times va |
---|
1521 | REAL(WP) :: va6 !< 6 times va |
---|
1522 | REAL(WP) :: part_d_tmrt !< Mean radiant temp. related part of the reg. |
---|
1523 | REAL(WP) :: d_tmrt2 !< 2 times d_tmrt |
---|
1524 | REAL(WP) :: d_tmrt3 !< 3 times d_tmrt |
---|
1525 | REAL(WP) :: d_tmrt4 !< 4 times d_tmrt |
---|
1526 | REAL(WP) :: d_tmrt5 !< 5 times d_tmrt |
---|
1527 | REAL(WP) :: d_tmrt6 !< 6 times d_tmrt |
---|
1528 | REAL(WP) :: part_pa !< Air pressure related part of the regression |
---|
1529 | REAL(WP) :: pa2 !< 2 times pa |
---|
1530 | REAL(WP) :: pa3 !< 3 times pa |
---|
1531 | REAL(WP) :: pa4 !< 4 times pa |
---|
1532 | REAL(WP) :: pa5 !< 5 times pa |
---|
1533 | REAL(WP) :: pa6 !< 6 times pa |
---|
1534 | REAL(WP) :: part_pa2 !< Air pressure^2 related part of the regression |
---|
1535 | REAL(WP) :: part_pa3 !< Air pressure^3 related part of the regression |
---|
1536 | REAL(WP) :: part_pa46 !< Air pressure^4-6 related part of the regression |
---|
1537 | |
---|
1538 | ! |
---|
1539 | !-- Initialize |
---|
1540 | offset = 0._wp |
---|
1541 | ta = ta_in |
---|
1542 | d_tmrt = tmrt - ta_in |
---|
1543 | ! |
---|
1544 | !-- Use vapour pressure in kpa |
---|
1545 | pa = vp / 10.0_wp |
---|
1546 | ! |
---|
1547 | !-- Wind altitude correction from hag to 10m after Broede et al. (2012), eq.3 |
---|
1548 | ! z(0) is set to 0.01 according to UTCI profile definition |
---|
1549 | va = ws_hag * log ( 10.0_wp / 0.01_wp ) / log ( hag / 0.01_wp ) |
---|
1550 | ! |
---|
1551 | !-- Check if input values in range after Broede et al. (2012) |
---|
1552 | IF ( ( d_tmrt > 70._wp ) .OR. ( d_tmrt < -30._wp ) .OR. & |
---|
1553 | ( vp >= 50._wp ) ) THEN |
---|
1554 | utci_ij = bio_fill_value |
---|
1555 | RETURN |
---|
1556 | ENDIF |
---|
1557 | ! |
---|
1558 | !-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds |
---|
1559 | IF ( ta > 50._wp ) THEN |
---|
1560 | offset = ta - 50._wp |
---|
1561 | ta = 50._wp |
---|
1562 | ENDIF |
---|
1563 | IF ( ta < -50._wp ) THEN |
---|
1564 | offset = ta + 50._wp |
---|
1565 | ta = -50._wp |
---|
1566 | ENDIF |
---|
1567 | ! |
---|
1568 | !-- For routine application. For wind speeds and relative |
---|
1569 | ! humidity values below 0.5 m/s or 5%, respectively, the |
---|
1570 | ! user is advised to use the lower bounds for the calculations. |
---|
1571 | IF ( va < 0.5_wp ) va = 0.5_wp |
---|
1572 | IF ( va > 17._wp ) va = 17._wp |
---|
1573 | |
---|
1574 | ! |
---|
1575 | !-- Pre-calculate multiples of input parameters to save time later |
---|
1576 | |
---|
1577 | ta2 = ta * ta |
---|
1578 | ta3 = ta2 * ta |
---|
1579 | ta4 = ta3 * ta |
---|
1580 | ta5 = ta4 * ta |
---|
1581 | ta6 = ta5 * ta |
---|
1582 | |
---|
1583 | va2 = va * va |
---|
1584 | va3 = va2 * va |
---|
1585 | va4 = va3 * va |
---|
1586 | va5 = va4 * va |
---|
1587 | va6 = va5 * va |
---|
1588 | |
---|
1589 | d_tmrt2 = d_tmrt * d_tmrt |
---|
1590 | d_tmrt3 = d_tmrt2 * d_tmrt |
---|
1591 | d_tmrt4 = d_tmrt3 * d_tmrt |
---|
1592 | d_tmrt5 = d_tmrt4 * d_tmrt |
---|
1593 | d_tmrt6 = d_tmrt5 * d_tmrt |
---|
1594 | |
---|
1595 | pa2 = pa * pa |
---|
1596 | pa3 = pa2 * pa |
---|
1597 | pa4 = pa3 * pa |
---|
1598 | pa5 = pa4 * pa |
---|
1599 | pa6 = pa5 * pa |
---|
1600 | |
---|
1601 | ! |
---|
1602 | !-- Pre-calculate parts of the regression equation |
---|
1603 | part_ta = ( 6.07562052e-01_wp ) + & |
---|
1604 | ( -2.27712343e-02_wp ) * ta + & |
---|
1605 | ( 8.06470249e-04_wp ) * ta2 + & |
---|
1606 | ( -1.54271372e-04_wp ) * ta3 + & |
---|
1607 | ( -3.24651735e-06_wp ) * ta4 + & |
---|
1608 | ( 7.32602852e-08_wp ) * ta5 + & |
---|
1609 | ( 1.35959073e-09_wp ) * ta6 |
---|
1610 | |
---|
1611 | part_va = ( -2.25836520e+00_wp ) * va + & |
---|
1612 | ( 8.80326035e-02_wp ) * ta * va + & |
---|
1613 | ( 2.16844454e-03_wp ) * ta2 * va + & |
---|
1614 | ( -1.53347087e-05_wp ) * ta3 * va + & |
---|
1615 | ( -5.72983704e-07_wp ) * ta4 * va + & |
---|
1616 | ( -2.55090145e-09_wp ) * ta5 * va + & |
---|
1617 | ( -7.51269505e-01_wp ) * va2 + & |
---|
1618 | ( -4.08350271e-03_wp ) * ta * va2 + & |
---|
1619 | ( -5.21670675e-05_wp ) * ta2 * va2 + & |
---|
1620 | ( 1.94544667e-06_wp ) * ta3 * va2 + & |
---|
1621 | ( 1.14099531e-08_wp ) * ta4 * va2 + & |
---|
1622 | ( 1.58137256e-01_wp ) * va3 + & |
---|
1623 | ( -6.57263143e-05_wp ) * ta * va3 + & |
---|
1624 | ( 2.22697524e-07_wp ) * ta2 * va3 + & |
---|
1625 | ( -4.16117031e-08_wp ) * ta3 * va3 + & |
---|
1626 | ( -1.27762753e-02_wp ) * va4 + & |
---|
1627 | ( 9.66891875e-06_wp ) * ta * va4 + & |
---|
1628 | ( 2.52785852e-09_wp ) * ta2 * va4 + & |
---|
1629 | ( 4.56306672e-04_wp ) * va5 + & |
---|
1630 | ( -1.74202546e-07_wp ) * ta * va5 + & |
---|
1631 | ( -5.91491269e-06_wp ) * va6 |
---|
1632 | |
---|
1633 | part_d_tmrt = ( 3.98374029e-01_wp ) * d_tmrt + & |
---|
1634 | ( 1.83945314e-04_wp ) * ta * d_tmrt + & |
---|
1635 | ( -1.73754510e-04_wp ) * ta2 * d_tmrt + & |
---|
1636 | ( -7.60781159e-07_wp ) * ta3 * d_tmrt + & |
---|
1637 | ( 3.77830287e-08_wp ) * ta4 * d_tmrt + & |
---|
1638 | ( 5.43079673e-10_wp ) * ta5 * d_tmrt + & |
---|
1639 | ( -2.00518269e-02_wp ) * va * d_tmrt + & |
---|
1640 | ( 8.92859837e-04_wp ) * ta * va * d_tmrt + & |
---|
1641 | ( 3.45433048e-06_wp ) * ta2 * va * d_tmrt + & |
---|
1642 | ( -3.77925774e-07_wp ) * ta3 * va * d_tmrt + & |
---|
1643 | ( -1.69699377e-09_wp ) * ta4 * va * d_tmrt + & |
---|
1644 | ( 1.69992415e-04_wp ) * va2 * d_tmrt + & |
---|
1645 | ( -4.99204314e-05_wp ) * ta * va2 * d_tmrt + & |
---|
1646 | ( 2.47417178e-07_wp ) * ta2 * va2 * d_tmrt + & |
---|
1647 | ( 1.07596466e-08_wp ) * ta3 * va2 * d_tmrt + & |
---|
1648 | ( 8.49242932e-05_wp ) * va3 * d_tmrt + & |
---|
1649 | ( 1.35191328e-06_wp ) * ta * va3 * d_tmrt + & |
---|
1650 | ( -6.21531254e-09_wp ) * ta2 * va3 * d_tmrt + & |
---|
1651 | ( -4.99410301e-06_wp ) * va4 * d_tmrt + & |
---|
1652 | ( -1.89489258e-08_wp ) * ta * va4 * d_tmrt + & |
---|
1653 | ( 8.15300114e-08_wp ) * va5 * d_tmrt + & |
---|
1654 | ( 7.55043090e-04_wp ) * d_tmrt2 + & |
---|
1655 | ( -5.65095215e-05_wp ) * ta * d_tmrt2 + & |
---|
1656 | ( -4.52166564e-07_wp ) * ta2 * d_tmrt2 + & |
---|
1657 | ( 2.46688878e-08_wp ) * ta3 * d_tmrt2 + & |
---|
1658 | ( 2.42674348e-10_wp ) * ta4 * d_tmrt2 + & |
---|
1659 | ( 1.54547250e-04_wp ) * va * d_tmrt2 + & |
---|
1660 | ( 5.24110970e-06_wp ) * ta * va * d_tmrt2 + & |
---|
1661 | ( -8.75874982e-08_wp ) * ta2 * va * d_tmrt2 + & |
---|
1662 | ( -1.50743064e-09_wp ) * ta3 * va * d_tmrt2 + & |
---|
1663 | ( -1.56236307e-05_wp ) * va2 * d_tmrt2 + & |
---|
1664 | ( -1.33895614e-07_wp ) * ta * va2 * d_tmrt2 + & |
---|
1665 | ( 2.49709824e-09_wp ) * ta2 * va2 * d_tmrt2 + & |
---|
1666 | ( 6.51711721e-07_wp ) * va3 * d_tmrt2 + & |
---|
1667 | ( 1.94960053e-09_wp ) * ta * va3 * d_tmrt2 + & |
---|
1668 | ( -1.00361113e-08_wp ) * va4 * d_tmrt2 + & |
---|
1669 | ( -1.21206673e-05_wp ) * d_tmrt3 + & |
---|
1670 | ( -2.18203660e-07_wp ) * ta * d_tmrt3 + & |
---|
1671 | ( 7.51269482e-09_wp ) * ta2 * d_tmrt3 + & |
---|
1672 | ( 9.79063848e-11_wp ) * ta3 * d_tmrt3 + & |
---|
1673 | ( 1.25006734e-06_wp ) * va * d_tmrt3 + & |
---|
1674 | ( -1.81584736e-09_wp ) * ta * va * d_tmrt3 + & |
---|
1675 | ( -3.52197671e-10_wp ) * ta2 * va * d_tmrt3 + & |
---|
1676 | ( -3.36514630e-08_wp ) * va2 * d_tmrt3 + & |
---|
1677 | ( 1.35908359e-10_wp ) * ta * va2 * d_tmrt3 + & |
---|
1678 | ( 4.17032620e-10_wp ) * va3 * d_tmrt3 + & |
---|
1679 | ( -1.30369025e-09_wp ) * d_tmrt4 + & |
---|
1680 | ( 4.13908461e-10_wp ) * ta * d_tmrt4 + & |
---|
1681 | ( 9.22652254e-12_wp ) * ta2 * d_tmrt4 + & |
---|
1682 | ( -5.08220384e-09_wp ) * va * d_tmrt4 + & |
---|
1683 | ( -2.24730961e-11_wp ) * ta * va * d_tmrt4 + & |
---|
1684 | ( 1.17139133e-10_wp ) * va2 * d_tmrt4 + & |
---|
1685 | ( 6.62154879e-10_wp ) * d_tmrt5 + & |
---|
1686 | ( 4.03863260e-13_wp ) * ta * d_tmrt5 + & |
---|
1687 | ( 1.95087203e-12_wp ) * va * d_tmrt5 + & |
---|
1688 | ( -4.73602469e-12_wp ) * d_tmrt6 |
---|
1689 | |
---|
1690 | part_pa = ( 5.12733497e+00_wp ) * pa + & |
---|
1691 | ( -3.12788561e-01_wp ) * ta * pa + & |
---|
1692 | ( -1.96701861e-02_wp ) * ta2 * pa + & |
---|
1693 | ( 9.99690870e-04_wp ) * ta3 * pa + & |
---|
1694 | ( 9.51738512e-06_wp ) * ta4 * pa + & |
---|
1695 | ( -4.66426341e-07_wp ) * ta5 * pa + & |
---|
1696 | ( 5.48050612e-01_wp ) * va * pa + & |
---|
1697 | ( -3.30552823e-03_wp ) * ta * va * pa + & |
---|
1698 | ( -1.64119440e-03_wp ) * ta2 * va * pa + & |
---|
1699 | ( -5.16670694e-06_wp ) * ta3 * va * pa + & |
---|
1700 | ( 9.52692432e-07_wp ) * ta4 * va * pa + & |
---|
1701 | ( -4.29223622e-02_wp ) * va2 * pa + & |
---|
1702 | ( 5.00845667e-03_wp ) * ta * va2 * pa + & |
---|
1703 | ( 1.00601257e-06_wp ) * ta2 * va2 * pa + & |
---|
1704 | ( -1.81748644e-06_wp ) * ta3 * va2 * pa + & |
---|
1705 | ( -1.25813502e-03_wp ) * va3 * pa + & |
---|
1706 | ( -1.79330391e-04_wp ) * ta * va3 * pa + & |
---|
1707 | ( 2.34994441e-06_wp ) * ta2 * va3 * pa + & |
---|
1708 | ( 1.29735808e-04_wp ) * va4 * pa + & |
---|
1709 | ( 1.29064870e-06_wp ) * ta * va4 * pa + & |
---|
1710 | ( -2.28558686e-06_wp ) * va5 * pa + & |
---|
1711 | ( -3.69476348e-02_wp ) * d_tmrt * pa + & |
---|
1712 | ( 1.62325322e-03_wp ) * ta * d_tmrt * pa + & |
---|
1713 | ( -3.14279680e-05_wp ) * ta2 * d_tmrt * pa + & |
---|
1714 | ( 2.59835559e-06_wp ) * ta3 * d_tmrt * pa + & |
---|
1715 | ( -4.77136523e-08_wp ) * ta4 * d_tmrt * pa + & |
---|
1716 | ( 8.64203390e-03_wp ) * va * d_tmrt * pa + & |
---|
1717 | ( -6.87405181e-04_wp ) * ta * va * d_tmrt * pa + & |
---|
1718 | ( -9.13863872e-06_wp ) * ta2 * va * d_tmrt * pa + & |
---|
1719 | ( 5.15916806e-07_wp ) * ta3 * va * d_tmrt * pa + & |
---|
1720 | ( -3.59217476e-05_wp ) * va2 * d_tmrt * pa + & |
---|
1721 | ( 3.28696511e-05_wp ) * ta * va2 * d_tmrt * pa + & |
---|
1722 | ( -7.10542454e-07_wp ) * ta2 * va2 * d_tmrt * pa + & |
---|
1723 | ( -1.24382300e-05_wp ) * va3 * d_tmrt * pa + & |
---|
1724 | ( -7.38584400e-09_wp ) * ta * va3 * d_tmrt * pa + & |
---|
1725 | ( 2.20609296e-07_wp ) * va4 * d_tmrt * pa + & |
---|
1726 | ( -7.32469180e-04_wp ) * d_tmrt2 * pa + & |
---|
1727 | ( -1.87381964e-05_wp ) * ta * d_tmrt2 * pa + & |
---|
1728 | ( 4.80925239e-06_wp ) * ta2 * d_tmrt2 * pa + & |
---|
1729 | ( -8.75492040e-08_wp ) * ta3 * d_tmrt2 * pa + & |
---|
1730 | ( 2.77862930e-05_wp ) * va * d_tmrt2 * pa + & |
---|
1731 | ( -5.06004592e-06_wp ) * ta * va * d_tmrt2 * pa + & |
---|
1732 | ( 1.14325367e-07_wp ) * ta2 * va * d_tmrt2 * pa + & |
---|
1733 | ( 2.53016723e-06_wp ) * va2 * d_tmrt2 * pa + & |
---|
1734 | ( -1.72857035e-08_wp ) * ta * va2 * d_tmrt2 * pa + & |
---|
1735 | ( -3.95079398e-08_wp ) * va3 * d_tmrt2 * pa + & |
---|
1736 | ( -3.59413173e-07_wp ) * d_tmrt3 * pa + & |
---|
1737 | ( 7.04388046e-07_wp ) * ta * d_tmrt3 * pa + & |
---|
1738 | ( -1.89309167e-08_wp ) * ta2 * d_tmrt3 * pa + & |
---|
1739 | ( -4.79768731e-07_wp ) * va * d_tmrt3 * pa + & |
---|
1740 | ( 7.96079978e-09_wp ) * ta * va * d_tmrt3 * pa + & |
---|
1741 | ( 1.62897058e-09_wp ) * va2 * d_tmrt3 * pa + & |
---|
1742 | ( 3.94367674e-08_wp ) * d_tmrt4 * pa + & |
---|
1743 | ( -1.18566247e-09_wp ) * ta * d_tmrt4 * pa + & |
---|
1744 | ( 3.34678041e-10_wp ) * va * d_tmrt4 * pa + & |
---|
1745 | ( -1.15606447e-10_wp ) * d_tmrt5 * pa |
---|
1746 | |
---|
1747 | part_pa2 = ( -2.80626406e+00_wp ) * pa2 + & |
---|
1748 | ( 5.48712484e-01_wp ) * ta * pa2 + & |
---|
1749 | ( -3.99428410e-03_wp ) * ta2 * pa2 + & |
---|
1750 | ( -9.54009191e-04_wp ) * ta3 * pa2 + & |
---|
1751 | ( 1.93090978e-05_wp ) * ta4 * pa2 + & |
---|
1752 | ( -3.08806365e-01_wp ) * va * pa2 + & |
---|
1753 | ( 1.16952364e-02_wp ) * ta * va * pa2 + & |
---|
1754 | ( 4.95271903e-04_wp ) * ta2 * va * pa2 + & |
---|
1755 | ( -1.90710882e-05_wp ) * ta3 * va * pa2 + & |
---|
1756 | ( 2.10787756e-03_wp ) * va2 * pa2 + & |
---|
1757 | ( -6.98445738e-04_wp ) * ta * va2 * pa2 + & |
---|
1758 | ( 2.30109073e-05_wp ) * ta2 * va2 * pa2 + & |
---|
1759 | ( 4.17856590e-04_wp ) * va3 * pa2 + & |
---|
1760 | ( -1.27043871e-05_wp ) * ta * va3 * pa2 + & |
---|
1761 | ( -3.04620472e-06_wp ) * va4 * pa2 + & |
---|
1762 | ( 5.14507424e-02_wp ) * d_tmrt * pa2 + & |
---|
1763 | ( -4.32510997e-03_wp ) * ta * d_tmrt * pa2 + & |
---|
1764 | ( 8.99281156e-05_wp ) * ta2 * d_tmrt * pa2 + & |
---|
1765 | ( -7.14663943e-07_wp ) * ta3 * d_tmrt * pa2 + & |
---|
1766 | ( -2.66016305e-04_wp ) * va * d_tmrt * pa2 + & |
---|
1767 | ( 2.63789586e-04_wp ) * ta * va * d_tmrt * pa2 + & |
---|
1768 | ( -7.01199003e-06_wp ) * ta2 * va * d_tmrt * pa2 + & |
---|
1769 | ( -1.06823306e-04_wp ) * va2 * d_tmrt * pa2 + & |
---|
1770 | ( 3.61341136e-06_wp ) * ta * va2 * d_tmrt * pa2 + & |
---|
1771 | ( 2.29748967e-07_wp ) * va3 * d_tmrt * pa2 + & |
---|
1772 | ( 3.04788893e-04_wp ) * d_tmrt2 * pa2 + & |
---|
1773 | ( -6.42070836e-05_wp ) * ta * d_tmrt2 * pa2 + & |
---|
1774 | ( 1.16257971e-06_wp ) * ta2 * d_tmrt2 * pa2 + & |
---|
1775 | ( 7.68023384e-06_wp ) * va * d_tmrt2 * pa2 + & |
---|
1776 | ( -5.47446896e-07_wp ) * ta * va * d_tmrt2 * pa2 + & |
---|
1777 | ( -3.59937910e-08_wp ) * va2 * d_tmrt2 * pa2 + & |
---|
1778 | ( -4.36497725e-06_wp ) * d_tmrt3 * pa2 + & |
---|
1779 | ( 1.68737969e-07_wp ) * ta * d_tmrt3 * pa2 + & |
---|
1780 | ( 2.67489271e-08_wp ) * va * d_tmrt3 * pa2 + & |
---|
1781 | ( 3.23926897e-09_wp ) * d_tmrt4 * pa2 |
---|
1782 | |
---|
1783 | part_pa3 = ( -3.53874123e-02_wp ) * pa3 + & |
---|
1784 | ( -2.21201190e-01_wp ) * ta * pa3 + & |
---|
1785 | ( 1.55126038e-02_wp ) * ta2 * pa3 + & |
---|
1786 | ( -2.63917279e-04_wp ) * ta3 * pa3 + & |
---|
1787 | ( 4.53433455e-02_wp ) * va * pa3 + & |
---|
1788 | ( -4.32943862e-03_wp ) * ta * va * pa3 + & |
---|
1789 | ( 1.45389826e-04_wp ) * ta2 * va * pa3 + & |
---|
1790 | ( 2.17508610e-04_wp ) * va2 * pa3 + & |
---|
1791 | ( -6.66724702e-05_wp ) * ta * va2 * pa3 + & |
---|
1792 | ( 3.33217140e-05_wp ) * va3 * pa3 + & |
---|
1793 | ( -2.26921615e-03_wp ) * d_tmrt * pa3 + & |
---|
1794 | ( 3.80261982e-04_wp ) * ta * d_tmrt * pa3 + & |
---|
1795 | ( -5.45314314e-09_wp ) * ta2 * d_tmrt * pa3 + & |
---|
1796 | ( -7.96355448e-04_wp ) * va * d_tmrt * pa3 + & |
---|
1797 | ( 2.53458034e-05_wp ) * ta * va * d_tmrt * pa3 + & |
---|
1798 | ( -6.31223658e-06_wp ) * va2 * d_tmrt * pa3 + & |
---|
1799 | ( 3.02122035e-04_wp ) * d_tmrt2 * pa3 + & |
---|
1800 | ( -4.77403547e-06_wp ) * ta * d_tmrt2 * pa3 + & |
---|
1801 | ( 1.73825715e-06_wp ) * va * d_tmrt2 * pa3 + & |
---|
1802 | ( -4.09087898e-07_wp ) * d_tmrt3 * pa3 |
---|
1803 | |
---|
1804 | part_pa46 = ( 6.14155345e-01_wp ) * pa4 + & |
---|
1805 | ( -6.16755931e-02_wp ) * ta * pa4 + & |
---|
1806 | ( 1.33374846e-03_wp ) * ta2 * pa4 + & |
---|
1807 | ( 3.55375387e-03_wp ) * va * pa4 + & |
---|
1808 | ( -5.13027851e-04_wp ) * ta * va * pa4 + & |
---|
1809 | ( 1.02449757e-04_wp ) * va2 * pa4 + & |
---|
1810 | ( -1.48526421e-03_wp ) * d_tmrt * pa4 + & |
---|
1811 | ( -4.11469183e-05_wp ) * ta * d_tmrt * pa4 + & |
---|
1812 | ( -6.80434415e-06_wp ) * va * d_tmrt * pa4 + & |
---|
1813 | ( -9.77675906e-06_wp ) * d_tmrt2 * pa4 + & |
---|
1814 | ( 8.82773108e-02_wp ) * pa5 + & |
---|
1815 | ( -3.01859306e-03_wp ) * ta * pa5 + & |
---|
1816 | ( 1.04452989e-03_wp ) * va * pa5 + & |
---|
1817 | ( 2.47090539e-04_wp ) * d_tmrt * pa5 + & |
---|
1818 | ( 1.48348065e-03_wp ) * pa6 |
---|
1819 | ! |
---|
1820 | !-- Calculate 6th order polynomial as approximation |
---|
1821 | utci_ij = ta + part_ta + part_va + part_d_tmrt + part_pa + part_pa2 + & |
---|
1822 | part_pa3 + part_pa46 |
---|
1823 | ! |
---|
1824 | !-- Consider offset in result |
---|
1825 | utci_ij = utci_ij + offset |
---|
1826 | |
---|
1827 | END SUBROUTINE calculate_utci_static |
---|
1828 | |
---|
1829 | |
---|
1830 | |
---|
1831 | |
---|
1832 | !------------------------------------------------------------------------------! |
---|
1833 | ! Description: |
---|
1834 | ! ------------ |
---|
1835 | !> calculate_perct_static: Estimation of perceived temperature (PT, degree_C) |
---|
1836 | !> Value of perct is the Perceived Temperature, degree centigrade |
---|
1837 | !------------------------------------------------------------------------------! |
---|
1838 | SUBROUTINE calculate_perct_static( ta, vp, ws, tmrt, pair, clo, perct_ij ) |
---|
1839 | |
---|
1840 | IMPLICIT NONE |
---|
1841 | ! |
---|
1842 | !-- Type of input of the argument list |
---|
1843 | REAL(wp), INTENT ( IN ) :: ta !< Local air temperature (degC) |
---|
1844 | REAL(wp), INTENT ( IN ) :: vp !< Local vapour pressure (hPa) |
---|
1845 | REAL(wp), INTENT ( IN ) :: tmrt !< Local mean radiant temperature (degC) |
---|
1846 | REAL(wp), INTENT ( IN ) :: ws !< Local wind velocitry (m/s) |
---|
1847 | REAL(wp), INTENT ( IN ) :: pair !< Local barometric air pressure (hPa) |
---|
1848 | ! |
---|
1849 | !-- Type of output of the argument list |
---|
1850 | REAL(wp), INTENT ( OUT ) :: perct_ij !< Perceived temperature (degC) |
---|
1851 | REAL(wp), INTENT ( OUT ) :: clo !< Clothing index (dimensionless) |
---|
1852 | ! |
---|
1853 | !-- Parameters for standard "Klima-Michel" |
---|
1854 | REAL(wp), PARAMETER :: eta = 0._wp !< Mechanical work efficiency for walking on flat ground (compare to Fanger (1972) pp 24f) |
---|
1855 | REAL(wp), PARAMETER :: actlev = 134.6862_wp !< Workload by activity per standardized surface (A_Du) |
---|
1856 | ! |
---|
1857 | !-- Type of program variables |
---|
1858 | REAL(wp), PARAMETER :: eps = 0.0005 !< Accuracy in clothing insulation (clo) for evaluation the root of Fanger's PMV (pmva=0) |
---|
1859 | REAL(wp) :: sclo !< summer clothing insulation |
---|
1860 | REAL(wp) :: wclo !< winter clothing insulation |
---|
1861 | REAL(wp) :: d_pmv !< PMV deviation (dimensionless --> PMV) |
---|
1862 | REAL(wp) :: svp_ta !< saturation vapor pressure (hPa) |
---|
1863 | REAL(wp) :: sult_lim !< threshold for sultrieness (hPa) |
---|
1864 | REAL(wp) :: dgtcm !< Mean deviation dependent on perct |
---|
1865 | REAL(wp) :: dgtcstd !< Mean deviation plus its standard deviation |
---|
1866 | REAL(wp) :: clon !< clo for neutral conditions (clo) |
---|
1867 | REAL(wp) :: ireq_minimal !< Minimal required clothing insulation (clo) |
---|
1868 | ! REAL(wp) :: clo_fanger !< clo for fanger subroutine, unused |
---|
1869 | REAL(wp) :: pmv_w !< Fangers predicted mean vote for winter clothing |
---|
1870 | REAL(wp) :: pmv_s !< Fangers predicted mean vote for summer clothing |
---|
1871 | REAL(wp) :: pmva !< adjusted predicted mean vote |
---|
1872 | REAL(wp) :: ptc !< perceived temp. for cold conditions (degree_C) |
---|
1873 | REAL(wp) :: d_std !< factor to threshold for sultriness |
---|
1874 | REAL(wp) :: pmvs !< pred. mean vote considering sultrieness |
---|
1875 | REAL(wp) :: top !< Gagge's operative temperatures (degree_C) |
---|
1876 | |
---|
1877 | INTEGER(iwp) :: ncount !< running index |
---|
1878 | INTEGER(iwp) :: nerr_cold !< error number (cold conditions) |
---|
1879 | INTEGER(iwp) :: nerr !< error number |
---|
1880 | |
---|
1881 | LOGICAL :: sultrieness |
---|
1882 | ! |
---|
1883 | !-- Initialise |
---|
1884 | perct_ij = bio_fill_value |
---|
1885 | |
---|
1886 | nerr = 0_iwp |
---|
1887 | ncount = 0_iwp |
---|
1888 | sultrieness = .FALSE. |
---|
1889 | ! |
---|
1890 | !-- Tresholds: clothing insulation (account for model inaccuracies) |
---|
1891 | ! |
---|
1892 | ! summer clothing |
---|
1893 | sclo = 0.44453_wp |
---|
1894 | ! |
---|
1895 | ! winter clothing |
---|
1896 | wclo = 1.76267_wp |
---|
1897 | ! |
---|
1898 | !-- decision: firstly calculate for winter or summer clothing |
---|
1899 | IF ( ta <= 10._wp ) THEN |
---|
1900 | ! |
---|
1901 | !-- First guess: winter clothing insulation: cold stress |
---|
1902 | clo = wclo |
---|
1903 | CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top ) |
---|
1904 | pmv_w = pmva |
---|
1905 | |
---|
1906 | IF ( pmva > 0._wp ) THEN |
---|
1907 | ! |
---|
1908 | !-- Case summer clothing insulation: heat load ? |
---|
1909 | clo = sclo |
---|
1910 | CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, & |
---|
1911 | top ) |
---|
1912 | pmv_s = pmva |
---|
1913 | IF ( pmva <= 0._wp ) THEN |
---|
1914 | ! |
---|
1915 | !-- Case: comfort achievable by varying clothing insulation |
---|
1916 | ! Between winter and summer set values |
---|
1917 | CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, & |
---|
1918 | pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo ) |
---|
1919 | IF ( ncount < 0_iwp ) THEN |
---|
1920 | nerr = -1_iwp |
---|
1921 | RETURN |
---|
1922 | ENDIF |
---|
1923 | ELSE IF ( pmva > 0.06_wp ) THEN |
---|
1924 | clo = 0.5_wp |
---|
1925 | CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, & |
---|
1926 | pmva, top ) |
---|
1927 | ENDIF |
---|
1928 | ELSE IF ( pmva < -0.11_wp ) THEN |
---|
1929 | clo = 1.75_wp |
---|
1930 | CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, & |
---|
1931 | top ) |
---|
1932 | ENDIF |
---|
1933 | ELSE |
---|
1934 | ! |
---|
1935 | !-- First guess: summer clothing insulation: heat load |
---|
1936 | clo = sclo |
---|
1937 | CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top ) |
---|
1938 | pmv_s = pmva |
---|
1939 | |
---|
1940 | IF ( pmva < 0._wp ) THEN |
---|
1941 | ! |
---|
1942 | !-- Case winter clothing insulation: cold stress ? |
---|
1943 | clo = wclo |
---|
1944 | CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, & |
---|
1945 | top ) |
---|
1946 | pmv_w = pmva |
---|
1947 | |
---|
1948 | IF ( pmva >= 0._wp ) THEN |
---|
1949 | ! |
---|
1950 | !-- Case: comfort achievable by varying clothing insulation |
---|
1951 | ! between winter and summer set values |
---|
1952 | CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, & |
---|
1953 | pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo ) |
---|
1954 | IF ( ncount < 0_iwp ) THEN |
---|
1955 | nerr = -1_iwp |
---|
1956 | RETURN |
---|
1957 | ENDIF |
---|
1958 | ELSE IF ( pmva < -0.11_wp ) THEN |
---|
1959 | clo = 1.75_wp |
---|
1960 | CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, & |
---|
1961 | pmva, top ) |
---|
1962 | ENDIF |
---|
1963 | ELSE IF ( pmva > 0.06_wp ) THEN |
---|
1964 | clo = 0.5_wp |
---|
1965 | CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, & |
---|
1966 | top ) |
---|
1967 | ENDIF |
---|
1968 | |
---|
1969 | ENDIF |
---|
1970 | ! |
---|
1971 | !-- Determine perceived temperature by regression equation + adjustments |
---|
1972 | pmvs = pmva |
---|
1973 | CALL perct_regression ( pmva, clo, perct_ij ) |
---|
1974 | ptc = perct_ij |
---|
1975 | IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN |
---|
1976 | ! |
---|
1977 | !-- Adjust for cold conditions according to Gagge 1986 |
---|
1978 | CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv ) |
---|
1979 | IF ( nerr_cold > 0_iwp ) nerr = -5_iwp |
---|
1980 | pmvs = pmva - d_pmv |
---|
1981 | IF ( pmvs > -0.11_wp ) THEN |
---|
1982 | d_pmv = 0._wp |
---|
1983 | pmvs = -0.11_wp |
---|
1984 | ENDIF |
---|
1985 | CALL perct_regression ( pmvs, clo, perct_ij ) |
---|
1986 | ENDIF |
---|
1987 | ! clo_fanger = clo |
---|
1988 | clon = clo |
---|
1989 | IF ( clo > 0.5_wp .AND. perct_ij <= 8.73_wp ) THEN |
---|
1990 | ! |
---|
1991 | !-- Required clothing insulation (ireq) is exclusively defined for |
---|
1992 | ! operative temperatures (top) less 10 (C) for a |
---|
1993 | ! reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s |
---|
1994 | clon = ireq_neutral ( perct_ij, ireq_minimal, nerr ) |
---|
1995 | clo = clon |
---|
1996 | ENDIF |
---|
1997 | CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim ) |
---|
1998 | sultrieness = .FALSE. |
---|
1999 | d_std = -99._wp |
---|
2000 | IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN |
---|
2001 | ! |
---|
2002 | !-- Adjust for warm/humid conditions according to Gagge 1986 |
---|
2003 | CALL saturation_vapor_pressure ( ta, svp_ta ) |
---|
2004 | d_pmv = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr ) |
---|
2005 | pmvs = pmva + d_pmv |
---|
2006 | CALL perct_regression ( pmvs, clo, perct_ij ) |
---|
2007 | IF ( sult_lim < 99._wp ) THEN |
---|
2008 | IF ( (perct_ij - ptc) > sult_lim ) sultrieness = .TRUE. |
---|
2009 | ! |
---|
2010 | !-- Set factor to threshold for sultriness |
---|
2011 | IF ( dgtcstd /= 0_iwp ) THEN |
---|
2012 | d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd |
---|
2013 | ENDIF |
---|
2014 | ENDIF |
---|
2015 | ENDIF |
---|
2016 | |
---|
2017 | END SUBROUTINE calculate_perct_static |
---|
2018 | |
---|
2019 | !------------------------------------------------------------------------------! |
---|
2020 | ! Description: |
---|
2021 | ! ------------ |
---|
2022 | !> The SUBROUTINE calculates the saturation water vapour pressure |
---|
2023 | !> (hPa = hecto Pascal) for a given temperature ta (degC). |
---|
2024 | !> For example, ta can be the air temperature or the dew point temperature. |
---|
2025 | !------------------------------------------------------------------------------! |
---|
2026 | SUBROUTINE saturation_vapor_pressure( ta, svp_ta ) |
---|
2027 | |
---|
2028 | IMPLICIT NONE |
---|
2029 | |
---|
2030 | REAL(wp), INTENT ( IN ) :: ta !< ambient air temperature (degC) |
---|
2031 | REAL(wp), INTENT ( OUT ) :: svp_ta !< saturation water vapour pressure (hPa) |
---|
2032 | |
---|
2033 | REAL(wp) :: b |
---|
2034 | REAL(wp) :: c |
---|
2035 | |
---|
2036 | |
---|
2037 | IF ( ta < 0._wp ) THEN |
---|
2038 | ! |
---|
2039 | !-- ta < 0 (degC): saturation water vapour pressure over ice |
---|
2040 | b = 17.84362_wp |
---|
2041 | c = 245.425_wp |
---|
2042 | ELSE |
---|
2043 | ! |
---|
2044 | !-- ta >= 0 (degC): saturation water vapour pressure over water |
---|
2045 | b = 17.08085_wp |
---|
2046 | c = 234.175_wp |
---|
2047 | ENDIF |
---|
2048 | ! |
---|
2049 | !-- Saturation water vapour pressure |
---|
2050 | svp_ta = 6.1078_wp * EXP ( b * ta / ( c + ta ) ) |
---|
2051 | |
---|
2052 | END SUBROUTINE saturation_vapor_pressure |
---|
2053 | |
---|
2054 | !------------------------------------------------------------------------------! |
---|
2055 | ! Description: |
---|
2056 | ! ------------ |
---|
2057 | !> Find the clothing insulation value clo_res (clo) to make Fanger's Predicted |
---|
2058 | !> Mean Vote (PMV) equal comfort (pmva=0) for actual meteorological conditions |
---|
2059 | !> (ta,tmrt, vp, ws, pair) and values of individual's activity level |
---|
2060 | !------------------------------------------------------------------------------! |
---|
2061 | SUBROUTINE iso_ridder( ta, tmrt, vp, ws, pair, actlev, eta, sclo, & |
---|
2062 | pmv_s, wclo, pmv_w, eps, pmva, top, nerr, & |
---|
2063 | clo_res ) |
---|
2064 | |
---|
2065 | IMPLICIT NONE |
---|
2066 | ! |
---|
2067 | !-- Input variables of argument list: |
---|
2068 | REAL(wp), INTENT ( IN ) :: ta !< Ambient temperature (degC) |
---|
2069 | REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (degC) |
---|
2070 | REAL(wp), INTENT ( IN ) :: vp !< Water vapour pressure (hPa) |
---|
2071 | REAL(wp), INTENT ( IN ) :: ws !< Wind speed (m/s) 1 m above ground |
---|
2072 | REAL(wp), INTENT ( IN ) :: pair !< Barometric pressure (hPa) |
---|
2073 | REAL(wp), INTENT ( IN ) :: actlev !< Individuals activity level per unit surface area (W/m2) |
---|
2074 | REAL(wp), INTENT ( IN ) :: eta !< Individuals work efficiency (dimensionless) |
---|
2075 | REAL(wp), INTENT ( IN ) :: sclo !< Lower threshold of bracketing clothing insulation (clo) |
---|
2076 | REAL(wp), INTENT ( IN ) :: wclo !< Upper threshold of bracketing clothing insulation (clo) |
---|
2077 | REAL(wp), INTENT ( IN ) :: eps !< (0.05) accuracy in clothing insulation (clo) for |
---|
2078 | ! evaluation the root of Fanger's PMV (pmva=0) |
---|
2079 | REAL(wp), INTENT ( IN ) :: pmv_w !< Fanger's PMV corresponding to wclo |
---|
2080 | REAL(wp), INTENT ( IN ) :: pmv_s !< Fanger's PMV corresponding to sclo |
---|
2081 | ! |
---|
2082 | ! Output variables of argument list: |
---|
2083 | REAL(wp), INTENT ( OUT ) :: pmva !< 0 (set to zero, because clo is evaluated for comfort) |
---|
2084 | REAL(wp), INTENT ( OUT ) :: top !< Operative temperature (degC) at found root of Fanger's PMV |
---|
2085 | REAL(wp), INTENT ( OUT ) :: clo_res !< Resulting clothing insulation value (clo) |
---|
2086 | INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag |
---|
2087 | ! nerr >= 0, o.k., and nerr is the number of iterations for |
---|
2088 | ! convergence |
---|
2089 | ! nerr = -1: error = malfunction of Ridder's convergence method |
---|
2090 | ! nerr = -2: error = maximum iterations (max_iteration) exceeded |
---|
2091 | ! nerr = -3: error = root not bracketed between sclo and wclo |
---|
2092 | ! |
---|
2093 | !-- Type of program variables |
---|
2094 | INTEGER(iwp), PARAMETER :: max_iteration = 15_iwp !< max number of iterations |
---|
2095 | REAL(wp), PARAMETER :: guess_0 = -1.11e30_wp !< initial guess |
---|
2096 | REAL(wp) :: x_ridder !< current guess for clothing insulation (clo) |
---|
2097 | REAL(wp) :: clo_lower !< lower limit of clothing insulation (clo) |
---|
2098 | REAL(wp) :: clo_upper !< upper limit of clothing insulation (clo) |
---|
2099 | REAL(wp) :: x_lower !< lower guess for clothing insulation (clo) |
---|
2100 | REAL(wp) :: x_upper !< upper guess for clothing insulation (clo) |
---|
2101 | REAL(wp) :: x_average !< average of x_lower and x_upper (clo) |
---|
2102 | REAL(wp) :: x_new !< preliminary result for clothing insulation (clo) |
---|
2103 | REAL(wp) :: y_lower !< predicted mean vote for summer clothing |
---|
2104 | REAL(wp) :: y_upper !< predicted mean vote for winter clothing |
---|
2105 | REAL(wp) :: y_average !< average of y_lower and y_upper |
---|
2106 | REAL(wp) :: y_new !< preliminary result for pred. mean vote |
---|
2107 | REAL(wp) :: sroot !< sqrt of PMV-guess |
---|
2108 | INTEGER(iwp) :: j !< running index |
---|
2109 | ! |
---|
2110 | !-- Initialise |
---|
2111 | nerr = 0_iwp |
---|
2112 | ! |
---|
2113 | !-- Set pmva = 0 (comfort): Root of PMV depending on clothing insulation |
---|
2114 | x_ridder = bio_fill_value |
---|
2115 | pmva = 0._wp |
---|
2116 | clo_lower = sclo |
---|
2117 | y_lower = pmv_s |
---|
2118 | clo_upper = wclo |
---|
2119 | y_upper = pmv_w |
---|
2120 | IF ( ( y_lower > 0._wp .AND. y_upper < 0._wp ) .OR. & |
---|
2121 | ( y_lower < 0._wp .AND. y_upper > 0._wp ) ) THEN |
---|
2122 | x_lower = clo_lower |
---|
2123 | x_upper = clo_upper |
---|
2124 | x_ridder = guess_0 |
---|
2125 | |
---|
2126 | DO j = 1_iwp, max_iteration |
---|
2127 | x_average = 0.5_wp * ( x_lower + x_upper ) |
---|
2128 | CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta, & |
---|
2129 | y_average, top ) |
---|
2130 | sroot = SQRT ( y_average**2 - y_lower * y_upper ) |
---|
2131 | IF ( sroot == 0._wp ) THEN |
---|
2132 | clo_res = x_average |
---|
2133 | nerr = j |
---|
2134 | RETURN |
---|
2135 | ENDIF |
---|
2136 | x_new = x_average + ( x_average - x_lower ) * & |
---|
2137 | ( SIGN ( 1._wp, y_lower - y_upper ) * y_average / sroot ) |
---|
2138 | IF ( ABS ( x_new - x_ridder ) <= eps ) THEN |
---|
2139 | clo_res = x_ridder |
---|
2140 | nerr = j |
---|
2141 | RETURN |
---|
2142 | ENDIF |
---|
2143 | x_ridder = x_new |
---|
2144 | CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta, & |
---|
2145 | y_new, top ) |
---|
2146 | IF ( y_new == 0._wp ) THEN |
---|
2147 | clo_res = x_ridder |
---|
2148 | nerr = j |
---|
2149 | RETURN |
---|
2150 | ENDIF |
---|
2151 | IF ( SIGN ( y_average, y_new ) /= y_average ) THEN |
---|
2152 | x_lower = x_average |
---|
2153 | y_lower = y_average |
---|
2154 | x_upper = x_ridder |
---|
2155 | y_upper = y_new |
---|
2156 | ELSE IF ( SIGN ( y_lower, y_new ) /= y_lower ) THEN |
---|
2157 | x_upper = x_ridder |
---|
2158 | y_upper = y_new |
---|
2159 | ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper ) THEN |
---|
2160 | x_lower = x_ridder |
---|
2161 | y_lower = y_new |
---|
2162 | ELSE |
---|
2163 | ! |
---|
2164 | !-- Never get here in x_ridder: singularity in y |
---|
2165 | nerr = -1_iwp |
---|
2166 | clo_res = x_ridder |
---|
2167 | RETURN |
---|
2168 | ENDIF |
---|
2169 | IF ( ABS ( x_upper - x_lower ) <= eps ) THEN |
---|
2170 | clo_res = x_ridder |
---|
2171 | nerr = j |
---|
2172 | RETURN |
---|
2173 | ENDIF |
---|
2174 | ENDDO |
---|
2175 | ! |
---|
2176 | !-- x_ridder exceed maximum iterations |
---|
2177 | nerr = -2_iwp |
---|
2178 | clo_res = y_new |
---|
2179 | RETURN |
---|
2180 | ELSE IF ( y_lower == 0. ) THEN |
---|
2181 | x_ridder = clo_lower |
---|
2182 | ELSE IF ( y_upper == 0. ) THEN |
---|
2183 | x_ridder = clo_upper |
---|
2184 | ELSE |
---|
2185 | ! |
---|
2186 | !-- x_ridder not bracketed by u_clo and o_clo |
---|
2187 | nerr = -3_iwp |
---|
2188 | clo_res = x_ridder |
---|
2189 | RETURN |
---|
2190 | ENDIF |
---|
2191 | |
---|
2192 | END SUBROUTINE iso_ridder |
---|
2193 | |
---|
2194 | !------------------------------------------------------------------------------! |
---|
2195 | ! Description: |
---|
2196 | ! ------------ |
---|
2197 | !> Regression relations between perceived temperature (perct) and (adjusted) |
---|
2198 | !> PMV. The regression presumes the Klima-Michel settings for reference |
---|
2199 | !> individual and reference environment. |
---|
2200 | !------------------------------------------------------------------------------! |
---|
2201 | SUBROUTINE perct_regression( pmv, clo, perct_ij ) |
---|
2202 | |
---|
2203 | IMPLICIT NONE |
---|
2204 | |
---|
2205 | REAL(wp), INTENT ( IN ) :: pmv !< Fangers predicted mean vote (dimensionless) |
---|
2206 | REAL(wp), INTENT ( IN ) :: clo !< clothing insulation index (clo) |
---|
2207 | |
---|
2208 | REAL(wp), INTENT ( OUT ) :: perct_ij !< perct (degC) corresponding to given PMV / clo |
---|
2209 | |
---|
2210 | IF ( pmv <= -0.11_wp ) THEN |
---|
2211 | perct_ij = 5.805_wp + 12.6784_wp * pmv |
---|
2212 | ELSE |
---|
2213 | IF ( pmv >= + 0.01_wp ) THEN |
---|
2214 | perct_ij = 16.826_wp + 6.163_wp * pmv |
---|
2215 | ELSE |
---|
2216 | perct_ij = 21.258_wp - 9.558_wp * clo |
---|
2217 | ENDIF |
---|
2218 | ENDIF |
---|
2219 | |
---|
2220 | END SUBROUTINE perct_regression |
---|
2221 | |
---|
2222 | !------------------------------------------------------------------------------! |
---|
2223 | ! Description: |
---|
2224 | ! ------------ |
---|
2225 | !> FANGER.F90 |
---|
2226 | !> |
---|
2227 | !> SI-VERSION: ACTLEV W m-2, DAMPFDRUCK hPa |
---|
2228 | !> Berechnet das aktuelle Predicted Mean Vote nach Fanger |
---|
2229 | !> |
---|
2230 | !> The case of free convection (ws < 0.1 m/s) is dealt with ws = 0.1 m/s |
---|
2231 | !------------------------------------------------------------------------------! |
---|
2232 | SUBROUTINE fanger( ta, tmrt, pa, in_ws, pair, in_clo, actlev, eta, pmva, top ) |
---|
2233 | |
---|
2234 | IMPLICIT NONE |
---|
2235 | ! |
---|
2236 | !-- Input variables of argument list: |
---|
2237 | REAL(wp), INTENT ( IN ) :: ta !< Ambient air temperature (degC) |
---|
2238 | REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (degC) |
---|
2239 | REAL(wp), INTENT ( IN ) :: pa !< Water vapour pressure (hPa) |
---|
2240 | REAL(wp), INTENT ( IN ) :: pair !< Barometric pressure (hPa) at site |
---|
2241 | REAL(wp), INTENT ( IN ) :: in_ws !< Wind speed (m/s) 1 m above ground |
---|
2242 | REAL(wp), INTENT ( IN ) :: in_clo !< Clothing insulation (clo) |
---|
2243 | REAL(wp), INTENT ( IN ) :: actlev !< Individuals activity level per unit surface area (W/m2) |
---|
2244 | REAL(wp), INTENT ( IN ) :: eta !< Individuals mechanical work efficiency (dimensionless) |
---|
2245 | ! |
---|
2246 | !-- Output variables of argument list: |
---|
2247 | REAL(wp), INTENT ( OUT ) :: pmva !< Actual Predicted Mean Vote (PMV, |
---|
2248 | ! dimensionless) according to Fanger corresponding to meteorological |
---|
2249 | ! (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta) |
---|
2250 | REAL(wp), INTENT ( OUT ) :: top !< operative temperature (degC) |
---|
2251 | ! |
---|
2252 | !-- Internal variables |
---|
2253 | REAL(wp) :: f_cl !< Increase in surface due to clothing (factor) |
---|
2254 | REAL(wp) :: heat_convection !< energy loss by autocnvection (W) |
---|
2255 | REAL(wp) :: activity !< persons activity (must stay == actlev, W) |
---|
2256 | REAL(wp) :: t_skin_aver !< average skin temperature (degree_C) |
---|
2257 | REAL(wp) :: bc !< preliminary result storage |
---|
2258 | REAL(wp) :: cc !< preliminary result storage |
---|
2259 | REAL(wp) :: dc !< preliminary result storage |
---|
2260 | REAL(wp) :: ec !< preliminary result storage |
---|
2261 | REAL(wp) :: gc !< preliminary result storage |
---|
2262 | REAL(wp) :: t_clothing !< clothing temperature (degree_C) |
---|
2263 | REAL(wp) :: hr !< radiational heat resistence |
---|
2264 | REAL(wp) :: clo !< clothing insulation index (clo) |
---|
2265 | REAL(wp) :: ws !< wind speed (m/s) |
---|
2266 | REAL(wp) :: z1 !< Empiric factor for the adaption of the heat |
---|
2267 | ! ballance equation to the psycho-physical scale (Equ. 40 in FANGER) |
---|
2268 | REAL(wp) :: z2 !< Water vapour diffution through the skin |
---|
2269 | REAL(wp) :: z3 !< Sweat evaporation from the skin surface |
---|
2270 | REAL(wp) :: z4 !< Loss of latent heat through respiration |
---|
2271 | REAL(wp) :: z5 !< Loss of radiational heat |
---|
2272 | REAL(wp) :: z6 !< Heat loss through forced convection |
---|
2273 | INTEGER(iwp) :: i !< running index |
---|
2274 | ! |
---|
2275 | !-- Clo must be > 0. to avoid div. by 0! |
---|
2276 | clo = in_clo |
---|
2277 | IF ( clo <= 0._wp ) clo = .001_wp |
---|
2278 | ! |
---|
2279 | !-- f_cl = Increase in surface due to clothing |
---|
2280 | f_cl = 1._wp + .15_wp * clo |
---|
2281 | ! |
---|
2282 | !-- Case of free convection (ws < 0.1 m/s ) not considered |
---|
2283 | ws = in_ws |
---|
2284 | IF ( ws < .1_wp ) THEN |
---|
2285 | ws = .1_wp |
---|
2286 | ENDIF |
---|
2287 | ! |
---|
2288 | !-- Heat_convection = forced convection |
---|
2289 | heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp ) |
---|
2290 | ! |
---|
2291 | !-- Activity = inner heat produktion per standardized surface |
---|
2292 | activity = actlev * ( 1._wp - eta ) |
---|
2293 | ! |
---|
2294 | !-- T_skin_aver = average skin temperature |
---|
2295 | t_skin_aver = 35.7_wp - .0275_wp * activity |
---|
2296 | ! |
---|
2297 | !-- Calculation of constants for evaluation below |
---|
2298 | bc = .155_wp * clo * 3.96_wp * 10._wp**( -8 ) * f_cl |
---|
2299 | cc = f_cl * heat_convection |
---|
2300 | ec = .155_wp * clo |
---|
2301 | dc = ( 1._wp + ec * cc ) / bc |
---|
2302 | gc = ( t_skin_aver + bc * ( tmrt + degc_to_k )**4 + ec * cc * ta ) / bc |
---|
2303 | ! |
---|
2304 | !-- Calculation of clothing surface temperature (t_clothing) based on |
---|
2305 | ! Newton-approximation with air temperature as initial guess |
---|
2306 | t_clothing = ta |
---|
2307 | DO i = 1, 3 |
---|
2308 | t_clothing = t_clothing - ( ( t_clothing + degc_to_k )**4 + t_clothing & |
---|
2309 | * dc - gc ) / ( 4._wp * ( t_clothing + degc_to_k )**3 + dc ) |
---|
2310 | ENDDO |
---|
2311 | ! |
---|
2312 | !-- Empiric factor for the adaption of the heat ballance equation |
---|
2313 | ! to the psycho-physical scale (Equ. 40 in FANGER) |
---|
2314 | z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp ) |
---|
2315 | ! |
---|
2316 | !-- Water vapour diffution through the skin |
---|
2317 | z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa ) |
---|
2318 | ! |
---|
2319 | !-- Sweat evaporation from the skin surface |
---|
2320 | z3 = .42_wp * ( activity - 58._wp ) |
---|
2321 | ! |
---|
2322 | !-- Loss of latent heat through respiration |
---|
2323 | z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev * & |
---|
2324 | ( 34._wp - ta ) |
---|
2325 | ! |
---|
2326 | !-- Loss of radiational heat |
---|
2327 | z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt + & |
---|
2328 | degc_to_k )**4 ) |
---|
2329 | IF ( ABS ( t_clothing - tmrt ) > 0._wp ) THEN |
---|
2330 | hr = z5 / f_cl / ( t_clothing - tmrt ) |
---|
2331 | ELSE |
---|
2332 | hr = 0._wp |
---|
2333 | ENDIF |
---|
2334 | ! |
---|
2335 | !-- Heat loss through forced convection cc*(t_clothing-TT) |
---|
2336 | z6 = cc * ( t_clothing - ta ) |
---|
2337 | ! |
---|
2338 | !-- Predicted Mean Vote |
---|
2339 | pmva = z1 * ( activity - z2 - z3 - z4 - z5 - z6 ) |
---|
2340 | ! |
---|
2341 | !-- Operative temperatur |
---|
2342 | top = ( hr * tmrt + heat_convection * ta ) / ( hr + heat_convection ) |
---|
2343 | |
---|
2344 | END SUBROUTINE fanger |
---|
2345 | |
---|
2346 | !------------------------------------------------------------------------------! |
---|
2347 | ! Description: |
---|
2348 | ! ------------ |
---|
2349 | !> For pmva > 0 and clo =0.5 the increment (deltapmv) is calculated |
---|
2350 | !> that converts pmva into Gagge's et al. (1986) PMV*. |
---|
2351 | !------------------------------------------------------------------------------! |
---|
2352 | REAL(wp) FUNCTION deltapmv( pmva, ta, vp, svp_ta, tmrt, ws, nerr ) |
---|
2353 | |
---|
2354 | IMPLICIT NONE |
---|
2355 | ! |
---|
2356 | !-- Input variables of argument list: |
---|
2357 | REAL(wp), INTENT ( IN ) :: pmva !< Actual Predicted Mean Vote (PMV) according to Fanger |
---|
2358 | REAL(wp), INTENT ( IN ) :: ta !< Ambient temperature (degC) at screen level |
---|
2359 | REAL(wp), INTENT ( IN ) :: vp !< Water vapour pressure (hPa) at screen level |
---|
2360 | REAL(wp), INTENT ( IN ) :: svp_ta !< Saturation water vapour pressure (hPa) at ta |
---|
2361 | REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (degC) at screen level |
---|
2362 | REAL(wp), INTENT ( IN ) :: ws !< Wind speed (m/s) 1 m above ground |
---|
2363 | ! |
---|
2364 | !-- Output variables of argument list: |
---|
2365 | INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag |
---|
2366 | ! 0 = o.k. |
---|
2367 | ! -2 = pmva outside valid regression range |
---|
2368 | ! -3 = rel. humidity set to 5 % or 95 %, respectively |
---|
2369 | ! -4 = deltapmv set to avoid pmvs < 0 |
---|
2370 | ! |
---|
2371 | !-- Internal variable types: |
---|
2372 | REAL(wp) :: pmv !< |
---|
2373 | REAL(wp) :: pa_p50 !< |
---|
2374 | REAL(wp) :: pa !< |
---|
2375 | REAL(wp) :: apa !< |
---|
2376 | REAL(wp) :: dapa !< |
---|
2377 | REAL(wp) :: sqvel !< |
---|
2378 | REAL(wp) :: dtmrt !< |
---|
2379 | REAL(wp) :: p10 !< |
---|
2380 | REAL(wp) :: p95 !< |
---|
2381 | REAL(wp) :: gew !< |
---|
2382 | REAL(wp) :: gew2 !< |
---|
2383 | REAL(wp) :: dpmv_1 !< |
---|
2384 | REAL(wp) :: dpmv_2 !< |
---|
2385 | REAL(wp) :: pmvs !< |
---|
2386 | REAL(wp) :: bpmv(0:7) !< |
---|
2387 | REAL(wp) :: bpa_p50(0:7) !< |
---|
2388 | REAL(wp) :: bpa(0:7) !< |
---|
2389 | REAL(wp) :: bapa(0:7) !< |
---|
2390 | REAL(wp) :: bdapa(0:7) !< |
---|
2391 | REAL(wp) :: bsqvel(0:7) !< |
---|
2392 | REAL(wp) :: bta(0:7) !< |
---|
2393 | REAL(wp) :: bdtmrt(0:7) !< |
---|
2394 | REAL(wp) :: aconst(0:7) !< |
---|
2395 | INTEGER(iwp) :: nreg !< |
---|
2396 | |
---|
2397 | DATA bpmv / & |
---|
2398 | -0.0556602_wp, -0.1528680_wp, -0.2336104_wp, -0.2789387_wp, -0.3551048_wp,& |
---|
2399 | -0.4304076_wp, -0.4884961_wp, -0.4897495_wp / |
---|
2400 | DATA bpa_p50 / & |
---|
2401 | -0.1607154_wp, -0.4177296_wp, -0.4120541_wp, -0.0886564_wp, +0.4285938_wp,& |
---|
2402 | +0.6281256_wp, +0.5067361_wp, +0.3965169_wp / |
---|
2403 | DATA bpa / & |
---|
2404 | +0.0580284_wp, +0.0836264_wp, +0.1009919_wp, +0.1020777_wp, +0.0898681_wp,& |
---|
2405 | +0.0839116_wp, +0.0853258_wp, +0.0866589_wp / |
---|
2406 | DATA bapa / & |
---|
2407 | -1.7838788_wp, -2.9306231_wp, -1.6350334_wp, +0.6211547_wp, +3.3918083_wp,& |
---|
2408 | +5.5521025_wp, +8.4897418_wp, +16.6265851_wp / |
---|
2409 | DATA bdapa / & |
---|
2410 | +1.6752720_wp, +2.7379504_wp, +1.2940526_wp, -1.0985759_wp, -3.9054732_wp,& |
---|
2411 | -6.0403012_wp, -8.9437119_wp, -17.0671201_wp / |
---|
2412 | DATA bsqvel / & |
---|
2413 | -0.0315598_wp, -0.0286272_wp, -0.0009228_wp, +0.0483344_wp, +0.0992366_wp,& |
---|
2414 | +0.1491379_wp, +0.1951452_wp, +0.2133949_wp / |
---|
2415 | DATA bta / & |
---|
2416 | +0.0953986_wp, +0.1524760_wp, +0.0564241_wp, -0.0893253_wp, -0.2398868_wp,& |
---|
2417 | -0.3515237_wp, -0.5095144_wp, -0.9469258_wp / |
---|
2418 | DATA bdtmrt / & |
---|
2419 | -0.0004672_wp, -0.0000514_wp, -0.0018037_wp, -0.0049440_wp, -0.0069036_wp,& |
---|
2420 | -0.0075844_wp, -0.0079602_wp, -0.0089439_wp / |
---|
2421 | DATA aconst / & |
---|
2422 | +1.8686215_wp, +3.4260713_wp, +2.0116185_wp, -0.7777552_wp, -4.6715853_wp,& |
---|
2423 | -7.7314281_wp, -11.7602578_wp, -23.5934198_wp / |
---|
2424 | ! |
---|
2425 | !-- Test for compliance with regression range |
---|
2426 | IF ( pmva < -1.0_wp .OR. pmva > 7.0_wp ) THEN |
---|
2427 | nerr = -2_iwp |
---|
2428 | ELSE |
---|
2429 | nerr = 0_iwp |
---|
2430 | ENDIF |
---|
2431 | ! |
---|
2432 | !-- Initialise classic PMV |
---|
2433 | pmv = pmva |
---|
2434 | ! |
---|
2435 | !-- Water vapour pressure of air |
---|
2436 | p10 = 0.05_wp * svp_ta |
---|
2437 | p95 = 1.00_wp * svp_ta |
---|
2438 | IF ( vp >= p10 .AND. vp <= p95 ) THEN |
---|
2439 | pa = vp |
---|
2440 | ELSE |
---|
2441 | nerr = -3_iwp |
---|
2442 | IF ( vp < p10 ) THEN |
---|
2443 | ! |
---|
2444 | !-- Due to conditions of regression: r.H. >= 5 % |
---|
2445 | pa = p10 |
---|
2446 | ELSE |
---|
2447 | ! |
---|
2448 | !-- Due to conditions of regression: r.H. <= 95 % |
---|
2449 | pa = p95 |
---|
2450 | ENDIF |
---|
2451 | ENDIF |
---|
2452 | IF ( pa > 0._wp ) THEN |
---|
2453 | ! |
---|
2454 | !-- Natural logarithm of pa |
---|
2455 | apa = LOG ( pa ) |
---|
2456 | ELSE |
---|
2457 | apa = -5._wp |
---|
2458 | ENDIF |
---|
2459 | ! |
---|
2460 | !-- Ratio actual water vapour pressure to that of a r.H. of 50 % |
---|
2461 | pa_p50 = 0.5_wp * svp_ta |
---|
2462 | IF ( pa_p50 > 0._wp .AND. pa > 0._wp ) THEN |
---|
2463 | dapa = apa - LOG ( pa_p50 ) |
---|
2464 | pa_p50 = pa / pa_p50 |
---|
2465 | ELSE |
---|
2466 | dapa = -5._wp |
---|
2467 | pa_p50 = 0._wp |
---|
2468 | ENDIF |
---|
2469 | ! |
---|
2470 | !-- Square root of wind velocity |
---|
2471 | IF ( ws >= 0._wp ) THEN |
---|
2472 | sqvel = SQRT ( ws ) |
---|
2473 | ELSE |
---|
2474 | sqvel = 0._wp |
---|
2475 | ENDIF |
---|
2476 | ! |
---|
2477 | !-- Air temperature |
---|
2478 | ! ta = ta |
---|
2479 | !-- Difference mean radiation to air temperature |
---|
2480 | dtmrt = tmrt - ta |
---|
2481 | ! |
---|
2482 | !-- Select the valid regression coefficients |
---|
2483 | nreg = INT ( pmv ) |
---|
2484 | IF ( nreg < 0_iwp ) THEN |
---|
2485 | ! |
---|
2486 | !-- value of the FUNCTION in the case pmv <= -1 |
---|
2487 | deltapmv = 0._wp |
---|
2488 | RETURN |
---|
2489 | ENDIF |
---|
2490 | gew = MOD ( pmv, 1._wp ) |
---|
2491 | IF ( gew < 0._wp ) gew = 0._wp |
---|
2492 | IF ( nreg > 5_iwp ) THEN |
---|
2493 | ! nreg=6 |
---|
2494 | nreg = 5_iwp |
---|
2495 | gew = pmv - 5._wp |
---|
2496 | gew2 = pmv - 6._wp |
---|
2497 | IF ( gew2 > 0_iwp ) THEN |
---|
2498 | gew = ( gew - gew2 ) / gew |
---|
2499 | ENDIF |
---|
2500 | ENDIF |
---|
2501 | ! |
---|
2502 | !-- Regression valid for 0. <= pmv <= 6. |
---|
2503 | dpmv_1 = & |
---|
2504 | + bpa ( nreg ) * pa & |
---|
2505 | + bpmv ( nreg ) * pmv & |
---|
2506 | + bapa ( nreg ) * apa & |
---|
2507 | + bta ( nreg ) * ta & |
---|
2508 | + bdtmrt ( nreg ) * dtmrt & |
---|
2509 | + bdapa ( nreg ) * dapa & |
---|
2510 | + bsqvel ( nreg ) * sqvel & |
---|
2511 | + bpa_p50 ( nreg ) * pa_p50 & |
---|
2512 | + aconst ( nreg ) |
---|
2513 | |
---|
2514 | dpmv_2 = 0._wp |
---|
2515 | IF ( nreg < 6_iwp ) THEN |
---|
2516 | dpmv_2 = & |
---|
2517 | + bpa ( nreg + 1_iwp ) * pa & |
---|
2518 | + bpmv ( nreg + 1_iwp ) * pmv & |
---|
2519 | + bapa ( nreg + 1_iwp ) * apa & |
---|
2520 | + bta ( nreg + 1_iwp ) * ta & |
---|
2521 | + bdtmrt ( nreg + 1_iwp ) * dtmrt & |
---|
2522 | + bdapa ( nreg + 1_iwp ) * dapa & |
---|
2523 | + bsqvel ( nreg + 1_iwp ) * sqvel & |
---|
2524 | + bpa_p50 ( nreg + 1_iwp ) * pa_p50 & |
---|
2525 | + aconst ( nreg + 1_iwp ) |
---|
2526 | ENDIF |
---|
2527 | ! |
---|
2528 | !-- Calculate pmv modification |
---|
2529 | deltapmv = ( 1._wp - gew ) * dpmv_1 + gew * dpmv_2 |
---|
2530 | pmvs = pmva + deltapmv |
---|
2531 | IF ( ( pmvs ) < 0._wp ) THEN |
---|
2532 | ! |
---|
2533 | !-- Prevent negative pmv* due to problems with clothing insulation |
---|
2534 | nerr = -4_iwp |
---|
2535 | IF ( pmvs > -0.11_wp ) THEN |
---|
2536 | ! |
---|
2537 | !-- Threshold from FUNCTION perct_regression for winter clothing insulation |
---|
2538 | deltapmv = deltapmv + 0.11_wp |
---|
2539 | ELSE |
---|
2540 | ! |
---|
2541 | !-- Set pmvs to "0" for compliance with summer clothing insulation |
---|
2542 | deltapmv = -1._wp * pmva |
---|
2543 | ENDIF |
---|
2544 | ENDIF |
---|
2545 | |
---|
2546 | END FUNCTION deltapmv |
---|
2547 | |
---|
2548 | !------------------------------------------------------------------------------! |
---|
2549 | ! Description: |
---|
2550 | ! ------------ |
---|
2551 | !> The subroutine "calc_sultr" returns a threshold value to perceived |
---|
2552 | !> temperature allowing to decide whether the actual perceived temperature |
---|
2553 | !> is linked to perecption of sultriness. The threshold values depends |
---|
2554 | !> on the Fanger's classical PMV, expressed here as perceived temperature |
---|
2555 | !> perct. |
---|
2556 | !------------------------------------------------------------------------------! |
---|
2557 | SUBROUTINE calc_sultr( perct_ij, dperctm, dperctstd, sultr_res ) |
---|
2558 | |
---|
2559 | IMPLICIT NONE |
---|
2560 | ! |
---|
2561 | !-- Input of the argument list: |
---|
2562 | REAL(wp), INTENT ( IN ) :: perct_ij !< Classical perceived temperature: Base is Fanger's PMV |
---|
2563 | ! |
---|
2564 | !-- Additional output variables of argument list: |
---|
2565 | REAL(wp), INTENT ( OUT ) :: dperctm !< Mean deviation perct (classical gt) to gt* (rational gt |
---|
2566 | ! calculated based on Gagge's rational PMV*) |
---|
2567 | REAL(wp), INTENT ( OUT ) :: dperctstd !< dperctm plus its standard deviation times a factor |
---|
2568 | ! determining the significance to perceive sultriness |
---|
2569 | REAL(wp), INTENT ( OUT ) :: sultr_res |
---|
2570 | ! |
---|
2571 | !-- Types of coefficients mean deviation: third order polynomial |
---|
2572 | REAL(wp), PARAMETER :: dperctka = +7.5776086_wp |
---|
2573 | REAL(wp), PARAMETER :: dperctkb = -0.740603_wp |
---|
2574 | REAL(wp), PARAMETER :: dperctkc = +0.0213324_wp |
---|
2575 | REAL(wp), PARAMETER :: dperctkd = -0.00027797237_wp |
---|
2576 | ! |
---|
2577 | !-- Types of coefficients mean deviation plus standard deviation |
---|
2578 | ! regression coefficients: third order polynomial |
---|
2579 | REAL(wp), PARAMETER :: dperctsa = +0.0268918_wp |
---|
2580 | REAL(wp), PARAMETER :: dperctsb = +0.0465957_wp |
---|
2581 | REAL(wp), PARAMETER :: dperctsc = -0.00054709752_wp |
---|
2582 | REAL(wp), PARAMETER :: dperctsd = +0.0000063714823_wp |
---|
2583 | ! |
---|
2584 | !-- Factor to mean standard deviation defining SIGNificance for |
---|
2585 | ! sultriness |
---|
2586 | REAL(wp), PARAMETER :: faktor = 1._wp |
---|
2587 | ! |
---|
2588 | !-- Initialise |
---|
2589 | sultr_res = +99._wp |
---|
2590 | dperctm = 0._wp |
---|
2591 | dperctstd = 999999._wp |
---|
2592 | |
---|
2593 | IF ( perct_ij < 16.826_wp .OR. perct_ij > 56._wp ) THEN |
---|
2594 | ! |
---|
2595 | !-- Unallowed value of classical perct! |
---|
2596 | RETURN |
---|
2597 | ENDIF |
---|
2598 | ! |
---|
2599 | !-- Mean deviation dependent on perct |
---|
2600 | dperctm = dperctka + dperctkb * perct_ij + dperctkc * perct_ij**2._wp + & |
---|
2601 | dperctkd * perct_ij**3._wp |
---|
2602 | ! |
---|
2603 | !-- Mean deviation plus its standard deviation |
---|
2604 | dperctstd = dperctsa + dperctsb * perct_ij + dperctsc * perct_ij**2._wp + & |
---|
2605 | dperctsd * perct_ij**3._wp |
---|
2606 | ! |
---|
2607 | !-- Value of the FUNCTION |
---|
2608 | sultr_res = dperctm + faktor * dperctstd |
---|
2609 | IF ( ABS ( sultr_res ) > 99._wp ) sultr_res = +99._wp |
---|
2610 | |
---|
2611 | END SUBROUTINE calc_sultr |
---|
2612 | |
---|
2613 | !------------------------------------------------------------------------------! |
---|
2614 | ! Description: |
---|
2615 | ! ------------ |
---|
2616 | !> Multiple linear regression to calculate an increment delta_cold, |
---|
2617 | !> to adjust Fanger's classical PMV (pmva) by Gagge's 2 node model, |
---|
2618 | !> applying Fanger's convective heat transfer coefficient, hcf. |
---|
2619 | !> Wind velocitiy of the reference environment is 0.10 m/s |
---|
2620 | !------------------------------------------------------------------------------! |
---|
2621 | SUBROUTINE dpmv_cold( pmva, ta, ws, tmrt, nerr, dpmv_cold_res ) |
---|
2622 | |
---|
2623 | IMPLICIT NONE |
---|
2624 | ! |
---|
2625 | !-- Type of input arguments |
---|
2626 | REAL(wp), INTENT ( IN ) :: pmva !< Fanger's classical predicted mean vote |
---|
2627 | REAL(wp), INTENT ( IN ) :: ta !< Air temperature 2 m above ground (degC) |
---|
2628 | REAL(wp), INTENT ( IN ) :: ws !< Relative wind velocity 1 m above ground (m/s) |
---|
2629 | REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (degC) |
---|
2630 | ! |
---|
2631 | !-- Type of output argument |
---|
2632 | INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error indicator: 0 = o.k., +1 = denominator for intersection = 0 |
---|
2633 | REAL(wp), INTENT ( OUT ) :: dpmv_cold_res !< Increment to adjust pmva according to the results of Gagge's |
---|
2634 | ! 2 node model depending on the input |
---|
2635 | ! |
---|
2636 | !-- Type of program variables |
---|
2637 | REAL(wp) :: delta_cold(3) |
---|
2638 | REAL(wp) :: pmv_cross(2) |
---|
2639 | REAL(wp) :: reg_a(3) |
---|
2640 | REAL(wp) :: coeff(3,5) |
---|
2641 | REAL(wp) :: r_nenner |
---|
2642 | REAL(wp) :: pmvc |
---|
2643 | REAL(wp) :: dtmrt |
---|
2644 | REAL(wp) :: sqrt_ws |
---|
2645 | INTEGER(iwp) :: i |
---|
2646 | ! INTEGER(iwp) :: j |
---|
2647 | INTEGER(iwp) :: i_bin |
---|
2648 | ! |
---|
2649 | !-- Coefficient of the 3 regression lines |
---|
2650 | ! 1:const 2:*pmvc 3:*ta 4:*sqrt_ws 5:*dtmrt |
---|
2651 | coeff(1,1:5) = & |
---|
2652 | (/ +0.161_wp, +0.130_wp, -1.125E-03_wp, +1.106E-03_wp, -4.570E-04_wp /) |
---|
2653 | coeff(2,1:5) = & |
---|
2654 | (/ 0.795_wp, 0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp/) |
---|
2655 | coeff(3,1:5) = & |
---|
2656 | (/ +0.05761_wp, +0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp /) |
---|
2657 | ! |
---|
2658 | !-- Initialise |
---|
2659 | nerr = 0_iwp |
---|
2660 | dpmv_cold_res = 0._wp |
---|
2661 | pmvc = pmva |
---|
2662 | dtmrt = tmrt - ta |
---|
2663 | sqrt_ws = ws |
---|
2664 | IF ( sqrt_ws < 0.10_wp ) THEN |
---|
2665 | sqrt_ws = 0.10_wp |
---|
2666 | ELSE |
---|
2667 | sqrt_ws = SQRT ( sqrt_ws ) |
---|
2668 | ENDIF |
---|
2669 | |
---|
2670 | DO i = 1, 3 |
---|
2671 | delta_cold (i) = 0._wp |
---|
2672 | IF ( i < 3 ) THEN |
---|
2673 | pmv_cross (i) = pmvc |
---|
2674 | ENDIF |
---|
2675 | ENDDO |
---|
2676 | |
---|
2677 | DO i = 1, 3 |
---|
2678 | ! |
---|
2679 | !-- Regression constant for given meteorological variables |
---|
2680 | reg_a(i) = coeff(i, 1) + coeff(i, 3) * ta + coeff(i, 4) * & |
---|
2681 | sqrt_ws + coeff(i,5)*dtmrt |
---|
2682 | delta_cold(i) = reg_a(i) + coeff(i, 2) * pmvc |
---|
2683 | ENDDO |
---|
2684 | ! |
---|
2685 | !-- Intersection points of regression lines in terms of Fanger's PMV |
---|
2686 | DO i = 1, 2 |
---|
2687 | r_nenner = coeff (i, 2) - coeff (i + 1, 2) |
---|
2688 | IF ( ABS ( r_nenner ) > 0.00001_wp ) THEN |
---|
2689 | pmv_cross(i) = ( reg_a (i + 1) - reg_a (i) ) / r_nenner |
---|
2690 | ELSE |
---|
2691 | nerr = 1_iwp |
---|
2692 | RETURN |
---|
2693 | ENDIF |
---|
2694 | ENDDO |
---|
2695 | |
---|
2696 | i_bin = 3 |
---|
2697 | DO i = 1, 2 |
---|
2698 | IF ( pmva > pmv_cross (i) ) THEN |
---|
2699 | i_bin = i |
---|
2700 | EXIT |
---|
2701 | ENDIF |
---|
2702 | ENDDO |
---|
2703 | ! |
---|
2704 | !-- Adjust to operative temperature scaled according |
---|
2705 | ! to classical PMV (Fanger) |
---|
2706 | dpmv_cold_res = delta_cold(i_bin) - dpmv_adj(pmva) |
---|
2707 | |
---|
2708 | END SUBROUTINE dpmv_cold |
---|
2709 | |
---|
2710 | !------------------------------------------------------------------------------! |
---|
2711 | ! Description: |
---|
2712 | ! ------------ |
---|
2713 | !> Calculates the summand dpmv_adj adjusting to the operative temperature |
---|
2714 | !> scaled according to classical PMV (Fanger) |
---|
2715 | !> Reference environment: v_1m = 0.10 m/s, dTMRT = 0 K, r.h. = 50 % |
---|
2716 | !------------------------------------------------------------------------------! |
---|
2717 | REAL(wp) FUNCTION dpmv_adj( pmva ) |
---|
2718 | |
---|
2719 | IMPLICIT NONE |
---|
2720 | |
---|
2721 | REAL(wp), INTENT ( IN ) :: pmva |
---|
2722 | INTEGER(iwp), PARAMETER :: n_bin = 3 |
---|
2723 | INTEGER(iwp), PARAMETER :: n_ca = 0 |
---|
2724 | INTEGER(iwp), PARAMETER :: n_ce = 3 |
---|
2725 | REAL(wp), dimension (n_bin, n_ca:n_ce) :: coef |
---|
2726 | |
---|
2727 | REAL(wp) :: pmv |
---|
2728 | INTEGER(iwp) :: i, i_bin |
---|
2729 | |
---|
2730 | ! range_1 range_2 range_3 |
---|
2731 | DATA (coef(i, 0), i = 1, n_bin) /0.0941540_wp, -0.1506620_wp, -0.0871439_wp/ |
---|
2732 | DATA (coef(i, 1), i = 1, n_bin) /0.0783162_wp, -1.0612651_wp, 0.1695040_wp/ |
---|
2733 | DATA (coef(i, 2), i = 1, n_bin) /0.1350144_wp, -1.0049144_wp, -0.0167627_wp/ |
---|
2734 | DATA (coef(i, 3), i = 1, n_bin) /0.1104037_wp, -0.2005277_wp, -0.0003230_wp/ |
---|
2735 | |
---|
2736 | IF ( pmva <= -2.1226_wp ) THEN |
---|
2737 | i_bin = 3_iwp |
---|
2738 | ELSE IF ( pmva <= -1.28_wp ) THEN |
---|
2739 | i_bin = 2_iwp |
---|
2740 | ELSE |
---|
2741 | i_bin = 1_iwp |
---|
2742 | ENDIF |
---|
2743 | |
---|
2744 | dpmv_adj = coef( i_bin, n_ca ) |
---|
2745 | pmv = 1._wp |
---|
2746 | |
---|
2747 | DO i = n_ca + 1, n_ce |
---|
2748 | pmv = pmv * pmva |
---|
2749 | dpmv_adj = dpmv_adj + coef(i_bin, i) * pmv |
---|
2750 | ENDDO |
---|
2751 | RETURN |
---|
2752 | END FUNCTION dpmv_adj |
---|
2753 | |
---|
2754 | !------------------------------------------------------------------------------! |
---|
2755 | ! Description: |
---|
2756 | ! ------------ |
---|
2757 | !> Based on perceived temperature (perct) as input, ireq_neutral determines |
---|
2758 | !> the required clothing insulation (clo) for thermally neutral conditions |
---|
2759 | !> (neither body cooling nor body heating). It is related to the Klima- |
---|
2760 | !> Michel activity level (134.682 W/m2). IREQ_neutral is only defined |
---|
2761 | !> for perct < 10 (degC) |
---|
2762 | !------------------------------------------------------------------------------! |
---|
2763 | REAL(wp) FUNCTION ireq_neutral( perct_ij, ireq_minimal, nerr ) |
---|
2764 | |
---|
2765 | IMPLICIT NONE |
---|
2766 | ! |
---|
2767 | !-- Type declaration of arguments |
---|
2768 | REAL(wp), INTENT ( IN ) :: perct_ij |
---|
2769 | REAL(wp), INTENT ( OUT ) :: ireq_minimal |
---|
2770 | INTEGER(iwp), INTENT ( OUT ) :: nerr |
---|
2771 | ! |
---|
2772 | !-- Type declaration for internal varables |
---|
2773 | REAL(wp) :: top02 |
---|
2774 | ! |
---|
2775 | !-- Initialise |
---|
2776 | nerr = 0_iwp |
---|
2777 | ! |
---|
2778 | !-- Convert perceived temperature from basis 0.1 m/s to basis 0.2 m/s |
---|
2779 | top02 = 1.8788_wp + 0.9296_wp * perct_ij |
---|
2780 | ! |
---|
2781 | !-- IREQ neutral conditions (thermal comfort) |
---|
2782 | ireq_neutral = 1.62_wp - 0.0564_wp * top02 |
---|
2783 | ! |
---|
2784 | !-- Regression only defined for perct <= 10 (degC) |
---|
2785 | IF ( ireq_neutral < 0.5_wp ) THEN |
---|
2786 | IF ( ireq_neutral < 0.48_wp ) THEN |
---|
2787 | nerr = 1_iwp |
---|
2788 | ENDIF |
---|
2789 | ireq_neutral = 0.5_wp |
---|
2790 | ENDIF |
---|
2791 | ! |
---|
2792 | !-- Minimal required clothing insulation: maximal acceptable body cooling |
---|
2793 | ireq_minimal = 1.26_wp - 0.0588_wp * top02 |
---|
2794 | IF ( nerr > 0_iwp ) THEN |
---|
2795 | ireq_minimal = ireq_neutral |
---|
2796 | ENDIF |
---|
2797 | |
---|
2798 | RETURN |
---|
2799 | END FUNCTION ireq_neutral |
---|
2800 | |
---|
2801 | |
---|
2802 | !------------------------------------------------------------------------------! |
---|
2803 | ! Description: |
---|
2804 | ! ------------ |
---|
2805 | !> The SUBROUTINE surface area calculates the surface area of the individual |
---|
2806 | !> according to its height (m), weight (kg), and age (y) |
---|
2807 | !------------------------------------------------------------------------------! |
---|
2808 | SUBROUTINE surface_area ( height_cm, weight, age, surf ) |
---|
2809 | |
---|
2810 | IMPLICIT NONE |
---|
2811 | |
---|
2812 | REAL(wp) , INTENT(in) :: weight |
---|
2813 | REAL(wp) , INTENT(in) :: height_cm |
---|
2814 | INTEGER(iwp), INTENT(in) :: age |
---|
2815 | REAL(wp) , INTENT(out) :: surf |
---|
2816 | REAL(wp) :: height |
---|
2817 | |
---|
2818 | height = height_cm * 100._wp |
---|
2819 | ! |
---|
2820 | !-- According to Gehan-George, for children |
---|
2821 | IF ( age < 19_iwp ) THEN |
---|
2822 | IF ( age < 5_iwp ) THEN |
---|
2823 | surf = 0.02667_wp * height**0.42246_wp * weight**0.51456_wp |
---|
2824 | RETURN |
---|
2825 | END IF |
---|
2826 | surf = 0.03050_wp * height**0.35129_wp * weight**0.54375_wp |
---|
2827 | RETURN |
---|
2828 | END IF |
---|
2829 | ! |
---|
2830 | !-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if |
---|
2831 | ! height and weight be known. In: Arch. Int. Med.. 17, 1916, S. 863?871. |
---|
2832 | surf = 0.007184_wp * height**0.725_wp * weight**0.425_wp |
---|
2833 | RETURN |
---|
2834 | |
---|
2835 | END SUBROUTINE surface_area |
---|
2836 | |
---|
2837 | !------------------------------------------------------------------------------! |
---|
2838 | ! Description: |
---|
2839 | ! ------------ |
---|
2840 | !> The SUBROUTINE persdat calculates |
---|
2841 | !> - the total internal energy production = metabolic + workload, |
---|
2842 | !> - the total internal energy production for a standardized surface (actlev) |
---|
2843 | !> - the DuBois - area (a_surf [m2]) |
---|
2844 | !> from |
---|
2845 | !> - the persons age (years), |
---|
2846 | !> - weight (kg), |
---|
2847 | !> - height (m), |
---|
2848 | !> - sex (1 = male, 2 = female), |
---|
2849 | !> - work load (W) |
---|
2850 | !> for a sample human. |
---|
2851 | !------------------------------------------------------------------------------! |
---|
2852 | SUBROUTINE persdat ( age, weight, height, sex, work, a_surf, actlev ) |
---|
2853 | ! |
---|
2854 | IMPLICIT NONE |
---|
2855 | |
---|
2856 | REAL(wp), INTENT(in) :: age |
---|
2857 | REAL(wp), INTENT(in) :: weight |
---|
2858 | REAL(wp), INTENT(in) :: height |
---|
2859 | REAL(wp), INTENT(in) :: work |
---|
2860 | INTEGER(iwp), INTENT(in) :: sex |
---|
2861 | REAL(wp), INTENT(out) :: actlev |
---|
2862 | REAL(wp) :: a_surf |
---|
2863 | REAL(wp) :: energy_prod |
---|
2864 | REAL(wp) :: s |
---|
2865 | REAL(wp) :: factor |
---|
2866 | REAL(wp) :: basic_heat_prod |
---|
2867 | |
---|
2868 | CALL surface_area ( height, weight, INT( age ), a_surf ) |
---|
2869 | s = height * 100._wp / ( weight ** ( 1._wp / 3._wp ) ) |
---|
2870 | factor = 1._wp + .004_wp * ( 30._wp - age ) |
---|
2871 | basic_heat_prod = 0. |
---|
2872 | IF ( sex == 1_iwp ) THEN |
---|
2873 | basic_heat_prod = 3.45_wp * weight ** ( 3._wp / 4._wp ) * ( factor + & |
---|
2874 | .01_wp * ( s - 43.4_wp ) ) |
---|
2875 | ELSE IF ( sex == 2_iwp ) THEN |
---|
2876 | basic_heat_prod = 3.19_wp * weight ** ( 3._wp / 4._wp ) * ( factor + & |
---|
2877 | .018_wp * ( s - 42.1_wp ) ) |
---|
2878 | END IF |
---|
2879 | |
---|
2880 | energy_prod = work + basic_heat_prod |
---|
2881 | actlev = energy_prod / a_surf |
---|
2882 | |
---|
2883 | END SUBROUTINE persdat |
---|
2884 | |
---|
2885 | |
---|
2886 | !------------------------------------------------------------------------------! |
---|
2887 | ! Description: |
---|
2888 | ! ------------ |
---|
2889 | !> SUBROUTINE ipt_init |
---|
2890 | !> initializes the instationary perceived temperature |
---|
2891 | !------------------------------------------------------------------------------! |
---|
2892 | |
---|
2893 | SUBROUTINE ipt_init (age, weight, height, sex, work, actlev, clo, & |
---|
2894 | ta, vp, ws, tmrt, pair, dt, storage, t_clothing, & |
---|
2895 | ipt ) |
---|
2896 | |
---|
2897 | IMPLICIT NONE |
---|
2898 | ! |
---|
2899 | !-- Input parameters |
---|
2900 | REAL(wp), INTENT(in) :: age !< Persons age (years) |
---|
2901 | REAL(wp), INTENT(in) :: weight !< Persons weight (kg) |
---|
2902 | REAL(wp), INTENT(in) :: height !< Persons height (m) |
---|
2903 | REAL(wp), INTENT(in) :: work !< Current workload (W) |
---|
2904 | REAL(wp), INTENT(in) :: ta !< Air Temperature (degree_C) |
---|
2905 | REAL(wp), INTENT(in) :: vp !< Vapor pressure (hPa) |
---|
2906 | REAL(wp), INTENT(in) :: ws !< Wind speed in approx. 1.1m (m/s) |
---|
2907 | REAL(wp), INTENT(in) :: tmrt !< Mean radiant temperature (degree_C) |
---|
2908 | REAL(wp), INTENT(in) :: pair !< Air pressure (hPa) |
---|
2909 | REAL(wp), INTENT(in) :: dt !< Timestep (s) |
---|
2910 | INTEGER(iwp), INTENT(in) :: sex !< Persons sex (1 = male, 2 = female) |
---|
2911 | ! |
---|
2912 | !-- Output parameters |
---|
2913 | REAL(wp), INTENT(out) :: actlev |
---|
2914 | REAL(wp), INTENT(out) :: clo |
---|
2915 | REAL(wp), INTENT(out) :: storage |
---|
2916 | REAL(wp), INTENT(out) :: t_clothing |
---|
2917 | REAL(wp), INTENT(out) :: ipt |
---|
2918 | ! |
---|
2919 | !-- Internal variables |
---|
2920 | REAL(wp), PARAMETER :: eps = 0.0005_wp |
---|
2921 | REAL(wp), PARAMETER :: eta = 0._wp |
---|
2922 | REAL(wp) :: sclo |
---|
2923 | REAL(wp) :: wclo |
---|
2924 | REAL(wp) :: d_pmv |
---|
2925 | REAL(wp) :: svp_ta |
---|
2926 | REAL(wp) :: sult_lim |
---|
2927 | REAL(wp) :: dgtcm |
---|
2928 | REAL(wp) :: dgtcstd |
---|
2929 | REAL(wp) :: clon |
---|
2930 | REAL(wp) :: ireq_minimal |
---|
2931 | ! REAL(wp) :: clo_fanger |
---|
2932 | REAL(wp) :: pmv_w |
---|
2933 | REAL(wp) :: pmv_s |
---|
2934 | REAL(wp) :: pmva |
---|
2935 | REAL(wp) :: ptc |
---|
2936 | REAL(wp) :: d_std |
---|
2937 | REAL(wp) :: pmvs |
---|
2938 | REAL(wp) :: top |
---|
2939 | REAL(wp) :: a_surf |
---|
2940 | ! REAL(wp) :: acti |
---|
2941 | INTEGER(iwp) :: ncount |
---|
2942 | INTEGER(iwp) :: nerr_cold |
---|
2943 | INTEGER(iwp) :: nerr |
---|
2944 | |
---|
2945 | LOGICAL :: sultrieness |
---|
2946 | |
---|
2947 | storage = 0._wp |
---|
2948 | CALL persdat ( age, weight, height, sex, work, a_surf, actlev ) |
---|
2949 | ! |
---|
2950 | !-- Initialise |
---|
2951 | t_clothing = bio_fill_value |
---|
2952 | ipt = bio_fill_value |
---|
2953 | nerr = 0_wp |
---|
2954 | ncount = 0_wp |
---|
2955 | sultrieness = .FALSE. |
---|
2956 | ! |
---|
2957 | !-- Tresholds: clothing insulation (account for model inaccuracies) |
---|
2958 | !-- Summer clothing |
---|
2959 | sclo = 0.44453_wp |
---|
2960 | !-- Winter clothing |
---|
2961 | wclo = 1.76267_wp |
---|
2962 | ! |
---|
2963 | !-- Decision: firstly calculate for winter or summer clothing |
---|
2964 | IF ( ta <= 10._wp ) THEN |
---|
2965 | ! |
---|
2966 | !-- First guess: winter clothing insulation: cold stress |
---|
2967 | clo = wclo |
---|
2968 | t_clothing = bio_fill_value ! force initial run |
---|
2969 | CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, & |
---|
2970 | t_clothing, storage, dt, pmva ) |
---|
2971 | pmv_w = pmva |
---|
2972 | |
---|
2973 | IF ( pmva > 0._wp ) THEN |
---|
2974 | ! |
---|
2975 | !-- Case summer clothing insulation: heat load ? |
---|
2976 | clo = sclo |
---|
2977 | t_clothing = bio_fill_value ! force initial run |
---|
2978 | CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, & |
---|
2979 | t_clothing, storage, dt, pmva ) |
---|
2980 | pmv_s = pmva |
---|
2981 | IF ( pmva <= 0._wp ) THEN |
---|
2982 | ! |
---|
2983 | !-- Case: comfort achievable by varying clothing insulation |
---|
2984 | ! between winter and summer set values |
---|
2985 | CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo, & |
---|
2986 | pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo ) |
---|
2987 | IF ( ncount < 0_iwp ) THEN |
---|
2988 | nerr = -1_iwp |
---|
2989 | RETURN |
---|
2990 | ENDIF |
---|
2991 | ELSE IF ( pmva > 0.06_wp ) THEN |
---|
2992 | clo = 0.5_wp |
---|
2993 | t_clothing = bio_fill_value |
---|
2994 | CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, & |
---|
2995 | t_clothing, storage, dt, pmva ) |
---|
2996 | ENDIF |
---|
2997 | ELSE IF ( pmva < -0.11_wp ) THEN |
---|
2998 | clo = 1.75_wp |
---|
2999 | t_clothing = bio_fill_value |
---|
3000 | CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, & |
---|
3001 | t_clothing, storage, dt, pmva ) |
---|
3002 | ENDIF |
---|
3003 | |
---|
3004 | ELSE |
---|
3005 | ! |
---|
3006 | !-- First guess: summer clothing insulation: heat load |
---|
3007 | clo = sclo |
---|
3008 | t_clothing = bio_fill_value |
---|
3009 | CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, & |
---|
3010 | t_clothing, storage, dt, pmva ) |
---|
3011 | pmv_s = pmva |
---|
3012 | |
---|
3013 | IF ( pmva < 0._wp ) THEN |
---|
3014 | ! |
---|
3015 | !-- Case winter clothing insulation: cold stress ? |
---|
3016 | clo = wclo |
---|
3017 | t_clothing = bio_fill_value |
---|
3018 | CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, & |
---|
3019 | t_clothing, storage, dt, pmva ) |
---|
3020 | pmv_w = pmva |
---|
3021 | |
---|
3022 | IF ( pmva >= 0._wp ) THEN |
---|
3023 | ! |
---|
3024 | !-- Case: comfort achievable by varying clothing insulation |
---|
3025 | ! between winter and summer set values |
---|
3026 | CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, & |
---|
3027 | pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo ) |
---|
3028 | IF ( ncount < 0_wp ) THEN |
---|
3029 | nerr = -1_iwp |
---|
3030 | RETURN |
---|
3031 | ENDIF |
---|
3032 | ELSE IF ( pmva < -0.11_wp ) THEN |
---|
3033 | clo = 1.75_wp |
---|
3034 | t_clothing = bio_fill_value |
---|
3035 | CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, & |
---|
3036 | t_clothing, storage, dt, pmva ) |
---|
3037 | ENDIF |
---|
3038 | ELSE IF ( pmva > 0.06_wp ) THEN |
---|
3039 | clo = 0.5_wp |
---|
3040 | t_clothing = bio_fill_value |
---|
3041 | CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, & |
---|
3042 | t_clothing, storage, dt, pmva ) |
---|
3043 | ENDIF |
---|
3044 | |
---|
3045 | ENDIF |
---|
3046 | ! |
---|
3047 | !-- Determine perceived temperature by regression equation + adjustments |
---|
3048 | pmvs = pmva |
---|
3049 | CALL perct_regression ( pmva, clo, ipt ) |
---|
3050 | ptc = ipt |
---|
3051 | IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN |
---|
3052 | ! |
---|
3053 | !-- Adjust for cold conditions according to Gagge 1986 |
---|
3054 | CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv ) |
---|
3055 | IF ( nerr_cold > 0_iwp ) nerr = -5_iwp |
---|
3056 | pmvs = pmva - d_pmv |
---|
3057 | IF ( pmvs > -0.11_wp ) THEN |
---|
3058 | d_pmv = 0._wp |
---|
3059 | pmvs = -0.11_wp |
---|
3060 | ENDIF |
---|
3061 | CALL perct_regression ( pmvs, clo, ipt ) |
---|
3062 | ENDIF |
---|
3063 | ! clo_fanger = clo |
---|
3064 | clon = clo |
---|
3065 | IF ( clo > 0.5_wp .AND. ipt <= 8.73_wp ) THEN |
---|
3066 | ! |
---|
3067 | !-- Required clothing insulation (ireq) is exclusively defined for |
---|
3068 | ! operative temperatures (top) less 10 (C) for a |
---|
3069 | ! reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s |
---|
3070 | clon = ireq_neutral ( ipt, ireq_minimal, nerr ) |
---|
3071 | clo = clon |
---|
3072 | ENDIF |
---|
3073 | CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim ) |
---|
3074 | sultrieness = .FALSE. |
---|
3075 | d_std = -99._wp |
---|
3076 | IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN |
---|
3077 | ! |
---|
3078 | !-- Adjust for warm/humid conditions according to Gagge 1986 |
---|
3079 | CALL saturation_vapor_pressure ( ta, svp_ta ) |
---|
3080 | d_pmv = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr ) |
---|
3081 | pmvs = pmva + d_pmv |
---|
3082 | CALL perct_regression ( pmvs, clo, ipt ) |
---|
3083 | IF ( sult_lim < 99._wp ) THEN |
---|
3084 | IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE. |
---|
3085 | ENDIF |
---|
3086 | ENDIF |
---|
3087 | |
---|
3088 | |
---|
3089 | END SUBROUTINE ipt_init |
---|
3090 | |
---|
3091 | !------------------------------------------------------------------------------! |
---|
3092 | ! Description: |
---|
3093 | ! ------------ |
---|
3094 | !> SUBROUTINE ipt_cycle |
---|
3095 | !> Calculates one timestep for the instationary version of perceived |
---|
3096 | !> temperature (iPT, degree_C) for |
---|
3097 | !> - standard measured/predicted meteorological values and TMRT |
---|
3098 | !> as input; |
---|
3099 | !> - regressions for determination of PT; |
---|
3100 | !> - adjustment to Gagge's PMV* (2-node-model, 1986) as base of PT |
---|
3101 | !> under warm/humid conditions (Icl= 0.50 clo) and under cold |
---|
3102 | !> conditions (Icl= 1.75 clo) |
---|
3103 | !> |
---|
3104 | !------------------------------------------------------------------------------! |
---|
3105 | SUBROUTINE ipt_cycle( ta, vp, ws, tmrt, pair, dt, storage, t_clothing, clo, & |
---|
3106 | actlev, work, ipt ) |
---|
3107 | |
---|
3108 | IMPLICIT NONE |
---|
3109 | ! |
---|
3110 | !-- Type of input of the argument list |
---|
3111 | REAL(wp), INTENT ( IN ) :: ta !< Air temperature (degree_C) |
---|
3112 | REAL(wp), INTENT ( IN ) :: vp !< Vapor pressure (hPa) |
---|
3113 | REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (degree_C) |
---|
3114 | REAL(wp), INTENT ( IN ) :: ws !< Wind speed (m/s) |
---|
3115 | REAL(wp), INTENT ( IN ) :: pair !< Air pressure (hPa) |
---|
3116 | REAL(wp), INTENT ( IN ) :: dt !< Timestep (s) |
---|
3117 | REAL(wp), INTENT ( IN ) :: clo !< Clothing index (no dim) |
---|
3118 | REAL(wp), INTENT ( IN ) :: actlev !< Internal heat production (W) |
---|
3119 | REAL(wp), INTENT ( IN ) :: work !< Mechanical work load (W) |
---|
3120 | ! |
---|
3121 | !-- In and output parameters |
---|
3122 | REAL(wp), INTENT (INOUT) :: storage !< Heat storage (W) |
---|
3123 | REAL(wp), INTENT (INOUT) :: t_clothing !< Clothig temperature (degree_C) |
---|
3124 | ! |
---|
3125 | !-- Type of output of the argument list |
---|
3126 | REAL(wp), INTENT ( OUT ) :: ipt !< Instationary perceived temperature (degree_C) |
---|
3127 | ! |
---|
3128 | !-- Type of internal variables |
---|
3129 | REAL(wp) :: d_pmv |
---|
3130 | REAL(wp) :: svp_ta |
---|
3131 | REAL(wp) :: sult_lim |
---|
3132 | REAL(wp) :: dgtcm |
---|
3133 | REAL(wp) :: dgtcstd |
---|
3134 | REAL(wp) :: pmva |
---|
3135 | REAL(wp) :: ptc |
---|
3136 | REAL(wp) :: d_std |
---|
3137 | REAL(wp) :: pmvs |
---|
3138 | INTEGER(iwp) :: nerr_cold |
---|
3139 | INTEGER(iwp) :: nerr |
---|
3140 | |
---|
3141 | LOGICAL :: sultrieness |
---|
3142 | ! |
---|
3143 | !-- Initialise |
---|
3144 | ipt = bio_fill_value |
---|
3145 | |
---|
3146 | nerr = 0_iwp |
---|
3147 | sultrieness = .FALSE. |
---|
3148 | ! |
---|
3149 | !-- Determine pmv_adjusted for current conditions |
---|
3150 | CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, & |
---|
3151 | t_clothing, storage, dt, pmva ) |
---|
3152 | ! |
---|
3153 | !-- Determine perceived temperature by regression equation + adjustments |
---|
3154 | CALL perct_regression ( pmva, clo, ipt ) |
---|
3155 | |
---|
3156 | IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN |
---|
3157 | ! |
---|
3158 | !-- Adjust for cold conditions according to Gagge 1986 |
---|
3159 | CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv ) |
---|
3160 | IF ( nerr_cold > 0_iwp ) nerr = -5_iwp |
---|
3161 | pmvs = pmva - d_pmv |
---|
3162 | IF ( pmvs > -0.11_wp ) THEN |
---|
3163 | d_pmv = 0._wp |
---|
3164 | pmvs = -0.11_wp |
---|
3165 | ENDIF |
---|
3166 | CALL perct_regression ( pmvs, clo, ipt ) |
---|
3167 | ENDIF |
---|
3168 | ! |
---|
3169 | !-- Consider sultriness if appropriate |
---|
3170 | ptc = ipt |
---|
3171 | CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim ) |
---|
3172 | sultrieness = .FALSE. |
---|
3173 | d_std = -99._wp |
---|
3174 | IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN |
---|
3175 | ! |
---|
3176 | !-- Adjust for warm/humid conditions according to Gagge 1986 |
---|
3177 | CALL saturation_vapor_pressure ( ta, svp_ta ) |
---|
3178 | d_pmv = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr ) |
---|
3179 | pmvs = pmva + d_pmv |
---|
3180 | CALL perct_regression ( pmvs, clo, ipt ) |
---|
3181 | IF ( sult_lim < 99._wp ) THEN |
---|
3182 | IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE. |
---|
3183 | ENDIF |
---|
3184 | ENDIF |
---|
3185 | |
---|
3186 | END SUBROUTINE ipt_cycle |
---|
3187 | |
---|
3188 | !------------------------------------------------------------------------------! |
---|
3189 | ! Description: |
---|
3190 | ! ------------ |
---|
3191 | !> SUBROUTINE fanger_s calculates the |
---|
3192 | !> actual Predicted Mean Vote (dimensionless) according |
---|
3193 | !> to Fanger corresponding to meteorological (ta,tmrt,pa,ws,pair) |
---|
3194 | !> and individual variables (clo, actlev, eta) considering a storage |
---|
3195 | !> and clothing temperature for a given timestep. |
---|
3196 | !------------------------------------------------------------------------------! |
---|
3197 | SUBROUTINE fanger_s_acti( ta, tmrt, pa, in_ws, pair, in_clo, actlev, & |
---|
3198 | activity, t_cloth, s, dt, pmva ) |
---|
3199 | |
---|
3200 | IMPLICIT NONE |
---|
3201 | ! |
---|
3202 | !-- Input argument types |
---|
3203 | REAL(wp), INTENT ( IN ) :: ta !< Air temperature (degree_C) |
---|
3204 | REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (degree_C) |
---|
3205 | REAL(wp), INTENT ( IN ) :: pa !< Vapour pressure (hPa) |
---|
3206 | REAL(wp), INTENT ( IN ) :: pair !< Air pressure (hPa) |
---|
3207 | REAL(wp), INTENT ( IN ) :: in_ws !< Wind speed (m/s) |
---|
3208 | REAL(wp), INTENT ( IN ) :: actlev !< Metabolic + work energy (W/m²) |
---|
3209 | REAL(wp), INTENT ( IN ) :: dt !< Timestep (s) |
---|
3210 | REAL(wp), INTENT ( IN ) :: activity !< Work load (W/m²) |
---|
3211 | REAL(wp), INTENT ( IN ) :: in_clo !< Clothing index (clo) (no dim) |
---|
3212 | ! |
---|
3213 | !-- Output argument types |
---|
3214 | REAL(wp), INTENT ( OUT ) :: pmva !< actual Predicted Mean Vote (no dim) |
---|
3215 | |
---|
3216 | REAL(wp), INTENT (INOUT) :: s !< storage var. of energy balance (W/m2) |
---|
3217 | REAL(wp), INTENT (INOUT) :: t_cloth !< clothing temperature (degree_C) |
---|
3218 | ! |
---|
3219 | !-- Internal variables |
---|
3220 | REAL(wp), PARAMETER :: time_equil = 7200._wp |
---|
3221 | |
---|
3222 | REAL(wp) :: f_cl !< Increase in surface due to clothing (factor) |
---|
3223 | REAL(wp) :: heat_convection !< energy loss by autocnvection (W) |
---|
3224 | REAL(wp) :: t_skin_aver !< average skin temperature (degree_C) |
---|
3225 | REAL(wp) :: bc !< preliminary result storage |
---|
3226 | REAL(wp) :: cc !< preliminary result storage |
---|
3227 | REAL(wp) :: dc !< preliminary result storage |
---|
3228 | REAL(wp) :: ec !< preliminary result storage |
---|
3229 | REAL(wp) :: gc !< preliminary result storage |
---|
3230 | REAL(wp) :: t_clothing !< clothing temperature (degree_C) |
---|
3231 | ! REAL(wp) :: hr !< radiational heat resistence |
---|
3232 | REAL(wp) :: clo !< clothing insulation index (clo) |
---|
3233 | REAL(wp) :: ws !< wind speed (m/s) |
---|
3234 | REAL(wp) :: z1 !< Empiric factor for the adaption of the heat |
---|
3235 | ! ballance equation to the psycho-physical scale (Equ. 40 in FANGER) |
---|
3236 | REAL(wp) :: z2 !< Water vapour diffution through the skin |
---|
3237 | REAL(wp) :: z3 !< Sweat evaporation from the skin surface |
---|
3238 | REAL(wp) :: z4 !< Loss of latent heat through respiration |
---|
3239 | REAL(wp) :: z5 !< Loss of radiational heat |
---|
3240 | REAL(wp) :: z6 !< Heat loss through forced convection |
---|
3241 | REAL(wp) :: en !< Energy ballance (W) |
---|
3242 | REAL(wp) :: d_s !< Storage delta (W) |
---|
3243 | REAL(wp) :: adjustrate !< Max storage adjustment rate |
---|
3244 | REAL(wp) :: adjustrate_cloth !< max clothing temp. adjustment rate |
---|
3245 | |
---|
3246 | INTEGER(iwp) :: i !< running index |
---|
3247 | INTEGER(iwp) :: niter !< Running index |
---|
3248 | |
---|
3249 | ! |
---|
3250 | !-- Clo must be > 0. to avoid div. by 0! |
---|
3251 | clo = in_clo |
---|
3252 | IF ( clo < 001._wp ) clo = .001_wp |
---|
3253 | ! |
---|
3254 | !-- Increase in surface due to clothing |
---|
3255 | f_cl = 1._wp + .15_wp * clo |
---|
3256 | ! |
---|
3257 | !-- Case of free convection (ws < 0.1 m/s ) not considered |
---|
3258 | ws = in_ws |
---|
3259 | IF ( ws < .1_wp ) THEN |
---|
3260 | ws = .1_wp |
---|
3261 | ENDIF |
---|
3262 | ! |
---|
3263 | !-- Heat_convection = forced convection |
---|
3264 | heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp ) |
---|
3265 | ! |
---|
3266 | !-- Average skin temperature |
---|
3267 | t_skin_aver = 35.7_wp - .0275_wp * activity |
---|
3268 | ! |
---|
3269 | !-- Calculation of constants for evaluation below |
---|
3270 | bc = .155_wp * clo * 3.96_wp * 10._wp**( -8._wp ) * f_cl |
---|
3271 | cc = f_cl * heat_convection |
---|
3272 | ec = .155_wp * clo |
---|
3273 | dc = ( 1._wp + ec * cc ) / bc |
---|
3274 | gc = ( t_skin_aver + bc * ( tmrt + 273.2_wp )**4._wp + ec * cc * ta ) / bc |
---|
3275 | ! |
---|
3276 | !-- Calculation of clothing surface temperature (t_clothing) based on |
---|
3277 | ! newton-approximation with air temperature as initial guess |
---|
3278 | niter = INT( dt * 10._wp, KIND=iwp ) |
---|
3279 | IF ( niter < 1 ) niter = 1_iwp |
---|
3280 | adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt ) |
---|
3281 | IF ( adjustrate >= 1._wp ) adjustrate = 1._wp |
---|
3282 | adjustrate_cloth = adjustrate * 30._wp |
---|
3283 | t_clothing = t_cloth |
---|
3284 | |
---|
3285 | IF ( t_cloth <= -998.0_wp ) THEN ! If initial run |
---|
3286 | niter = 3_iwp |
---|
3287 | adjustrate = 1._wp |
---|
3288 | adjustrate_cloth = 1._wp |
---|
3289 | t_clothing = ta |
---|
3290 | ENDIF |
---|
3291 | |
---|
3292 | DO i = 1, niter |
---|
3293 | t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing + & |
---|
3294 | 273.2_wp )**4._wp + t_clothing * & |
---|
3295 | dc - gc ) / ( 4._wp * ( t_clothing + 273.2_wp )**3._wp + dc ) |
---|
3296 | ENDDO |
---|
3297 | ! |
---|
3298 | !-- Empiric factor for the adaption of the heat ballance equation |
---|
3299 | ! to the psycho-physical scale (Equ. 40 in FANGER) |
---|
3300 | z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp ) |
---|
3301 | ! |
---|
3302 | !-- Water vapour diffution through the skin |
---|
3303 | z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa ) |
---|
3304 | ! |
---|
3305 | !-- Sweat evaporation from the skin surface |
---|
3306 | z3 = .42_wp * ( activity - 58._wp ) |
---|
3307 | ! |
---|
3308 | !-- Loss of latent heat through respiration |
---|
3309 | z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev * & |
---|
3310 | ( 34._wp - ta ) |
---|
3311 | ! |
---|
3312 | !-- Loss of radiational heat |
---|
3313 | z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + 273.2_wp )**4 - ( tmrt + & |
---|
3314 | 273.2_wp )**4 ) |
---|
3315 | ! |
---|
3316 | !-- Heat loss through forced convection |
---|
3317 | z6 = cc * ( t_clothing - ta ) |
---|
3318 | ! |
---|
3319 | !-- Write together as energy ballance |
---|
3320 | en = activity - z2 - z3 - z4 - z5 - z6 |
---|
3321 | ! |
---|
3322 | !-- Manage storage |
---|
3323 | d_s = adjustrate * en + ( 1._wp - adjustrate ) * s |
---|
3324 | ! |
---|
3325 | !-- Predicted Mean Vote |
---|
3326 | pmva = z1 * d_s |
---|
3327 | ! |
---|
3328 | !-- Update storage |
---|
3329 | s = d_s |
---|
3330 | t_cloth = t_clothing |
---|
3331 | |
---|
3332 | END SUBROUTINE fanger_s_acti |
---|
3333 | |
---|
3334 | |
---|
3335 | |
---|
3336 | !------------------------------------------------------------------------------! |
---|
3337 | ! |
---|
3338 | ! Description: |
---|
3339 | ! ------------ |
---|
3340 | !> Physiologically Equivalent Temperature (PET), |
---|
3341 | !> stationary (calculated based on MEMI), |
---|
3342 | !> Subroutine based on PETBER vers. 1.5.1996 by P. Hoeppe |
---|
3343 | !------------------------------------------------------------------------------! |
---|
3344 | |
---|
3345 | SUBROUTINE calculate_pet_static( ta, vpa, v, tmrt, pair, pet_ij ) |
---|
3346 | |
---|
3347 | IMPLICIT NONE |
---|
3348 | ! |
---|
3349 | !-- Input arguments: |
---|
3350 | REAL(wp), INTENT( IN ) :: ta !< Air temperature (degree_C) |
---|
3351 | REAL(wp), INTENT( IN ) :: tmrt !< Mean radiant temperature (degree_C) |
---|
3352 | REAL(wp), INTENT( IN ) :: v !< Wind speed (m/s) |
---|
3353 | REAL(wp), INTENT( IN ) :: vpa !< Vapor pressure (hPa) |
---|
3354 | REAL(wp), INTENT( IN ) :: pair !< Air pressure (hPa) |
---|
3355 | ! |
---|
3356 | !-- Output arguments: |
---|
3357 | REAL(wp), INTENT ( OUT ) :: pet_ij !< PET (degree_C) |
---|
3358 | ! |
---|
3359 | !-- Internal variables: |
---|
3360 | REAL(wp) :: acl !< clothing area (m²) |
---|
3361 | REAL(wp) :: adu !< Du Bois area (m²) |
---|
3362 | REAL(wp) :: aeff !< effective area (m²) |
---|
3363 | REAL(wp) :: ere !< energy ballance (W) |
---|
3364 | REAL(wp) :: erel !< latent energy ballance (W) |
---|
3365 | REAL(wp) :: esw !< Energy-loss through sweat evap. (W) |
---|
3366 | REAL(wp) :: facl !< Surface area extension through clothing (factor) |
---|
3367 | REAL(wp) :: feff !< Surface modification by posture (factor) |
---|
3368 | REAL(wp) :: rdcl !< Diffusion resistence of clothing (factor) |
---|
3369 | REAL(wp) :: rdsk !< Diffusion resistence of skin (factor) |
---|
3370 | REAL(wp) :: rtv |
---|
3371 | REAL(wp) :: vpts !< Sat. vapor pressure over skin (hPa) |
---|
3372 | REAL(wp) :: tsk !< Skin temperature (degree_C) |
---|
3373 | REAL(wp) :: tcl !< Clothing temperature (degree_C) |
---|
3374 | REAL(wp) :: wetsk !< Fraction of wet skin (factor) |
---|
3375 | ! |
---|
3376 | !-- Variables: |
---|
3377 | REAL(wp) :: int_heat !< Internal heat (W) |
---|
3378 | ! |
---|
3379 | !-- MEMI configuration |
---|
3380 | REAL(wp) :: age !< Persons age (a) |
---|
3381 | REAL(wp) :: mbody !< Persons body mass (kg) |
---|
3382 | REAL(wp) :: ht !< Persons height (m) |
---|
3383 | REAL(wp) :: work !< Work load (W) |
---|
3384 | REAL(wp) :: eta !< Work efficiency (dimensionless) |
---|
3385 | REAL(wp) :: clo !< Clothing insulation index (clo) |
---|
3386 | REAL(wp) :: fcl !< Surface area modification by clothing (factor) |
---|
3387 | ! INTEGER(iwp) :: pos !< Posture: 1 = standing, 2 = sitting |
---|
3388 | ! INTEGER(iwp) :: sex !< Sex: 1 = male, 2 = female |
---|
3389 | ! |
---|
3390 | !-- Configuration, keep standard parameters! |
---|
3391 | age = 35._wp |
---|
3392 | mbody = 75._wp |
---|
3393 | ht = 1.75_wp |
---|
3394 | work = 80._wp |
---|
3395 | eta = 0._wp |
---|
3396 | clo = 0.9_wp |
---|
3397 | fcl = 1.15_wp |
---|
3398 | ! |
---|
3399 | !-- Call subfunctions |
---|
3400 | CALL in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, & |
---|
3401 | vpa, work ) |
---|
3402 | |
---|
3403 | CALL heat_exch ( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht, & |
---|
3404 | int_heat,mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, & |
---|
3405 | vpts, wetsk ) |
---|
3406 | |
---|
3407 | CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair, & |
---|
3408 | rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk ) |
---|
3409 | |
---|
3410 | |
---|
3411 | END SUBROUTINE calculate_pet_static |
---|
3412 | |
---|
3413 | |
---|
3414 | !------------------------------------------------------------------------------! |
---|
3415 | ! Description: |
---|
3416 | ! ------------ |
---|
3417 | !> Calculate internal energy ballance |
---|
3418 | !------------------------------------------------------------------------------! |
---|
3419 | SUBROUTINE in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, & |
---|
3420 | vpa, work ) |
---|
3421 | ! |
---|
3422 | !-- Input arguments: |
---|
3423 | REAL(wp), INTENT( IN ) :: pair !< air pressure (hPa) |
---|
3424 | REAL(wp), INTENT( IN ) :: ta !< air temperature (degree_C) |
---|
3425 | REAL(wp), INTENT( IN ) :: vpa !< vapor pressure (hPa) |
---|
3426 | REAL(wp), INTENT( IN ) :: age !< Persons age (a) |
---|
3427 | REAL(wp), INTENT( IN ) :: mbody !< Persons body mass (kg) |
---|
3428 | REAL(wp), INTENT( IN ) :: ht !< Persons height (m) |
---|
3429 | REAL(wp), INTENT( IN ) :: work !< Work load (W) |
---|
3430 | REAL(wp), INTENT( IN ) :: eta !< Work efficiency (dimensionless) |
---|
3431 | ! |
---|
3432 | !-- Output arguments: |
---|
3433 | REAL(wp), INTENT( OUT ) :: ere !< energy ballance (W) |
---|
3434 | REAL(wp), INTENT( OUT ) :: erel !< latent energy ballance (W) |
---|
3435 | REAL(wp), INTENT( OUT ) :: int_heat !< internal heat production (W) |
---|
3436 | REAL(wp), INTENT( OUT ) :: rtv !< respiratory volume |
---|
3437 | ! |
---|
3438 | !-- Constants: |
---|
3439 | ! REAL(wp), PARAMETER :: cair = 1010._wp !< replaced by c_p |
---|
3440 | ! REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp !< replaced by l_v |
---|
3441 | ! |
---|
3442 | !-- Internal variables: |
---|
3443 | REAL(wp) :: eres !< Sensible respiratory heat flux (W) |
---|
3444 | REAL(wp) :: met |
---|
3445 | REAL(wp) :: tex |
---|
3446 | REAL(wp) :: vpex |
---|
3447 | |
---|
3448 | met = 3.45_wp * mbody ** ( 3._wp / 4._wp ) * (1._wp + 0.004_wp * & |
---|
3449 | ( 30._wp - age) + 0.010_wp * ( ( ht * 100._wp / & |
---|
3450 | ( mbody ** ( 1._wp / 3._wp ) ) ) - 43.4_wp ) ) |
---|
3451 | |
---|
3452 | met = work + met |
---|
3453 | |
---|
3454 | int_heat = met * (1._wp - eta) |
---|
3455 | ! |
---|
3456 | !-- Sensible respiration energy |
---|
3457 | tex = 0.47_wp * ta + 21.0_wp |
---|
3458 | rtv = 1.44_wp * 10._wp ** (-6._wp) * met |
---|
3459 | eres = c_p * (ta - tex) * rtv |
---|
3460 | ! |
---|
3461 | !-- Latent respiration energy |
---|
3462 | vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) ) |
---|
3463 | erel = 0.623_wp * l_v / pair * ( vpa - vpex ) * rtv |
---|
3464 | ! |
---|
3465 | !-- Sum of the results |
---|
3466 | ere = eres + erel |
---|
3467 | |
---|
3468 | END SUBROUTINE in_body |
---|
3469 | |
---|
3470 | |
---|
3471 | !------------------------------------------------------------------------------! |
---|
3472 | ! Description: |
---|
3473 | ! ------------ |
---|
3474 | !> Calculate heat gain or loss |
---|
3475 | !------------------------------------------------------------------------------! |
---|
3476 | SUBROUTINE heat_exch ( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, & |
---|
3477 | ht, int_heat, mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, & |
---|
3478 | vpts, wetsk ) |
---|
3479 | |
---|
3480 | ! |
---|
3481 | !-- Input arguments: |
---|
3482 | REAL(wp), INTENT( IN ) :: ere !< Energy ballance (W) |
---|
3483 | REAL(wp), INTENT( IN ) :: erel !< Latent energy ballance (W) |
---|
3484 | REAL(wp), INTENT( IN ) :: int_heat !< internal heat production (W) |
---|
3485 | REAL(wp), INTENT( IN ) :: pair !< Air pressure (hPa) |
---|
3486 | REAL(wp), INTENT( IN ) :: ta !< Air temperature (degree_C) |
---|
3487 | REAL(wp), INTENT( IN ) :: tmrt !< Mean radiant temperature (degree_C) |
---|
3488 | REAL(wp), INTENT( IN ) :: v !< Wind speed (m/s) |
---|
3489 | REAL(wp), INTENT( IN ) :: vpa !< Vapor pressure (hPa) |
---|
3490 | REAL(wp), INTENT( IN ) :: mbody !< body mass (kg) |
---|
3491 | REAL(wp), INTENT( IN ) :: ht !< height (m) |
---|
3492 | REAL(wp), INTENT( IN ) :: clo !< clothing insulation (clo) |
---|
3493 | REAL(wp), INTENT( IN ) :: fcl !< factor for surface area increase by clothing |
---|
3494 | ! |
---|
3495 | !-- Output arguments: |
---|
3496 | REAL(wp), INTENT( OUT ) :: acl !< Clothing surface area (m²) |
---|
3497 | REAL(wp), INTENT( OUT ) :: adu !< Du-Bois area (m²) |
---|
3498 | REAL(wp), INTENT( OUT ) :: aeff !< Effective surface area (m²) |
---|
3499 | REAL(wp), INTENT( OUT ) :: esw !< Energy-loss through sweat evap. (W) |
---|
3500 | REAL(wp), INTENT( OUT ) :: facl !< Surface area extension through clothing (factor) |
---|
3501 | REAL(wp), INTENT( OUT ) :: feff !< Surface modification by posture (factor) |
---|
3502 | REAL(wp), INTENT( OUT ) :: rdcl !< Diffusion resistence of clothing (factor) |
---|
3503 | REAL(wp), INTENT( OUT ) :: rdsk !< Diffusion resistence of skin (factor) |
---|
3504 | REAL(wp), INTENT( OUT ) :: tcl !< Clothing temperature (degree_C) |
---|
3505 | REAL(wp), INTENT( OUT ) :: tsk !< Skin temperature (degree_C) |
---|
3506 | REAL(wp), INTENT( OUT ) :: vpts !< Sat. vapor pressure over skin (hPa) |
---|
3507 | REAL(wp), INTENT( OUT ) :: wetsk !< Fraction of wet skin (dimensionless) |
---|
3508 | ! |
---|
3509 | !-- Cconstants: |
---|
3510 | ! REAL(wp), PARAMETER :: cair = 1010._wp !< replaced by c_p |
---|
3511 | REAL(wp), PARAMETER :: cb = 3640._wp !< |
---|
3512 | REAL(wp), PARAMETER :: emcl = 0.95_wp !< Longwave emission coef. of cloth |
---|
3513 | REAL(wp), PARAMETER :: emsk = 0.99_wp !< Longwave emission coef. of skin |
---|
3514 | ! REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp !< replaced by l_v |
---|
3515 | REAL(wp), PARAMETER :: food = 0._wp !< Heat gain by food (W) |
---|
3516 | REAL(wp), PARAMETER :: po = 1013.25_wp !< Air pressure at sea level (hPa) |
---|
3517 | REAL(wp), PARAMETER :: rob = 1.06_wp !< |
---|
3518 | ! |
---|
3519 | !-- Internal variables |
---|
3520 | REAL(wp) :: c(0:10) !< Core temperature array (degree_C) |
---|
3521 | REAL(wp) :: cbare !< Convection through bare skin |
---|
3522 | REAL(wp) :: cclo !< Convection through clothing |
---|
3523 | REAL(wp) :: csum !< Convection in total |
---|
3524 | REAL(wp) :: di !< difference between r1 and r2 |
---|
3525 | REAL(wp) :: ed !< energy transfer by diffusion (W) |
---|
3526 | REAL(wp) :: enbal !< energy ballance (W) |
---|
3527 | REAL(wp) :: enbal2 !< energy ballance (storage, last cycle) |
---|
3528 | REAL(wp) :: eswdif !< difference between sweat production and evaporation potential |
---|
3529 | REAL(wp) :: eswphy !< sweat created by physiology |
---|
3530 | REAL(wp) :: eswpot !< potential sweat evaporation |
---|
3531 | REAL(wp) :: fec !< |
---|
3532 | REAL(wp) :: hc !< |
---|
3533 | REAL(wp) :: he !< |
---|
3534 | REAL(wp) :: htcl !< |
---|
3535 | REAL(wp) :: r1 !< |
---|
3536 | REAL(wp) :: r2 !< |
---|
3537 | REAL(wp) :: rbare !< Radiational loss of bare skin (W/m²) |
---|
3538 | REAL(wp) :: rcl !< |
---|
3539 | REAL(wp) :: rclo !< Radiational loss of clothing (W/m²) |
---|
3540 | REAL(wp) :: rclo2 !< Longwave radiation gain or loss (W/m²) |
---|
3541 | REAL(wp) :: rsum !< Radiational loss or gain (W/m²) |
---|
3542 | REAL(wp) :: sw !< |
---|
3543 | ! REAL(wp) :: swf !< female factor, currently unused |
---|
3544 | REAL(wp) :: swm !< |
---|
3545 | REAL(wp) :: tbody !< |
---|
3546 | REAL(wp) :: tcore(1:7) !< |
---|
3547 | REAL(wp) :: vb !< |
---|
3548 | REAL(wp) :: vb1 !< |
---|
3549 | REAL(wp) :: vb2 !< |
---|
3550 | REAL(wp) :: wd !< |
---|
3551 | REAL(wp) :: wr !< |
---|
3552 | REAL(wp) :: ws !< |
---|
3553 | REAL(wp) :: wsum !< |
---|
3554 | REAL(wp) :: xx !< modification step (K) |
---|
3555 | REAL(wp) :: y !< fraction of bare skin |
---|
3556 | INTEGER(iwp) :: count1 !< running index |
---|
3557 | INTEGER(iwp) :: count3 !< running index |
---|
3558 | INTEGER(iwp) :: j !< running index |
---|
3559 | INTEGER(iwp) :: i !< running index |
---|
3560 | LOGICAL :: skipIncreaseCount !< iteration control flag |
---|
3561 | |
---|
3562 | wetsk = 0._wp |
---|
3563 | adu = 0.203_wp * mbody ** 0.425_wp * ht ** 0.725_wp |
---|
3564 | |
---|
3565 | hc = 2.67_wp + ( 6.5_wp * v ** 0.67_wp) |
---|
3566 | hc = hc * (pair /po) ** 0.55_wp |
---|
3567 | feff = 0.725_wp !< Posture: 0.725 for stading |
---|
3568 | facl = (- 2.36_wp + 173.51_wp * clo - 100.76_wp * clo * clo + 19.28_wp & |
---|
3569 | * (clo ** 3._wp)) / 100._wp |
---|
3570 | |
---|
3571 | IF ( facl > 1._wp ) facl = 1._wp |
---|
3572 | rcl = ( clo / 6.45_wp ) / facl |
---|
3573 | IF ( clo >= 2._wp ) y = 1._wp |
---|
3574 | |
---|
3575 | IF ( ( clo > 0.6_wp ) .AND. ( clo < 2._wp ) ) y = ( ht - 0.2_wp ) / ht |
---|
3576 | IF ( ( clo <= 0.6_wp ) .AND. ( clo > 0.3_wp ) ) y = 0.5_wp |
---|
3577 | IF ( ( clo <= 0.3_wp ) .AND. ( clo > 0._wp ) ) y = 0.1_wp |
---|
3578 | |
---|
3579 | r2 = adu * (fcl - 1._wp + facl) / (2._wp * 3.14_wp * ht * y) |
---|
3580 | r1 = facl * adu / (2._wp * 3.14_wp * ht * y) |
---|
3581 | |
---|
3582 | di = r2 - r1 |
---|
3583 | ! |
---|
3584 | !-- Skin temperatur |
---|
3585 | DO j = 1, 7 |
---|
3586 | |
---|
3587 | tsk = 34._wp |
---|
3588 | count1 = 0_iwp |
---|
3589 | tcl = ( ta + tmrt + tsk ) / 3._wp |
---|
3590 | count3 = 1_iwp |
---|
3591 | enbal2 = 0._wp |
---|
3592 | |
---|
3593 | DO i = 1, 100 ! allow for 100 iterations max |
---|
3594 | acl = adu * facl + adu * ( fcl - 1._wp ) |
---|
3595 | rclo2 = emcl * sigma_sb * ( ( tcl + degc_to_k )**4._wp - & |
---|
3596 | ( tmrt + degc_to_k )** 4._wp ) * feff |
---|
3597 | htcl = 6.28_wp * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl ) |
---|
3598 | tsk = 1._wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl |
---|
3599 | ! |
---|
3600 | !-- Radiation saldo |
---|
3601 | aeff = adu * feff |
---|
3602 | rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb * & |
---|
3603 | ( ( tmrt + degc_to_k )** 4._wp - ( tsk + degc_to_k )** 4._wp ) |
---|
3604 | rclo = feff * acl * emcl * sigma_sb * & |
---|
3605 | ( ( tmrt + degc_to_k )** 4._wp - ( tcl + degc_to_k )** 4._wp ) |
---|
3606 | rsum = rbare + rclo |
---|
3607 | ! |
---|
3608 | !-- Convection |
---|
3609 | cbare = hc * ( ta - tsk ) * adu * ( 1._wp - facl ) |
---|
3610 | cclo = hc * ( ta - tcl ) * acl |
---|
3611 | csum = cbare + cclo |
---|
3612 | ! |
---|
3613 | !-- Core temperature |
---|
3614 | c(0) = int_heat + ere |
---|
3615 | c(1) = adu * rob * cb |
---|
3616 | c(2) = 18._wp - 0.5_wp * tsk |
---|
3617 | c(3) = 5.28_wp * adu * c(2) |
---|
3618 | c(4) = 0.0208_wp * c(1) |
---|
3619 | c(5) = 0.76075_wp * c(1) |
---|
3620 | c(6) = c(3) - c(5) - tsk * c(4) |
---|
3621 | c(7) = - c(0) * c(2) - tsk * c(3) + tsk * c(5) |
---|
3622 | c(8) = c(6) * c(6) - 4._wp * c(4) * c(7) |
---|
3623 | c(9) = 5.28_wp * adu - c(5) - c(4) * tsk |
---|
3624 | c(10) = c(9) * c(9) - 4._wp * c(4) * & |
---|
3625 | ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk ) |
---|
3626 | |
---|
3627 | IF ( tsk == 36._wp ) tsk = 36.01_wp |
---|
3628 | tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600._wp ) + tsk |
---|
3629 | tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600._wp ) / & |
---|
3630 | ( 1._wp + 0.5_wp * ( 34._wp - tsk ) ) ) + tsk |
---|
3631 | IF ( c(10) >= 0._wp ) THEN |
---|
3632 | tcore(6) = ( - c(9) - c(10)**0.5_wp ) / ( 2._wp * c(4) ) |
---|
3633 | tcore(1) = ( - c(9) + c(10)**0.5_wp ) / ( 2._wp * c(4) ) |
---|
3634 | END IF |
---|
3635 | |
---|
3636 | IF ( c(8) >= 0._wp ) THEN |
---|
3637 | tcore(2) = ( - c(6) + ABS( c(8) ) ** 0.5_wp ) / ( 2._wp * c(4) ) |
---|
3638 | tcore(5) = ( - c(6) - ABS( c(8) ) ** 0.5_wp ) / ( 2._wp * c(4) ) |
---|
3639 | tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1._wp / 40._wp ) + tsk |
---|
3640 | END IF |
---|
3641 | ! |
---|
3642 | !-- Transpiration |
---|
3643 | tbody = 0.1_wp * tsk + 0.9_wp * tcore(j) |
---|
3644 | swm = 304.94_wp * ( tbody - 36.6_wp ) * adu / 3600000._wp |
---|
3645 | vpts = 6.11_wp * 10._wp**( 7.45_wp * tsk / ( 235._wp + tsk ) ) |
---|
3646 | |
---|
3647 | IF ( tbody <= 36.6_wp ) swm = 0._wp !< no need for sweating |
---|
3648 | |
---|
3649 | sw = swm |
---|
3650 | eswphy = - sw * l_v |
---|
3651 | he = 0.633_wp * hc / ( pair * c_p ) |
---|
3652 | fec = 1._wp / ( 1._wp + 0.92_wp * hc * rcl ) |
---|
3653 | eswpot = he * ( vpa - vpts ) * adu * l_v * fec |
---|
3654 | wetsk = eswphy / eswpot |
---|
3655 | |
---|
3656 | IF ( wetsk > 1._wp ) wetsk = 1._wp |
---|
3657 | ! |
---|
3658 | !-- Sweat production > evaporation? |
---|
3659 | eswdif = eswphy - eswpot |
---|
3660 | |
---|
3661 | IF ( eswdif <= 0._wp ) esw = eswpot !< Limit is evaporation |
---|
3662 | IF ( eswdif > 0._wp ) esw = eswphy !< Limit is sweat production |
---|
3663 | IF ( esw > 0._wp ) esw = 0._wp !< Sweat can't be evaporated, no more cooling effect |
---|
3664 | ! |
---|
3665 | !-- Diffusion |
---|
3666 | rdsk = 0.79_wp * 10._wp ** 7._wp |
---|
3667 | rdcl = 0._wp |
---|
3668 | ed = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( vpa - & |
---|
3669 | vpts ) |
---|
3670 | ! |
---|
3671 | !-- Max vb |
---|
3672 | vb1 = 34._wp - tsk |
---|
3673 | vb2 = tcore(j) - 36.6_wp |
---|
3674 | |
---|
3675 | IF ( vb2 < 0._wp ) vb2 = 0._wp |
---|
3676 | IF ( vb1 < 0._wp ) vb1 = 0._wp |
---|
3677 | vb = ( 6.3_wp + 75._wp * vb2 ) / ( 1._wp + 0.5_wp * vb1 ) |
---|
3678 | ! |
---|
3679 | !-- Energy ballence |
---|
3680 | enbal = int_heat + ed + ere + esw + csum + rsum + food |
---|
3681 | ! |
---|
3682 | !-- Clothing temperature |
---|
3683 | xx = 0.001_wp |
---|
3684 | IF ( count1 == 0_iwp ) xx = 1._wp |
---|
3685 | IF ( count1 == 1_iwp ) xx = 0.1_wp |
---|
3686 | IF ( count1 == 2_iwp ) xx = 0.01_wp |
---|
3687 | IF ( count1 == 3_iwp ) xx = 0.001_wp |
---|
3688 | |
---|
3689 | IF ( enbal > 0._wp ) tcl = tcl + xx |
---|
3690 | IF ( enbal < 0._wp ) tcl = tcl - xx |
---|
3691 | |
---|
3692 | skipIncreaseCount = .FALSE. |
---|
3693 | IF ( ( (enbal <= 0._wp ) .AND. (enbal2 > 0._wp ) ) .OR. & |
---|
3694 | ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) ) THEN |
---|
3695 | skipIncreaseCount = .TRUE. |
---|
3696 | ELSE |
---|
3697 | enbal2 = enbal |
---|
3698 | count3 = count3 + 1_iwp |
---|
3699 | END IF |
---|
3700 | |
---|
3701 | IF ( ( count3 > 200_iwp ) .OR. skipIncreaseCount ) THEN |
---|
3702 | IF ( count1 < 3_iwp ) THEN |
---|
3703 | count1 = count1 + 1_iwp |
---|
3704 | enbal2 = 0._wp |
---|
3705 | ELSE |
---|
3706 | EXIT |
---|
3707 | END IF |
---|
3708 | END IF |
---|
3709 | END DO |
---|
3710 | |
---|
3711 | IF ( count1 == 3_iwp ) THEN |
---|
3712 | SELECT CASE ( j ) |
---|
3713 | CASE ( 2, 5) |
---|
3714 | IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND. & |
---|
3715 | ( tsk <= 34.050_wp ) ) ) CYCLE |
---|
3716 | CASE ( 6, 1 ) |
---|
3717 | IF ( c(10) < 0._wp ) CYCLE |
---|
3718 | IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND. & |
---|
3719 | ( tsk > 33.850_wp ) ) ) CYCLE |
---|
3720 | CASE ( 3 ) |
---|
3721 | IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND. & |
---|
3722 | ( tsk <= 34.000_wp ) ) ) CYCLE |
---|
3723 | CASE ( 7 ) |
---|
3724 | IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND. & |
---|
3725 | ( tsk > 34.000_wp ) ) ) CYCLE |
---|
3726 | CASE default |
---|
3727 | END SELECT |
---|
3728 | END IF |
---|
3729 | |
---|
3730 | IF ( ( j /= 4_iwp ) .AND. ( vb >= 91._wp ) ) CYCLE |
---|
3731 | IF ( ( j == 4_iwp ) .AND. ( vb < 89._wp ) ) CYCLE |
---|
3732 | IF ( vb > 90._wp ) vb = 90._wp |
---|
3733 | ! |
---|
3734 | !-- Loses by water |
---|
3735 | ws = sw * 3600._wp * 1000._wp |
---|
3736 | IF ( ws > 2000._wp ) ws = 2000._wp |
---|
3737 | wd = ed / l_v * 3600._wp * ( -1000._wp ) |
---|
3738 | wr = erel / l_v * 3600._wp * ( -1000._wp ) |
---|
3739 | |
---|
3740 | wsum = ws + wr + wd |
---|
3741 | |
---|
3742 | RETURN |
---|
3743 | END DO |
---|
3744 | END SUBROUTINE heat_exch |
---|
3745 | |
---|
3746 | !------------------------------------------------------------------------------! |
---|
3747 | ! Description: |
---|
3748 | ! ------------ |
---|
3749 | !> Calculate PET |
---|
3750 | !------------------------------------------------------------------------------! |
---|
3751 | SUBROUTINE pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair, & |
---|
3752 | rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk ) |
---|
3753 | ! |
---|
3754 | !-- Input arguments: |
---|
3755 | REAL(wp), INTENT( IN ) :: acl !< clothing surface area (m²) |
---|
3756 | REAL(wp), INTENT( IN ) :: adu !< Du-Bois area (m²) |
---|
3757 | REAL(wp), INTENT( IN ) :: esw !< energy-loss through sweat evap. (W) |
---|
3758 | REAL(wp), INTENT( IN ) :: facl !< surface area extension through clothing (factor) |
---|
3759 | REAL(wp), INTENT( IN ) :: feff !< surface modification by posture (factor) |
---|
3760 | REAL(wp), INTENT( IN ) :: int_heat !< internal heat production (W) |
---|
3761 | REAL(wp), INTENT( IN ) :: pair !< air pressure (hPa) |
---|
3762 | REAL(wp), INTENT( IN ) :: rdcl !< diffusion resistence of clothing (factor) |
---|
3763 | REAL(wp), INTENT( IN ) :: rdsk !< diffusion resistence of skin (factor) |
---|
3764 | REAL(wp), INTENT( IN ) :: rtv !< respiratory volume |
---|
3765 | REAL(wp), INTENT( IN ) :: ta !< air temperature (degree_C) |
---|
3766 | REAL(wp), INTENT( IN ) :: tcl !< clothing temperature (degree_C) |
---|
3767 | REAL(wp), INTENT( IN ) :: tsk !< skin temperature (degree_C) |
---|
3768 | REAL(wp), INTENT( IN ) :: vpts !< sat. vapor pressure over skin (hPa) |
---|
3769 | REAL(wp), INTENT( IN ) :: wetsk !< fraction of wet skin (dimensionless) |
---|
3770 | ! |
---|
3771 | !-- Output arguments: |
---|
3772 | REAL(wp), INTENT( OUT ) :: aeff !< effective surface area (m²) |
---|
3773 | REAL(wp), INTENT( OUT ) :: pet_ij !< PET (degree_C) |
---|
3774 | ! |
---|
3775 | !-- Cconstants: |
---|
3776 | REAL(wp), PARAMETER :: emcl = 0.95_wp !< Longwave emission coef. of cloth |
---|
3777 | REAL(wp), PARAMETER :: emsk = 0.99_wp !< Longwave emission coef. of skin |
---|
3778 | REAL(wp), PARAMETER :: po = 1013.25_wp !< Air pressure at sea level (hPa) |
---|
3779 | ! |
---|
3780 | !-- Internal variables |
---|
3781 | REAL ( wp ) :: cbare !< Convection through bare skin |
---|
3782 | REAL ( wp ) :: cclo !< Convection through clothing |
---|
3783 | REAL ( wp ) :: csum !< Convection in total |
---|
3784 | REAL ( wp ) :: ed !< Diffusion (W) |
---|
3785 | REAL ( wp ) :: enbal !< Energy ballance (W) |
---|
3786 | REAL ( wp ) :: enbal2 !< Energy ballance (last iteration cycle) |
---|
3787 | REAL ( wp ) :: ere !< Energy ballance result (W) |
---|
3788 | REAL ( wp ) :: erel !< Latent energy ballance (W) |
---|
3789 | REAL ( wp ) :: eres !< Sensible respiratory heat flux (W) |
---|
3790 | REAL ( wp ) :: hc !< |
---|
3791 | REAL ( wp ) :: rbare !< Radiational loss of bare skin (W/m²) |
---|
3792 | REAL ( wp ) :: rclo !< Radiational loss of clothing (W/m²) |
---|
3793 | REAL ( wp ) :: rsum !< Radiational loss or gain (W/m²) |
---|
3794 | REAL ( wp ) :: tex !< Temperat. of exhaled air (degree_C) |
---|
3795 | REAL ( wp ) :: vpex !< Vapor pressure of exhaled air (hPa) |
---|
3796 | REAL ( wp ) :: xx !< Delta PET per iteration (K) |
---|
3797 | |
---|
3798 | INTEGER ( iwp ) :: count1 !< running index |
---|
3799 | INTEGER ( iwp ) :: i !< running index |
---|
3800 | |
---|
3801 | pet_ij = ta |
---|
3802 | enbal2 = 0._wp |
---|
3803 | |
---|
3804 | DO count1 = 0, 3 |
---|
3805 | DO i = 1, 125 ! 500 / 4 |
---|
3806 | hc = 2.67_wp + 6.5_wp * 0.1_wp ** 0.67_wp |
---|
3807 | hc = hc * ( pair / po ) ** 0.55_wp |
---|
3808 | ! |
---|
3809 | !-- Radiation |
---|
3810 | aeff = adu * feff |
---|
3811 | rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb * & |
---|
3812 | ( ( pet_ij + degc_to_k ) ** 4._wp - ( tsk + degc_to_k ) ** 4._wp ) |
---|
3813 | rclo = feff * acl * emcl * sigma_sb * & |
---|
3814 | ( ( pet_ij + degc_to_k ) ** 4._wp - ( tcl + degc_to_k ) ** 4._wp ) |
---|
3815 | rsum = rbare + rclo |
---|
3816 | ! |
---|
3817 | !-- Covection |
---|
3818 | cbare = hc * ( pet_ij - tsk ) * adu * ( 1._wp - facl ) |
---|
3819 | cclo = hc * ( pet_ij - tcl ) * acl |
---|
3820 | csum = cbare + cclo |
---|
3821 | ! |
---|
3822 | !-- Diffusion |
---|
3823 | ed = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( 12._wp - & |
---|
3824 | vpts ) |
---|
3825 | ! |
---|
3826 | !-- Respiration |
---|
3827 | tex = 0.47_wp * pet_ij + 21._wp |
---|
3828 | eres = c_p * ( pet_ij - tex ) * rtv |
---|
3829 | vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) ) |
---|
3830 | erel = 0.623_wp * l_v / pair * ( 12._wp - vpex ) * rtv |
---|
3831 | ere = eres + erel |
---|
3832 | ! |
---|
3833 | !-- Energy ballance |
---|
3834 | enbal = int_heat + ed + ere + esw + csum + rsum |
---|
3835 | ! |
---|
3836 | !-- Iteration concerning ta |
---|
3837 | xx = 0.001_wp |
---|
3838 | IF ( count1 == 0_iwp ) xx = 1._wp |
---|
3839 | IF ( count1 == 1_iwp ) xx = 0.1_wp |
---|
3840 | IF ( count1 == 2_iwp ) xx = 0.01_wp |
---|
3841 | ! IF ( count1 == 3_iwp ) xx = 0.001_wp |
---|
3842 | IF ( enbal > 0._wp ) pet_ij = pet_ij - xx |
---|
3843 | IF ( enbal < 0._wp ) pet_ij = pet_ij + xx |
---|
3844 | IF ( ( enbal <= 0._wp ) .AND. ( enbal2 > 0._wp ) ) EXIT |
---|
3845 | IF ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) EXIT |
---|
3846 | |
---|
3847 | enbal2 = enbal |
---|
3848 | END DO |
---|
3849 | END DO |
---|
3850 | END SUBROUTINE pet_iteration |
---|
3851 | |
---|
3852 | ! |
---|
3853 | !-- UVEM specific subroutines |
---|
3854 | |
---|
3855 | !---------------------------------------------------------------------------------------------------------------------! |
---|
3856 | ! Description: |
---|
3857 | ! ------------ |
---|
3858 | !> Module-specific routine for new module |
---|
3859 | !---------------------------------------------------------------------------------------------------------------------! |
---|
3860 | SUBROUTINE uvem_solar_position |
---|
3861 | |
---|
3862 | USE date_and_time_mod, & |
---|
3863 | ONLY: calc_date_and_time, day_of_year, time_utc |
---|
3864 | |
---|
3865 | USE control_parameters, & |
---|
3866 | ONLY: latitude, longitude |
---|
3867 | |
---|
3868 | IMPLICIT NONE |
---|
3869 | |
---|
3870 | |
---|
3871 | REAL(wp) :: alpha = 0.0_wp !< solar azimuth angle in radiant |
---|
3872 | REAL(wp) :: doy_r = 0.0_wp !< real format of day_of_year |
---|
3873 | REAL(wp) :: declination = 0.0_wp !< declination |
---|
3874 | REAL(wp) :: dtor = 0.0_wp !< factor to convert degree to radiant |
---|
3875 | REAL(wp) :: js = 0.0_wp !< parameter for solar position calculation |
---|
3876 | REAL(wp) :: lat = 52.39_wp !< latitude |
---|
3877 | REAL(wp) :: lon = 9.7_wp !< longitude |
---|
3878 | REAL(wp) :: thetar = 0.0_wp !< angle for solar zenith angle calculation |
---|
3879 | REAL(wp) :: thetasr = 0.0_wp !< angle for solar azimuth angle calculation |
---|
3880 | REAL(wp) :: zgl = 0.0_wp !< calculated exposure by direct beam |
---|
3881 | REAL(wp) :: woz = 0.0_wp !< calculated exposure by diffuse radiation |
---|
3882 | REAL(wp) :: wsp = 0.0_wp !< calculated exposure by direct beam |
---|
3883 | |
---|
3884 | |
---|
3885 | CALL calc_date_and_time |
---|
3886 | doy_r = real(day_of_year) |
---|
3887 | dtor = pi / 180.0_wp |
---|
3888 | lat = latitude |
---|
3889 | lon = longitude |
---|
3890 | ! |
---|
3891 | !-- calculation of js, necessary for calculation of equation of time (zgl) : |
---|
3892 | js= 72.0_wp * ( doy_r + ( time_utc / 86400.0_wp ) ) / 73.0_wp |
---|
3893 | ! |
---|
3894 | !-- calculation of equation of time (zgl): |
---|
3895 | zgl = 0.0066_wp + 7.3525_wp * cos( ( js + 85.9_wp ) * dtor ) + 9.9359_wp * & |
---|
3896 | cos( ( 2.0_wp * js + 108.9_wp ) * dtor ) + 0.3387_wp * cos( ( 3 * js + 105.2_wp ) * dtor ) |
---|
3897 | ! |
---|
3898 | !-- calculation of apparent solar time woz: |
---|
3899 | woz = ( ( time_utc / 3600.0_wp ) - ( 4.0_wp * ( 15.0_wp - lon ) ) / 60.0_wp ) + ( zgl / 60.0_wp ) |
---|
3900 | ! |
---|
3901 | !-- calculation of hour angle (wsp): |
---|
3902 | wsp = ( woz - 12.0_wp ) * 15.0_wp |
---|
3903 | ! |
---|
3904 | !-- calculation of declination: |
---|
3905 | declination = 0.3948_wp - 23.2559_wp * cos( ( js + 9.1_wp ) * dtor ) - & |
---|
3906 | 0.3915_wp * cos( ( 2.0_wp * js + 5.4_wp ) * dtor ) - 0.1764_wp * cos( ( 3.0_wp * js + 26.0_wp ) * dtor ) |
---|
3907 | ! |
---|
3908 | !-- calculation of solar zenith angle |
---|
3909 | thetar = acos( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) * & |
---|
3910 | cos( lat * dtor ) * cos( declination * dtor ) ) |
---|
3911 | thetasr = asin( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) * & |
---|
3912 | cos( lat * dtor ) * cos( declination * dtor ) ) |
---|
3913 | sza = thetar / dtor |
---|
3914 | ! |
---|
3915 | !-- calculation of solar azimuth angle |
---|
3916 | IF (woz .LE. 12.0_wp) alpha = pi - acos( ( sin(thetasr) * sin( lat * dtor ) - & |
---|
3917 | sin( declination * dtor ) ) / ( cos(thetasr) * cos( lat * dtor ) ) ) |
---|
3918 | IF (woz .GT. 12.0_wp) alpha = pi + acos( ( sin(thetasr) * sin( lat * dtor ) - & |
---|
3919 | sin( declination * dtor ) ) / ( cos(thetasr) * cos( lat * dtor ) ) ) |
---|
3920 | saa = alpha / dtor |
---|
3921 | |
---|
3922 | END SUBROUTINE uvem_solar_position |
---|
3923 | |
---|
3924 | |
---|
3925 | !------------------------------------------------------------------------------! |
---|
3926 | ! Description: |
---|
3927 | ! ------------ |
---|
3928 | !> Module-specific routine for new module |
---|
3929 | !---------------------------------------------------------------------------------------------------------------------! |
---|
3930 | SUBROUTINE uvem_calc_exposure |
---|
3931 | |
---|
3932 | USE indices, & |
---|
3933 | ONLY: nys, nyn, nxl, nxr |
---|
3934 | |
---|
3935 | |
---|
3936 | IMPLICIT NONE |
---|
3937 | |
---|
3938 | INTEGER(iwp) :: i !< loop index in x direction |
---|
3939 | INTEGER(iwp) :: j !< loop index in y direction |
---|
3940 | INTEGER(iwp) :: szai !< loop index for different sza values |
---|
3941 | |
---|
3942 | CALL uvem_solar_position |
---|
3943 | |
---|
3944 | IF (sza .GE. 90) THEN |
---|
3945 | vitd3_exposure(:,:) = 0.0_wp |
---|
3946 | ELSE |
---|
3947 | |
---|
3948 | DO ai = 0, 35 |
---|
3949 | DO zi = 0, 9 |
---|
3950 | projection_area_lookup_table(ai,zi) = uvem_projarea_f%var(clothing,zi,ai) |
---|
3951 | ENDDO |
---|
3952 | ENDDO |
---|
3953 | DO ai = 0, 35 |
---|
3954 | DO zi = 0, 9 |
---|
3955 | integration_array(ai,zi) = uvem_integration_f%var(zi,ai) |
---|
3956 | ENDDO |
---|
3957 | ENDDO |
---|
3958 | DO ai = 0, 2 |
---|
3959 | DO zi = 0, 90 |
---|
3960 | irradiance_lookup_table(ai,zi) = uvem_irradiance_f%var(zi,ai) |
---|
3961 | ENDDO |
---|
3962 | ENDDO |
---|
3963 | DO ai = 0, 35 |
---|
3964 | DO zi = 0, 9 |
---|
3965 | DO szai = 0, 90 |
---|
3966 | radiance_lookup_table(ai,zi,szai) = uvem_radiance_f%var(szai,zi,ai) |
---|
3967 | ENDDO |
---|
3968 | ENDDO |
---|
3969 | ENDDO |
---|
3970 | |
---|
3971 | |
---|
3972 | |
---|
3973 | !-- rotate 3D-Model human to desired direction ----------------------------- |
---|
3974 | projection_area_temp( 0:35,:) = projection_area_lookup_table |
---|
3975 | projection_area_temp(36:71,:) = projection_area_lookup_table |
---|
3976 | IF ( .NOT. turn_to_sun ) startpos_human = orientation_angle / 10.0_wp |
---|
3977 | IF ( turn_to_sun ) startpos_human = saa / 10.0_wp |
---|
3978 | DO ai = 0, 35 |
---|
3979 | xfactor = ( startpos_human ) - INT( startpos_human ) |
---|
3980 | DO zi = 0, 9 |
---|
3981 | projection_area(ai,zi) = ( projection_area_temp( 36 - INT( startpos_human ) - 1 + ai , zi) * & |
---|
3982 | ( xfactor ) ) & |
---|
3983 | +( projection_area_temp( 36 - INT( startpos_human ) + ai , zi) * & |
---|
3984 | ( 1.0_wp - xfactor ) ) |
---|
3985 | ENDDO |
---|
3986 | ENDDO |
---|
3987 | ! |
---|
3988 | ! |
---|
3989 | !-- interpolate to accurate Solar Zenith Angle ------------------ |
---|
3990 | DO ai = 0, 35 |
---|
3991 | xfactor = (sza)-INT(sza) |
---|
3992 | DO zi = 0, 9 |
---|
3993 | radiance_array(ai,zi) = ( radiance_lookup_table(ai, zi, INT(sza) ) * ( 1.0_wp - xfactor) ) + & |
---|
3994 | ( radiance_lookup_table(ai,zi,INT(sza) + 1) * xfactor ) |
---|
3995 | ENDDO |
---|
3996 | ENDDO |
---|
3997 | DO iq = 0, 2 |
---|
3998 | irradiance(iq) = ( irradiance_lookup_table(iq, INT(sza) ) * ( 1.0_wp - xfactor)) + & |
---|
3999 | (irradiance_lookup_table(iq, INT(sza) + 1) * xfactor ) |
---|
4000 | ENDDO |
---|
4001 | ! |
---|
4002 | !-- interpolate to accurate Solar Azimuth Angle ------------------ |
---|
4003 | IF ( sun_in_south ) THEN |
---|
4004 | startpos_saa_float = 180.0_wp / 10.0_wp |
---|
4005 | ELSE |
---|
4006 | startpos_saa_float = saa / 10.0_wp |
---|
4007 | ENDIF |
---|
4008 | radiance_array_temp( 0:35,:) = radiance_array |
---|
4009 | radiance_array_temp(36:71,:) = radiance_array |
---|
4010 | xfactor = (startpos_saa_float) - INT(startpos_saa_float) |
---|
4011 | DO ai = 0, 35 |
---|
4012 | DO zi = 0, 9 |
---|
4013 | radiance_array(ai,zi) = ( radiance_array_temp( 36 - INT( startpos_saa_float ) - 1 + ai , zi ) * & |
---|
4014 | ( xfactor ) ) & |
---|
4015 | + ( radiance_array_temp( 36 - INT( startpos_saa_float ) + ai , zi ) & |
---|
4016 | * ( 1.0_wp - xfactor ) ) |
---|
4017 | ENDDO |
---|
4018 | ENDDO |
---|
4019 | ! |
---|
4020 | ! |
---|
4021 | !-- calculate Projectionarea for direct beam -----------------------------' |
---|
4022 | projection_area_direct_temp( 0:35,:) = projection_area |
---|
4023 | projection_area_direct_temp(36:71,:) = projection_area |
---|
4024 | yfactor = ( sza / 10.0_wp ) - INT( sza / 10.0_wp ) |
---|
4025 | xfactor = ( startpos_saa_float ) - INT( startpos_saa_float ) |
---|
4026 | projection_area_direct_beam = ( projection_area_direct_temp( INT(startpos_saa_float) ,INT(sza/10.0_wp) ) * & |
---|
4027 | ( 1.0_wp - xfactor ) * ( 1.0_wp - yfactor ) ) + & |
---|
4028 | ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp) ) * & |
---|
4029 | ( xfactor ) * ( 1.0_wp - yfactor ) ) + & |
---|
4030 | ( projection_area_direct_temp( INT(startpos_saa_float) ,INT(sza/10.0_wp)+1) * & |
---|
4031 | ( 1.0_wp - xfactor ) * ( yfactor ) ) + & |
---|
4032 | ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)+1) * & |
---|
4033 | ( xfactor ) * ( yfactor ) ) |
---|
4034 | ! |
---|
4035 | ! |
---|
4036 | ! |
---|
4037 | DO i = nxl, nxr |
---|
4038 | DO j = nys, nyn |
---|
4039 | ! |
---|
4040 | ! !-- extract obstruction from IBSET-Integer_Array ------------------' |
---|
4041 | IF (consider_obstructions ) THEN |
---|
4042 | obstruction_temp1 = building_obstruction_f%var_3d(:,j,i) |
---|
4043 | IF (obstruction_temp1(0) .NE. 9) THEN |
---|
4044 | DO pobi = 0, 44 |
---|
4045 | DO bi = 0, 7 |
---|
4046 | IF ( btest( obstruction_temp1(pobi), bi ) .EQV. .TRUE.) THEN |
---|
4047 | obstruction_temp2( ( pobi * 8 ) + bi ) = 1 |
---|
4048 | ELSE |
---|
4049 | obstruction_temp2( ( pobi * 8 ) + bi ) = 0 |
---|
4050 | ENDIF |
---|
4051 | ENDDO |
---|
4052 | ENDDO |
---|
4053 | DO zi = 0, 9 |
---|
4054 | obstruction(:,zi) = obstruction_temp2( zi * 36 :( zi * 36) + 35 ) |
---|
4055 | ENDDO |
---|
4056 | ELSE |
---|
4057 | obstruction(:,:) = 0 |
---|
4058 | ENDIF |
---|
4059 | ENDIF |
---|
4060 | ! |
---|
4061 | ! !-- calculated human exposure ------------------' |
---|
4062 | diffuse_exposure = SUM( radiance_array * projection_area * integration_array * obstruction ) |
---|
4063 | |
---|
4064 | obstruction_direct_beam = obstruction( nint(startpos_saa_float), nint( sza / 10.0_wp ) ) |
---|
4065 | IF (sza .GE. 89.99_wp) THEN |
---|
4066 | sza = 89.99999_wp |
---|
4067 | ENDIF |
---|
4068 | ! |
---|
4069 | !-- calculate direct normal irradiance (direct beam) ------------------' |
---|
4070 | direct_exposure = ( irradiance(1) / cos( pi * sza / 180.0_wp ) ) * & |
---|
4071 | projection_area_direct_beam * obstruction_direct_beam |
---|
4072 | |
---|
4073 | vitd3_exposure(j,i) = ( diffuse_exposure + direct_exposure ) / 1000.0_wp * 70.97_wp |
---|
4074 | ! unit = international units vitamin D per second |
---|
4075 | ENDDO |
---|
4076 | ENDDO |
---|
4077 | ENDIF |
---|
4078 | |
---|
4079 | END SUBROUTINE uvem_calc_exposure |
---|
4080 | |
---|
4081 | END MODULE biometeorology_mod |
---|