Home

Context Navigation

biometeorology_mod.f90 @ 3685

Last change on this file since 3685 was 3685, checked in by knoop, 5 years ago

Some interface calls moved to module_interface + cleanup

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

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

Context Navigation

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

Download in other formats: