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