Home

Context Navigation

biometeorology_mod.f90 @ 4858

Last change on this file since 4858 was 4843, checked in by raasch, 4 years ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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-4492
/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: 198.1 KB

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

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

Context Navigation

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

Download in other formats: