Home

Context Navigation

biometeorology_mod.f90 @ 3749

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

biometeorology_mod.f90:
(B) Added addittional safety meassures to bio_calculate_thermal_index_maps.
(B) Replaced several REAL (un-)equality comparisons.

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

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

Context Navigation

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

Download in other formats: