Home

Context Navigation

biometeorology_mod.f90 @ 4177

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

Replace get_topography_top_index functions by pre-calculated arrays in order to save computational resources

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

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

Context Navigation

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

Download in other formats: