Home

Context Navigation

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

Last change on this file since 4830 was 4828, checked in by Giersch, 4 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |