Home

Context Navigation

biometeorology_mod.f90 @ 4127

Last change on this file since 4127 was 4127, checked in by suehring, 5 years ago

Merge with branch resler: biomet- output of bio_mrt added; plant_canopy - separate vertical dimension for 3D output (to save disk space); radiation - remove unused plant canopy variables; urban-surface model - do not add anthropogenic heat during wall spin-up

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-3320,3337-3474
/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: 182.5 KB

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

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

Context Navigation

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

Download in other formats: