Home

Context Navigation

biometeorology_mod.f90 @ 3609

Last change on this file since 3609 was 3593, checked in by kanani, 5 years ago

Bugfix for missing array allocation (biometeorology_mod), remove degree symbol (biometeorology_mod, indoor_model_mod, multi_agent_system_mod, surface_mod, wind_turbine_model_mod)

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

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

Context Navigation

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

Download in other formats: