Home

Context Navigation

biometeorology_mod.f90 @ 4583

Last change on this file since 4583 was 4577, checked in by raasch, 4 years ago

further re-formatting to follow the PALM coding standard

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

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

Context Navigation

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

Download in other formats: