Home

Context Navigation

biometeorology_mod.f90 @ 4777

Last change on this file since 4777 was 4768, checked in by suehring, 4 years ago

Enable 3D data output also with 64-bit precision

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

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

Context Navigation

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

Download in other formats: