Home

Context Navigation

biometeorology_mod.f90 @ 3739

Last change on this file since 3739 was 3739, checked in by dom_dwd_user, 6 years ago

biometeorology_mod.f90:
(N) Auto-adjusting thermal_comfort flag if not set by user, but thermal_indices set as output quantities.
(C) Renamed flags "bio_<index>" to "do_calculate_<index>" for better readability
(C) Removed everything related to "time_bio_results" as this is never used.
(C) Moved humidity warning to check_data_output so it will also be triggered if thermal_comofrt flag was auto-set later.
(B) Fixed bug in mrt calculation introduced with my commit yesterday.

time_integration.f90:
(C) Removed everything related to "time_bio_results" as this is never used

module_interface.f90:
(C) Removed call to empty method bio_check_parameters.

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/chemistry/SOURCE/biometeorology_mod.f90	2047-3190,3218-3297
/palm/branches/resler/SOURCE/biometeorology_mod.f90	2023-3320,3337-3474
/palm/branches/salsa/SOURCE/biometeorology_mod.f90	2503-3581
/palm/trunk/SOURCE/biometeorology_mod.f90	merged	eligible
/palm/branches/forwind/SOURCE/biometeorology_mod.f90	1564-1913
/palm/branches/fricke/SOURCE/biometeorology_mod.f90	942-977
/palm/branches/hoffmann/SOURCE/biometeorology_mod.f90	989-1052
/palm/branches/letzel/masked_output/SOURCE/biometeorology_mod.f90	296-409
/palm/branches/palm4u/SOURCE/biometeorology_mod.f90	2540-2692
/palm/branches/rans/SOURCE/biometeorology_mod.f90	2078-3128
/palm/branches/suehring/biometeorology_mod.f90	423-666

File size: 177.1 KB

Line
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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: palm/trunk/SOURCE/biometeorology_mod.f90 @ 3739

Download in other formats: