Home

Context Navigation

biometeorology_mod.f90 @ 3550

Last change on this file since 3550 was 3525, checked in by kanani, 6 years ago

Changes related to clean-up of biometeorology (by dom_dwd_user)

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

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

Context Navigation

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

Download in other formats: