Home

Context Navigation

biometeorology_mod.f90 @ 3712

Last change on this file since 3712 was 3711, checked in by knoop, 6 years ago

Introduced module_interface_init_checks for post-init checks

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

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

Context Navigation

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

Download in other formats: