Home

Context Navigation

biometeorology_mod.f90 @ 3735

Last change on this file since 3735 was 3735, checked in by dom_dwd_user, 5 years ago

biometeorology_mod.f90:
(N) Fixed auto-setting of thermal index calculation flags by output as
originally proposed by resler.
(C) removed bio_pet and outher configuration variables.
(C) Updated namelist.
(B) Forcing initialization of tmrt_av_grid to avoid mysterious mrt
values at i==0, j==0

module_interface_mod.f90:
(C) Receiving parameter j (averaging 0==.F./1==.T.) in
module_interface_check_data_output from check_parameters.f90.
(C) Passing j to bio_check_parameters.

check_parameters.f90:
(C) Passing j to module_interface_check_data_output

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

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

Context Navigation

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

Download in other formats: