Home

Context Navigation

biometeorology_mod.f90 @ 3924

Last change on this file since 3924 was 3885, checked in by kanani, 6 years ago

restructure/add location/debug messages

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

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

Context Navigation

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

Download in other formats: