Home

Context Navigation

biometeorology_mod.f90 @ 4397

Last change on this file since 4397 was 4392, checked in by pavelkrc, 5 years ago

Merge branch resler (updated documentation, optional flux tracing, bugfix)

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-4391
/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: 181.3 KB

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

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

Context Navigation

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

Download in other formats: