Home

Context Navigation

biometeorology_mod.f90 @ 4495

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

restart data handling with MPI-IO added, first part

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/chemistry/SOURCE/biometeorology_mod.f90	2047-3190,3218-3297
/palm/branches/resler/SOURCE/biometeorology_mod.f90	2023-4492
/palm/branches/salsa/SOURCE/biometeorology_mod.f90	2503-3581
/palm/trunk/SOURCE/biometeorology_mod.f90	merged	eligible
/palm/branches/forwind/SOURCE/biometeorology_mod.f90	1564-1913
/palm/branches/fricke/SOURCE/biometeorology_mod.f90	942-977
/palm/branches/hoffmann/SOURCE/biometeorology_mod.f90	989-1052
/palm/branches/letzel/masked_output/SOURCE/biometeorology_mod.f90	296-409
/palm/branches/palm4u/SOURCE/biometeorology_mod.f90	2540-2692
/palm/branches/rans/SOURCE/biometeorology_mod.f90	2078-3128
/palm/branches/suehring/biometeorology_mod.f90	423-666

File size: 184.1 KB

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

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

Context Navigation

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

Download in other formats: