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