Home

Context Navigation

biometeorology_mod.f90 @ 3582

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

Merge branch salsa with trunk

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

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

Context Navigation

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

Download in other formats: